Towards Resistant Geostatistics - Cressie2
Towards Resistant Geostatistics - Cressie2
Towards Resistant Geostatistics - Cressie2
Resources Characterization
Part 1
NATO ASI Series
Advanced Science Institutes Series
A series presenting the results of activities sponsored by the NATO Science Committee,
which aims at the dissemination of advanced scientific and technological knowledge,
with a view to strengthening links between scientific communities.
edited by
G. Verly
Department of Applied Earth Sciences, Stanford University, Stanford, U.S.A.
M. David
Departement de Chimie Mineral , Ecole Poly technique, Montreal, Canada
A. G. Journel
Department of Applied Earth Sciences, Stanford University, Stanford, U.S.A.
and
A. Marechal
S.N. Elf Aquitaine (Production), Pau, France
No part of the material protected by this copyright notice may be reproduced or utilized
in any form or by any means, electronic or mechanical, including photocopying, recording
or by any information storage and retrieval system, without written permission from the
copyright owner.
TABLE OF CONTENTS
FOREWORD xiii
ACKNOWLEDGEMENTS xiv
LIST OF PARTICIPANTS xv
PART I - VARIOGRAM
LECTURES
14 - POSITIVE KRIGING
R.J. Barnes. T.B. Johnson 231
15 - CORRECTING CONDITIONAL BIAS
K. Guertin 245
LECTURE
LECTURES
23 - RECOVERY ESTIMATION: A REVIEW OF MODELS AND
METHODS
A. Marechal 385
24 - THE SELECTIVITY OF THE DISTRIBUTIONS AND THE
"SECOND PRINCIPLE OF GEOSTATISTICS"
G. Matheron 421
CONTR I BUTI ONS
LECTURE
Since October 1975 and the first NATO ASl "Geostat - Roma
1975", there has not been an advanced workshop on geostatistics
where experts from throughout the world could meet, teach and
debate without the pressure of an ordinary professional confer-
ence. This second NATO ASI "Geostat - Tahoe 1983" was intended
as a high-level teaching activity yet opened to all new ideas and
contributions to the discipline of geostatistics. It was
expected that the institute would bridge the gap since "Geostat -
Roma 1975" and establish the state of the art of the discipline
as of 1983.
"Geostat - Tahoe 1983" fulfilled all expectations. The
institute, held in September 1983 at the Stanford Sierra Lodge
near Lake Tahoe, California, was attended by all major players in
the field, representing evenly the academy and the industry of 23
different countries. Twelve guest lecturers were backed by some
50 original contributions. Particularly important for the future
was the active participation of graduate students, glvlng evi-
dence of the dynamism of a still very young discipline.
A.G. Journel
Director of the ASI
ACKNO\{LEDGEMENTS
BRGM, France
CERCHAR, France
COGEMA, France
Elf Aquitaine (Production), France
Schlumberger E.P.S., France
TECMIN, France
Shell, Netherlands
George Verly
Michel David
Andre Journel
Alain Marechal
xi
LIST OF PARTICIPANTS
BAElE steve M.: Mining Engineer, u.s. Borax Corp., 3075 Wishire
Blvd., Los Angeles, California 90010, U.S.A.
BRYAN, Rex, C.: Director, Geostat Systems Inc., P.O. Box 1193,
Golden, Colorado 80402. U.S.A.
CHU, Darong: Student, 1126 C San Pablo Av., "3, Albany, Cali-
fornia 94706. U.S.A.
CHUNG, Chang-Jo r.: Research Scientist. Geological Survey of
Canada. 601 Booth Street. Ottawa, Ontario K1A OE8, Canada
RAYMOND. Gary F.: Sr. Mine Project Engineer. Cominco. Ltd .• P.O.
Box 2000. Kimberley. British Colombia V1A 2G3. Canada
YANCEY. James D.: Sr. Research Associate. Arco Oil & Gas Co .•
P.O. Box 2819. Dallas. Texas 75221. U.S.A.
Margaret ARMSTRONG
Centre de Geostatistique et Morphologie Mathematique
Ecole Nationale Superieure des Mines de Paris
35 rue Saint Honore, 77305 rONTAINEBLEAU (France)
1. INTRODUCTION
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1,1-19.
© 1984 by D. Reidel Publishing Company.
2 M.ARMSTRONG
~estion. Improvements have been made here and there but an over~
all approach has been lacking. It is to be hoped that the papers
presented in this section of the workshoD and the ensuing dis-
cussion will lead to a more unified view of the problems. With
this in mind. the objective of this paper is to provide a schematic
overview of the question. to indicate the main work that has al-
ready been done and to suggest points where improvements are
needed. Some specific proposals will be made.
~
:0
5
Cl
:0
TABLE 1
~
!J>
6 M. ARMSTRONG
Going further one could ask to what extent the kriging weights
and the kriging variance depend on these two factor~. Most of the
time, the results of kriging are fairly stable. However since
thi~ is not always so, it would be better to be able to quantify
the robustness of kriging to small perturbations. This would help
IMPROVING THE ESTIMATION AND MODELLING OF THE VARIOGRAM 9
EXAMPLE 2 : Kriging the same 1.0 x 1.0 square using the same
Gaussian model (with no nugget effectJ using three sample points.
Two sets of data locations considered.
y
r X = B (2)
Y
where X T (A 1 , .. ., AN,]J)
Y
B T Cy (x 1, V), ... , y(xN,V) , 1 )
Y
•• I •••• .; o
Similarly the system corresponding to g can be written as
GX =B •
g g
From (1),. integration over V gives, for i = 1, 2, ••• ,N,
whence
where nBll is some Up) vector norm (usually with p 1,2 or (0).
IMPROVING THE ESTIMATION AND MODELLING OF THE VARIOGRAM 11
On letting /l,X xy xg AB B
'(
B
g
and (II' f - G. it
can be shown that
o ~ Iti ,
1
(y,
1
- [' + [ ex p (- h , a ' ) J 2
1
( 7)
The fact that the parameter a (or a') appears in the expo-
nential makes these difficult to solve. A similar but even nastier
problem ari~es with the spherical model. In that case. the contri-
bution of any value to the sum of s~uares depends on whether the
distance h, is greater than the range or not. Consequently in
both these1cases it is possible to minimize the sum of squares
in a two stage procedure By solving the equations for a fixed
value of "a" and then letting "a" vary. but this a very unwieldy
procedure. If two or more nested structures are used one would
have great difficulty finding the minimum.
S Y(hJ se-shdh
o
with the observed quantity. In other words. the theoretical La-
place transform is equated with the observed one. From the geo-
statistical point of view the weighting function se- sh is ideal
as it downweights the tail of the variogram.
IMPROVING THE ESTIMATION AND MODELLING OF THE VARIOGRAM 15
Hence
L2 s2 - L1 S1
a Ll - L2
(12 )
L1 L2 (s 1 - S2)
C
Li' s 1 - L2 s2
y[h) = C ( -3h -h 3 )
(13 )
2a 2a3
is
,§£ [s) = r: ~ l
~ 2 as
- ~2(as) 3 + 3 e
[as)
~~2 (1 + ...1-)
as ~
t (14 )
Geost atistic ians are v'ell aware of the crucia l role of the
proce-
variogr am model in the overal l geosta tistica l estima tion
6 neighb ourhoo d around a variogr am merely
dure. The concep t of a
theore tical framew ork for the idea that simila r-look ing
provid es a
s from kriging It also
Variog rams ought to lead to simila r result
weight s
gives specif ic bounds to the relativ e change in the kriging
ce resulti ng from a slight pertur bation in
and the kriging varian
the variogr am model or the locatio n of the data points .
either
a
Implic itly it indica tes now fussy we need to be when fitting s
is Known up to a 6 neighb ourhoo d, the result
model. Once the model
r this impli-
of the kriging are not going to change much. Howeve
choice )
citly requir es a good knowledge [or rather , a suitab le
, partic ularly if the latter is small rela-
of the nugget effect
total sill. In the long run, these ideas also have
tive to the
on propos ed for fit-
implic ations for any goodne ss of fit criteri
want any such criteri on to be
ting variogr am models . One would
the sense that tWD differ ent models belong ing to
consis tent in
fit"
the same a neighb ourhoo d ought to have simila r "goodn ess of
e compar able result s in the
statis tics since they would produc
krigin g.
mi-
This notion of the sensit ivity of the kriging system to
one aspect of the wider questio n of
nor pertur bation s is just
taken
"robus tness" . Unfort unately most of the work to date has
have been made to
a narrow er vision of the problem . Sugges tions
the geosta tistica l estima tion
"robus tify" partic ular steps in
ures,
proced ure, genera lly by inserti ng robust statis tical proced
but an overal l approa ch to the questio n has been lackin g.
2) What do we need ?
[i) a robust estima te of the experim ental variogr am
REFERENCES
Noel Cressie
ABSTRACT
Implicit in many of the geostatistical techniques developed, is
faith that the data are Gaussian (normal) or can be conveniently
transformed to Gaussianity. A mining engineer knows however
that contamination is present all too often. This paper will
make an exploratory analysis of spatial data, graphing and sum-
marlzlng in a way that is resistant to that contamination. The
ultimate goal is variogram estimation and kriging, and the data
analytic techniques used reflect this. The pocket plot is a new
way of looking at contributions of small regions of points to
variogram estimation. The gridded data are thought of as a higher
way table and analysed by median polish. The residual table is
shown to contain useful information on the spatial relationships
between data. Coal ash measurements in Pennsylvania are analysed
from this resistant point of view.
INTRODUCTION
For some years now we have held the view that geostatistics
could profit considerably from the philosophy and methods of mod-
ern data analysis, as it currently pervades the mainstream of
applied statistics. A lot of time, effort and money is put into
geological exploration, but very little into statistical explora-
tion. Of course the data are always looked at and summarized,
21
G. Verly et al. (eds.). Geostatistics for Natural Resources Characterization, Part 1, 21-44.
© 1984 by D. Reidel Publishing Company.
22 N.CRESSIE
Data analytic methods that graph and summarize the way the
observations act and interact, should first of all be relatively
unaffected by the presence of contamination, and secondly be able
to highlight the contaminant from the rest. Such methods are
called resistant, and that part of any geostatistical analysis
which uses them, can safely be called resistant geostatistics.
There is a very important philosophical point to be made here,
namely that it is the underlying Gaussian component which is of
primary interest, and ultimately i t is its parameters which need
to be estimated.
The present paper will use coal ash data from Gomez and Hazen
(1970, Tables 19 and 20) to illustrate the techniques. The data
are discussed, graphed and summarized from an exploratory data
analysis (EDA) point of view (Tukey, 1977) in Section 2. Section
3 introduces more sophisticated data analytic techniques, such as
transformations and a new plot which detects small pockets of non-
stationarity (called a pocket plot). Further exploitation of the
spatial nature of the data allows simple adjustments to be made
for non-stationarity. Usual objections that the bias of the class-
ical variogram estimator becomes prohibitive, are countered by
appropriate first moment calculations. These two aspects are con-
sidered in Section 4. Section 5 compares the residual variogram
to the original variogram, and indicates what conclusions are
possible from the coal ash data.
Data on coal ash obtained from Gomez and Hazen (1970, Tables
19 and 20), for the Robena Mine Property in Greene County, Penn-
sylvania, serve to illustrate the methods of this paper. Other
data sets have received equal treatment with similar success; the
only requirement really is that the observations occur on a regular
or near to regular grid. We have used the 208 coal ash measurements
at core locations with west co-ordinates greater than 64,000 feet;
spatially this defines an approximately square grid, with 2,500
feet spacing, running roughly SW-NE and NW-SE. Figure 1 displays
the core locations. In order to talk easily and sensibly about
directions, we decided to reorient the grid so that it appears to
run in an E-W, N-S direction. Figure 2 shows the coal ash measure-
ments and their spatial location; comparison of Figures 1 and 2
will show how the grid was reoriented.
24 N.CRESSIE
24 ~
:<l
\:)
23 B. 59 900 11.86 8 '11 9 99 en
:<l
22 11.62 10.91 8. 76 8.89 9 10 762 'I. 65 m
en
21 10. 39 10.65 10.36 'I. 58 10.66 8. 92 7 80 7 84 'I 03 B. 60 ..,en
;.-
20 9. 79 9. 06 10.70 11. 21 8.'18 9.27 S.19 7. B8 7.61 8 20 B.77 Z
..,
19 10 74 12. 80 10.03 9.36 8 57 9 01 9. 04 7.28 9. 58 9. 69 9.96 9.91
~
18 11 21 9.89 10.34 8. 20 9. 82 10.06 8.58 8. 8'1 8.04 7.04 8.81 7.95 o
17 9 97 9. 70 9.B4 10.29 9.84 10.01 9 01 7.68 9.25 7.83 q, 14 7. 03 " 07 ~
10 II. 17 10. 14 9. 93 10, 27 10.21 1109 10.63 8. 82 10.18 9.34 B.61 ..,enr.;
15 '1.92 10.82 11.6S 8.'16 '1.88 8.90 10.18 'I. 34 10. 56 9.06 en
14 10.21 10. 73 'I. 46 9.35 9 78 10.38 9. 79 8.91 '9 22 11 43
3
10.'13 10.'14
'1.64
9. 53
9.52
10.01
10.00
10.27
12.05
9. 59
9. 63
9.82 7.81
w+, s
2 '1.29 B 75 8 96 8. 27 8. 14
Figure 2: Reoriented core locations and core values of coal ash data. N
V.
26 N.CRESSIE
61
6
7 003
7 66678888899
8 00111222222234
8 56666666788888899999999
9 000000001111122222223333333444444
9 555555666666666778888888888999999999
10 000000001111111222222333334444444
10 56666667777788888899999
11 0000011122222223344
11 5666689
12
12 568
13 11
13
14
14
15
15
16
16
171
17 6
Figure 3: Stem and leaf plot of coal ash data.
n sgn(Y. -~)
1 1
Y ~+-.
L (2.1)
n
i=l
TOWARDS RESISTANT GEOSTATISTICS 27
ROW SUMMARI ES
8 9 10 II 12
..• .-
x x x J( }Ii
• .. 1\
x x .x '" x ;I(. )C.
•
)(:XXXXX)(X~:x.
·lI
..
x x x x~ ~ x J( ~ ~ J(
~ ~
-.; x x
~ x ~
'lit
x x
lC It
~
.x
X
x
. '" ""
•
".
X X ~ ~ ~ ~ ~ ~ ~ eI'
x x x x x ••
."
~ ~ ~ ~ X
:to: J( )C. 'J{
•
.
x X ~ X x ~ x ~ ~ ~ x If;
•
"
x ~ x ~ ~ x )C. x w X X
)Ii
)
)C
X
J(
)(
~
~
"
~
~
x )( x
)C
~
~
J(
)C.
x x x
)II;
~ ~ •
·If
..
".
X K ~
x )( x x )C
t ~ ~
)(
)C.
)C.
)C.
ell
• iii
..
)C
)C )( )( )( x )( :IC. X
••
)C. )C. )( :It )(
• "
,+,
jI( ~ )t' )l )I(.
>< " x
"•
COLVMN SUMMARI ES
N
• •
.. S
....
•
•• •
... •
• • mean '\ ash
!If
• • median' ash
:.
l l~
'" •
•
n
1
y )1+-
n I (Y.
1
-)1) . (2.2)
q,.
i=l
1 ,gn('il ]
Then y-y n
, where E1 , .. ·,E n have zero mean.
-
i"'l 1 2f()1)
Now assume the E's to be Gaussian. Then it can be shown using the
methods of Section 4,
2
var(Y-Y'j = E-
n 2 - 1) '
(2:. ( 2.3 )
D = In (Y-Y)/(o/'"""O.....,.5=-=7'"="0"""8 ) (2.4)
Row 1 2 3 4 5 6 7 8 9
D -1.54 -0.40 6.12 -0.45 0.35 2.01 -0.56 -0.07 0.63
Row 10 11 12 13 14 15 16 17 18
D -0.18 2.12 0.80 0.46 0.78 0.10 -1.05 -0.60 1.05
Row 19 20 21 22 23
D 0.18 0.35 0.25 1.33 2.47
Column 1 2 3 4 5 6 7 8 9
D 1.11 -0.76 0.78 0.35 2.87 0.02 0.22 1.29 1.23
Column 10 11 12 13 14 15
D 1.03 -0.58 3.17 -1.24 1.39 1.48
(3.1)
(3.2)
C(h), (3.3)
Zt+l 1ir-------------------------------------------------------------,
16
l~
.
.'. .
•
12
• • ••• •
... . ··..'i··; . 'I
,.• ..~.,It",·:.
· •
10 .." ......
-'&"
JIIIf.-, I'.··
.. •
-. '- .. ,.
""
-
~",..
.... ....
• , • i""'"
.
10 12 16 z 18
t
Zt+1 18r-----------------------------------------------------~
16
I~
12
. ..... "
.• I'",,' it ·I'lOr".• •
.. ,t. ~~... .... ..
...'....
10 ~ ~'f-.,,,, ••
...... "•._,..
II .. ..
.' ",...1
v ,-',
•
~ ' ~:' ,
10 12 16
Most of the timewewill assume the C(o) function exists, and in-
deed use it instead of y(o) when convenient. The classical esti-
mator of the variogram proposed by Matheron (1963) is,
Nh
=.l... I (3.4)
Nh i=1
1
2y(h) (med{lz t .+h - Zt. 1~})4/(O.457 + 0.494/N h )· (3.6)
1. 1.
Let us compare these two plots. We all know what the box
plot of a roughly symmetric batch should look like (and besides
the mean should roughly be equal to the median). Therefore any
departures from this are easy to detect by eye. On the contrary
when the batches are inherently skewed, as they are in the vario-
gram cloud, it is difficult to gauge whether the large value is the
result of the skewness or the result of an atypical observation.
It is also clear now when comparing Figures 7 and 8, how the square
root differences do indeed lead to a more res~stant variogram esti-
mator; look at the influence of the outlier 17.61 for h = 6 in the
variogram cloud (Figure 7) and for h = 6 in the square root dif-
ferences cloud (Figure 8).
55
50
45
40
3)
30
25
15
5 6 7 9 11 ----
!3
Figure 7 : Variograrn cloud.
35
3.0
2.5
2.0
8~
1.5
r-
1.0
I-*- -"- ~ f-!!...
f-o.- f-o.- • ~ r.-
Q ........
0.5
10 11 12 13
Figure 8: Square root differences cloud.
34 N,CRESSIE
(3.7)
, •• f- T'
.. m ! 1'" I
!~ r ~ ~§
,... I I •• .. ~
'J !~ 8nJ j: r16 i g~ :I~ ·r r r !~ ·r ~ ~: I I
-0.2
-0.4
-0.• 6
I
2 3 10 11 12 13 14 1
Figure 9: Pocket plot in N-S 16. 17
d~lrectl0n. 18 19 20 21 22 23
..,
'"
36 N.CRESSIE
Y .. Y + (Y. - Y + (Y . - Y
oJ
+ (Y .. - Y. - Y . + Y ) ,
1J 1 0
1J 1 0
°J °°
A plot of the row values and of the column values would yield
a picture almost identical with Figure 4 (where the raw medians
of rows and of columns were considered), showing a definite trend
in the E-W direction, and much less non-stationarity in the N-S
direction. But most importantly, the table of residuals (and
their spatial locations) is available, and can be analysed in its
own right. Furthermore this table is the result of a resistant
analysis. Had the data been analysed by mean polish, we would
have seen the unusualness of the outliers at their spatial loca-
tion diluted, causing it to (confusingly) appear elsewhere in the
table.
-0.00 2.38 -0.74 -1.06 -1.34 -1.22 0.00 -197 0.76 O. :13 o. 23 I. 53 0.25 §
0.90 -0.10 0.00 -1.80 0.34 0.26 -0.03 0.07 0.25 -I. 69 -0. 49 0.00 -0. 18
0.27 -0.46 0.00 0.10 -0.00 0.68 -0.64 -0.78 0.58 -0.41 0.56 -0.17 O. 001-0. 33
~
0.53 -0.67 -0. 14 -0.26 0.00 0.53 0.42 -0.88 0.16 0.:11 -0 43 0.04
-0.09 1.78 0.38 2.05 0.65 -0.90 0.00 -1. 15 -0.07 -1 24 O. 09 -0.76
0.93 -0.23 -0.73 I. 43 -0.20 0.92 0.00 -0.31 1.73 o 00 0.60 -0.39
0.00 -0.50 -0.70 0.00 0.71 -1.39 0.00 0.83 -1.07 92 0.38
0.78 0 •. 95 0.21 0.67 O. 35 0.70 0.36 -0.16 O. 16 -1.03 -0. 82 -1.25 -0.91 -0.34 -1.69 -0. 42 9. 82
COLUMN ALL
Table 2: Median polish of coal ash data.
.......,
38 N.CRESSIE
This "plus one fit" (as it is called by Tukey, 1977 for the
coal ash data, showed k to be zero, so the purely additive decom-
position of (4.1) is reasonable.
1
.-'--.
x x x x X X (4.3)
E 1n~h n-h
I
t=1
(Zi - Z) (Zt+h - Z)
}
*
TOWARDS RESISTANT GEOSTATISTICS
,1
x x x x
(4.5)
x x x x
n-h
E
{
n~h t=1
L
I
C(h) L + C(h) } /n + 0(1/n2 ). (4.6)
k=1
Table 3 compares the O(1/n) term in (4.4) and (4.6) for various
choices of a. Clearly configuration (4.5), from which spatial
relationships are estimated, shows much less of a bias problem.
It is our contention therefore that bias problems, particularly
at higher lags, have beep over-dramatized due to a non-realistic
choice of data configuration (see Matheron, 1971, P.196).
40 N. CRESSIE
h 1 2 345
a =1 (4.4) 1.0 1.0 1.0 1.0 1.0
(4.6) 0.0 0.0 0.0 0.0 0.0
h 1234567
a =5 (4-:-4) 3.8 3.8 3.8 3.8 3.8 3.8 3.8
(4.6) 3.1 1.9 0.8 0.2 0.0 0.0 0.0
h 1 2 3 4 5 6 7 8 9 10
a = 10 (4-:-4) 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5 7.5
(4.6) 7.1 6.1 5.0 3.8 2.7 1.7 0.9 0.4 0.1 0.0
The second way in which our problem is not the one described
by (4.3) and (4.4), is that we are using resistant methods to
estimate the drift. Intuitively this should again amelioriate
the bias problem, since resistant methods do not give rise to
linear cons traints on residuals (e. g. >:: (Zt - Z) = 0).
t 1 n-h
C (h) = n-h L (Zt - med(Zt»(Zt+h - med(Zt»' (4.8)
t=l
(4.9)
f
C(h) - { (4;,) C(O) + 4 2 C (k)
k=l
(4.11)
a =1 a=5 a = 10
IO(l/n) in (4.4) 1 1.0 3.8 7.5
10(17;) in (4.11) 1 0.4 3.0 6.4
ACKNOWLEDGEMENTS
Department of Statistics,
Iowa State University,
Ames, Iowa 50011.
TOWARDS RESISTANT GEOST ATISTICS 43
o a Classical estimator 2Y
)( • Robust estimator if
o
o
o
o o
"
o
o o o
x x "
0
0 0
.
0 0
0 0
0
><
" "
0
" )(
O~--~----~----~--~----~----~----L---~----~--h
REFERENCES
Gomez, M., and Hazen, K. (1970). Evaluating sulfur a.nd ash dis-
tribution in coal seams by statistical response surface
regression analysis. U.S. Bureau of Mines Report, RI ?3??
MARCO ALFARO
Mining Dpt. University 01 Chile. Casilla 2777.
Santiago Chile
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 45-53.
© 1984 by D. Reidel Publishing Company.
46 M.ALFARO
semi ~variogram. a/
characterizes >fri h) - ~(/,) fluctuations and
plays an important role in the control of simulations of the
random function Z(x) • Actually, if UF2 is large, then the diffe~
rences >lJhJ-'({(h) , in average, will also be large.
The behavior pattern of semi-variograms for small values of
distance h is of Ireat importance for Geostatistics. Our p.1TI)ose
is to study the relative fluctuation variance, i.e. CF;/i[hY2 far
small distances: Ih("'o. Consider the following definition:
Definition: Let Z(x) , xe /Rn be an intrinsic random function
defined on a set S C /R". Its semi-variogram is microergodic near
the origin if:
Lim cr/ = 0 (3)
\h(-+o ~2
If the semi-variogram of Z(x) is microergodic, then for a small h
it is possible to determine the theoretical semi-variogram ~(h)
quite accurately by taking the local semi~variogram fram a single
spatial outcome of the random function.
Calculations of the fluguation variance of the local semi-vario~
gram demands that the 4 moments of the increments of the randan
function Z(x) be known. Actually, it follows from (1) and (2)
that: 2 --2 2
OF = E( )lR(/') ) - FhJ
Y<h)= C[ 1 - exp(_h2/o 2 )]
LIe = 5
1.0
5
2
o L
200
The above results are classics (cf. G. Matheron, 1978), and lean
on the Gaussian hipothesis. It is necessary to find more general
relations but the problem is very difficult.
II . Study of the Microergodici ty for the non ;..GaUssian Case.
In the geostatistical practice, simulation of a non-Gaussian
random function are commonly abtained by graphical anamorphosis
i.e transformation by a non-decreasing function.
Let Y(x) be a stationary Gaussian random function such that:
E(Y(x))=o J D Z (Y(,,)=1 1 f(h)= E(Y(x)Y(x+h)
STATISTICAL INFERENCE OF THE SEMIVARIOGRAM AND THE QUADRATIC MODEL 49
Y(~
-2
Z(x)=exp[Y(x))
50.
Lognormal
10.
Gauss
I !f'CX) I = constant, V x
Examine a non-rigurous demonstration of this theorem. Expression
(4) is written as:
C ~-------T----~~~-------
o h
IV. References.
M. Alfaro, 1979 Robustesse des simulations de fonctions
aleatores. Thesis. E.N.S. Mines,
Paris, France. 161 pp.
A. Journel
Ch. Huijbregts 1978 Mining Geostatistics. Academic Press,
N.Y. 600pp.
G. Matheron, 1965 Les variables regionalisees et leur
estimation Masson, Paris, France. 305 pp.
G. Matheron, 1978 Estimer et Choisir. E.N.S. Mines,
Paris, France, 202 pp.
G. Matheron, 1970 Recherche de Schemas polynamiaux E.N.~.
Mines, Paris, France IN.177).
USE OF THE JACKKNIFE METHOD TO ESTIMATE AUTOCORRELATION
FUNCTIONS (OR VARIOGRAMS)
C. F. Chung
Geological Survey of Canada
ABSTRACT
An application of the jackknife method to estimate
autocorrelation functions(or variograms) of stationary random
functions is discussed. Using the jackknife estimators and the
corresponding jackknife variances, the models of the
autocorrelation functions are fitted by the weighted least
squares method. The method is particularly effective to study
robustness of the estimators when the number of data points is
small. The technique is applied to simulated data sets with
known autocorrelation functions.
INTRODUCTION
To study spatial variation or "regionalized"
va ri ab 1es (Mathe ron, 1962) , the theory of stat i ona ry random
functions has been applied. In application, an important but
difficult problem is to postulate the underlying autocorrelation
model and to estimate the unknown parameters in the model from
the data. This paper does not discuss the problems of how to
postulate the model but deals with estimation of its parameters.
Contrary to the importance of the autocorrelation function,
only three studies, Agterberg(1970), Huijbregts(1971), a.nd
Cressie and Hawkins(1980) have been published on the estimation
problems, after their application had been introduced in
geosciences and mining engineering by Matheron(1962). However,
extensive statistical research have been performed on the
estimation of autocorrelation functions in time-series analysis
where the observations are made at regular intervals and the
55
G. Verly et al. (eds.). Geostatistics for Natural Resources Characterization, Part 1, 55-69.
© 1984 by D. Reidel Publishing Company.
56 C.F.CHUNG
( 2)
n
where z =.:;, 'Lz (td
;=1
However, no attempt was made to study the distribution Qf B(s)
because of its complexity. Kendall (1973) has also used B(s) to
compute correlogram in time-series analysis. Although the
distributions of sample correlation coefficients and serial
correlations of time-series data have been extensively studied
(Chapter 32, John~on and Kotz, 1970), the exact characteristics
of the statistic B(d) in (2) are not known.
( 3)
Assume that the nd values (Z(tk )z(t k,) : k = 1,2, ••• ,n and t k,
- tk =d} a re mutua 11.>; independent obse rvat ion s of the random
variable Z(t)Z(t+d), B(d) in (3) is the moment ·estimator of B(d)
and thus, it is a consistent estimator. In addition, if Z(t) is
assumed to represent a normal process, then
( ,)
~7 ~-----------------------------------------------------
~~ IJ\ M
!;rVJVV~
A. _ Ib
;;; ~ l-
___- _
_______
____ _ _ _ _ _ _ _ _ _ _ __
I .. 4 b R 1. ;. I.; 1h Ib ~I) ~4
X-A,rs
:>
<:: co
~ . .•.••• ::z:::..:. ~=--'~ _,
---__
]
0
'.
a: ~'"
a:
0
......
~
f
:.:
I . , ·0 .5 l.0 1.5--"2:0 '''', '>-,,1. 0
UISTAN 'E Js":"',\ .... 1-..
Figure 2. Correlograms (a) and ( b) computed uSlng eq. (4) based on the
realizations shown in Figures l.a and l.b , respectively and shown
as solid lines. The superimposed dotted lines are the models!
(a) B (5) : exp (-25) , (b) B (5) : exp ( _ 52) .
~ ~ j\."._.. . (a )
~o ~ . ~_ ___==~_
---
I]
° a;
.~
(b )
0:
0
U
0
f-
'.
"
0(
-.... - ..
.S 1.0 1. ~ 2. 0 • .5 J.J '" D .S 1.0 \. '> 2.u ).0
DI STANCE DISTANCe
..,~ ...
~'"
C
",0
f-
<
..J"
:>.
>:
~ ~ ~---------------+~--------------------------------
OJ ...
:>
..J
;::'"
..
c
",0
j ",
:> ,
>:
M ... L-____________- -__~----__- -__- - - -__- -__----~----__- -_
en • 0
l( 12 L4 16 19 20 22 2' 26
~-AXIS
Figure 6. Two sequences with random noises. These two realizations are
genera te d from the sequ e nces .i. n Fig. 1 by se lecting 25 out or the
250 values randomly and r ep l ac ing b y rand om numbers fro m
Normal (0 ,9 ) .
o
(a)
~~ I"""~"'"
~......
a; ...
a;
~-
o
g",
o . '. ....
(a)
... .
"'" o -------~-:.--"'>.......,-~r--·--..-'--·:.:.··"'·:,....·~~"'-=·..,....--.....-~~
~
j. . . . . . .
::2 ~ (b)
o
~ga; "';.I""-'~"::-"
'~_'" ~ -. ... - --
~~ ~ :::)
~
- --...
',~
"
.
Ib )
~
<0 ~~.=--
, , < 0 ........ 'I .'=....-.~:-===
.5 LO 1.5 2.0 2.5 3.0 0 .5 l.0 " ~. ~ ._-2~ ---l!:S -. 3.0
DISTI\NCE DISTANCE
o
( b)
o
'~r~'"
a= ••••
o .....
~"' .....
~ . ...................... ..
Co
o .; 1.0
""""'=-"
1.5 2. 0
---=::::::
2. 5 3.0 DlS1A"N~~ .. ,,, .. -
DISTANCE
0: .... ( a)
0:
0
(d)
~~
...; :>
,
'"
-~
0:
~~
0 ..~
............
:> ,
'"
0: ....
........ ..... o
( b)
" '"
.. .
(b)
o o
..... )Son.
)S'" ".
:> ...... :>
D~
Figure 11. Correlograms computed using Cressie - Hawkins ' estimators based on
Fi g. 6.
Figure 12 . Co rrel og r ams computed using jac kkni fed Cressie -H awkins '
estimato rs based on Fiq . 6.
USE OF THE JACKKNIFE METHOD 63
c.)
"'
:>
..J
, ['
:; v
~ ~--------------------------------------------------
'"
:>
• 1
<C
::>
c
"'<Cf- I'l · \, I . \ ",.
(b •
" (
.\.'l '1'.
f.
.'
..J
:>
:(
,
N
til ~
, 0 10 14 16 18 0 22 ~4 26
X-AXIS
Fiqure 13. Two realiza tions of 100 r a ndomly selected points from Fig. 6.
:---::.::.. - - - --
I
l ! J d~I·~~l
-
OrdiniH'y Ja ck kn! te d Ja c kkniled C lt: !. :> ie-
~::o,
ordlnary ere!. 5
I I
es t imdtors M- Ild"'~ III S · I., -
I
eq. 141 est illlato( ::; eq. (g) estillliltors estimat o r :; I lIawkins'
-Lsd
F iy. (14) F i <J. (1 S) riy. (16) fi':l. OJ) fj(J. (18) fig. (19) :
Table 1.
(a)
o
, ,I"" ......
·~r
(a)
Q: ",
C'"
'.~
>I
~ .
'" . ". '.
6
V
.. m
~ <
o
~ ..:j........... rb) .........
, .... . . . .
~
,
'\
(bl
~~ ............. . ......:.\
~o ~ ··r>i· . . . . ..
o .5 1":0 l~ ].0
DISTANCE
" -.,
r'igure 14. Correlograms computed using eg. (4) based on Fig . 13.
F i gure 15. Correlograms computed using jackknife method of eq. (4) based o n
Fig. 13.
USE OF THE JACKKNIFE METHOD 65
computed correlograms.
However, assume that the models of the B (s) are known as:
B (s) = exp( q s
(6)
B (s) = exp ( qs2
for the correlograms in Figures 14, 16 and 18, respectively.
Using the nonlinear least squares method, the q1s in (6) were
estimated (see Table 1). The computations were performed by
using SPSS subprogram NONLINEAR(Robinson, 1979).
On the same data set, the jackknife method was employed to
obtain the correlograms with the jackknife confidence intervals.
These are shown in Figures 15, 17 and 19. The parameters in (6)
were estimated by the weighted nonlinear least squares method
using the jackknife variances as weighting components(see also
Table 1).
From the Table 1, obviously the jackknifed M-estimators most
accurately estimates the original parameters as expected from the
Figures 14-19. The Cressie-Hawkins ' estimators underestimate the
parameters in both of the models. Although the estimators of the
parameters by the ordinary estimators and its jackknife
estimators in model 1 close to the true value q = -2.0, it may be
a just coincident.
(III) Within each of the original 250 realizations shown in
Figure 1, a set of 100 numbers were chosen at random for these
experiments. These series of the 100 values chosen are shown in
Figure 20. Using equation (4), the robust and the
Cressie-Hawkin's methods the correlograms were computed and are
shown in Figures 21,22 and 23, respectively. All three methods
result very similar correlograms.
Since the selections were made at random, the numbers of
pairs of observations used to compute the correlograms, are
similar to each other and these 100 points are identical 100
points as in (II).
CONCLUDING REMARKS
Much work remains to be done to understand the exact or an
approximate distribution of the 6(s) in (4). In addition further
study is required to improve upon the procedure of computing
correlogram by equation (4).
Several assumptions which are unrealistic in practice have
been made in this study. For instance, most spatial data in the
66 C. F.CHUNG
Co)
'"o
0: .~
Ib)
'"""
.5 3. 0
-1 ---' -
o
(a ) , " ' - ~ _"- - "',
~·l~
' o\\J\J~
Co )
.......
,
0:
0:
0
g
- ~
on
... -~~::~ ~-
,
,
( b)
~
:0
,
o+-----~--~----~~~~----~~-- ""
.5 1. 5 2 .0 2.5 1.0 .5 J.O
Fi gu re lB . Correlograms compu ted using Cressie - Hawk ins ' estimators based on
Fi9. 13 .
~ .. (al
:; N
.,..
.J
IJ·
~ o ~, I\' \
:5=> ",.
I r f.f-J~.
. .\
~~
~ "
>:
~ ~L- _____________________________________________________
"'V
::>
.J
(bl
:; '" ,.
.:5 ", - ., " .I \. .\ .
'\
o
",0 ro- '
,. . .\ '\ I .
::> .
>:
~ ~ ~-----------------------------------------------------
10 12 14 18 20 22 24 26
X-AXIS
Figure 20. Two real i zations of 100 randomly selected points from Fig. 1.
Figu r e 21. Correlograms computed using eg. (4) based on Fig . 20.
o
(b) "
n~ /
......
" ........ :....
.S L.O 1.~O 2: ; / J . O o .5 1. 0 1.5 2.0 2.5 J.O
Or STANCE ~ OISTA CE
Fi gure 22. Co rre log rams computed using M-estirnators based on Fig. 20.
(7)
( 8)
ABSTRACT
71
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 7J -89.
© 1984 by D. Reidel Publishing Company.
72 M. DAGBERT ET AL.
1 - INTRODUCTION
A more general method would use both the hanging wall and
the footwall limit of the stratum where point M is located (fi-
gure 10 again). If point M falls between two sections, the li-
near interpolation can also be extended to the interval between
the two sections (figure 11).
Thus the idea before computing any variogram using the usual
coordinate reference system is to estimate the orientation of the
"natural" reference system at any point in the deposit. One
vector of this system is on the stratigraphic axis. To determine
that vector, we can go back to the concept of stratigraphic score
(SS) which has already been used in the limited unfolding tech-
nique. An SS value can be assigned to any point in the deposit
with coordinates x, y, and z. This defines a regionalized varia-
ble SS (x,y,z). The stratigraphic axis at point M is parallel
to the gradient vector of the SS function at point M:
N=
z
78 M. DAGBERT ET AL.
-+ -+ -+
N N I\. N
Y z x
4 - CONCLUSIONS
Acknowledgements
/
--~:..--- DOC = 0
/
/
/
55: I
55: 2
55 = 3
55: 4
5S: 5
55: 6
55: 7
55: 8
DP
SS(M) = 4 + AM/AB
A2
SS(Ml) 4 + A1Ml/A1Bl
SS(M2) 4 + A2M2/A2B2
4
./
./
./
DP
SECTION
4
4.4
4.6
5
I
I
V> I
x DDC (M) I
<l
...---- -~----
M
W
:lo::: I DDC (M')
a::
l- 01
V> 4.4 ~ SS::; 4.6
~I
81
DOWN-DIP AXIS
SLICE MAP
~ I h)
--
sc
.3
()2=.24 . / . / - - AV
- - - - - - - - - -- - -;---- - - -
.2
(505)
" ",(1589)._---=_ _ ----D~ ./
...- --
'-~
./
.....--
.......- . / '
.1
.--- (1488)
(1130)
100 200 h
1m)
Xlh)
.3
SC
() 2 = 24 .- ---
- - - - - - - - -
-- - - -
-
~-
--
AV
2 ./' --.....
--..... ......... DC
.1
100 200 h
1m)
Figure 14 - Slice variograms derived from the same data as
in figure 13 but using a more conventional
reference system parallel to the average dip
plane of the stratum. SC is still the strike
variogram, DC is the variogram computed along
a direction dipping 45° to the south. AV is the
omnidirectional variogram in the surface defined
by the SC and DC axis. Numbers of pairs are
shown between parenthesis.
COMPUTING VARIOGRAMS 87
References
(2) Rendu, J.M., and Readdy, L., "Geology and the semivariogram.
A critical relationship", 1982, 17th APCOMP Symposium, AIME,
New York, pp. 771-783.
P.A. Dowd
G. Verly et al. (eds.), Geostatistics [or Natural Resources Characterization, Part 1, 91-106.
© 1984 by D. Reidel Publishing Company.
92 P.A.DOWD
arise in the same way as they do in real data sets. The major
problem in (i) is in determining a priori whether a real data set
is compatible with the model. Non-stationarity, or drift can
only be assessed from estimates based on discrete data and its
assumed overall absence on the scale required does not imply its
absence at all other scales. In practice, for any given value of
h the mean difference [z(xi) - z(xi+h)] will differ from zero
even though, on average, over all values of h it is zero as
illustrated in figure 1.
Figure 1
Note that under this assumption the bias correction for (3)
for q =0.5, for large numbers of samples, is straightforward and
is obtained by dividing y* by 0.455.
Their scale estimateqfor the differences is adapted from
Huber's scale estimator (Huber, 1964):
The variogram estimate is determined iteratively from:
1
n
2{
n
L
i=l \jJc
y.
\l2Y)<h) Bc }
where y. = z(x.) -z(Xi+h)
1 1
c x > c
Note that taking the median of the Y*h' (7) can be rewritten
as:
1 2 14
.457 {median[(z(xi)-z(xi+h»
J4}
1 n 2
2YD*(h) = {- I I z(x.) - z(x.+h) I} (8)
n i=l 1 1
or
y(h) is the median of the Yi(h) and MAD is the median
absolute deviation f rom the med ian:
} (lOb)
median [y. (h) - y(h)
1
I
Then n. (h) L [yo (h) -y(h)J 2 (1-wh 4
1 i El 1 1
2y i, (h) (lOa and b)
L
20 . ' -".
.04
.03
10 M£AN 8.77
Co 6.0
.02
VARIANCE lQ.08
C 13 .1
a =7500 N/ S
4000
h.103
10 15 20
the local geology but from regional geology long range structures
with good continuity were expected. After a laborious division
of the deposit into various combinations of horizontal "slices"
and calculation of associated variograms the cause of the random
effects was identified as several relatively thin layers of
supposedly weathered material. The analysis is summarised in
figure 4. Variograms in all but the "weathered" zones showed
exceptionally good spherical structures and a "global" spherical
model variogram was ultimately fitted with parameters as shown in
figure 2.
A summary of the performance of the estimators is given in
table 1.
uv:; 1 2 3 4 1 2 3 4
Y*ch (7) 14.5 15.0 15.0 18.4 18.6 13.4 17.3 16.1
y\ (lOa) k=6 14.5 13.9 15.5 15.2 18.6 14.0 15.8 17.4
v\ (lOa) k=9 13.0 12.7 14.0 13.8 16.8 12.9 14.2 14.3
y\ (lOb) k=6 9.8 12.7 11.4 15.1 14.1 12.0 14.4 16.0
Y\ (lOb) k=9 9.3 16.7 10.6 12.8 13.5 11.4 13.0 13.4
Table 1
Figure 4
Note:
(i) the similarity in results for Ych* and YL*
(ii) only YH* and YL* produce results which are significantly
different from y*
(iii) The best estimator is YH~' in 9b, i.e. the median absolute
deviation from the median, although none of the estimators
was particularly good in detecting the East West structure
~iv) the difference between the performance of (9a) and (9b) may
indicate that" localised" drift effec ts have an important
effect on estimates although this is not apparent with (lOa)
and (lOb).
Example 2: A coal deposit
The second set of data came from a southern African coal
deposit in the Karroo Series. The variable studied was % volatile
THE VARIOGRAM AND KRIGING 99
MEAN 45.8
An optimal
variogram model was
VARIANC! 128.5
determined by rigor-
ous cross-validation
40 by kriging each value
in the data set with
30
varying arrornts of da ta
20 removed. Variogram
parameters were
10
adjusted until all
.n
the criteria for an
optimal model were
Figure 5: Histogram of accumulations of satisfied i.e. global
volatl.le matter and conditional un-
biasedness and equal-
ity of experimental and theoretical estimation variances. The
optimal model is significantly different from the experimental
variogram and is shown in figure 6.
Each of the variogram estimators was used to estimate the
variogram and the results are summarised in table 2 and figure 7.
45
40
35
. .
.:...."!-
30 .,x
25 '~"~"'~\ "~" ~> o ... b
" Q
6 7
Figure 6 Figure 7
Taking the model as the criterion for performance the best overall
estimator is (9b), i.e. the median absolute deviation from the
100 P.A.DOWD
LAG INI'ERVAL 1 2 3 4 5 6 7 8 9
Model 10.6 14.0 15.4 16.4 17.2 18.4 19.2 19.6 20.0
Y* (2) 12.5 23.1 29.9 31.0 37.4 39.9 36.5 48.7 44.0
Ych* (7) 9.9 18.7 25.4 26.3 30.8 29.9 33.5 43.9 38.5
YH* (9a) 10.6 15.2 23.2 21.1 26.6 19.9 24.3 29.4 35.9
YH* (9b) 7.6 11. 9 17.2 18.8 22.1 20.2 26.9 27.4 25.3
YL* (lOa) k-6 9.9 19.7 29.5 29.2 34.9 31. 6 33.9 45.7 41. 8
YL* (lOa) k-9 10.3 20.6 29.1 29.3 35.4 33.7 33.3 45.4 41.0
,
yL*(lOb) k-" B.o 15.5 26.1 25.7 31. 4 30.4 32.7 41.1 36,6 ,
y L*(10b) k-9 8.9 17.5 26.5 26.9 32.0 32.2 31.8 41. 8 37.6
I
IV - KRIGING
IV-I Re s is tance
As a linear estimator kriging is not resistant to outliers:
a given data and estimation configuration will give the same
weights {Ai} regardless of the data values occurring at {x.}.
1
kriging equations.
To obtain the weighted median first sort the zi values into
ascending order and then interpret the kriging weights as
relative frequencies of occurrence. The weighted median 1S then
the zi value, or the interpolated value, with cumulative
frequency 0.5.
In Henleys example:
ordered data values 1.0 1.1 3.6 4.0 5.2 6.1 7.7
weights 0.1 0.3 0.05 0.15 0.1 0.25 0.05
cumulative weights 0.10 0.40 0.45 0.60 0.70 0.95 1.00
Linear interpolation gives the weighted median as 3.7.
There is, however, a difficulty in this interpretation when
negative kriging weights occur. A rough solution is to replace
each weight Ai by:
1 A.I
1
8.0".
x
w .~
J....
1.2'\
I ·rL
lL.=:J"""'.~.4.
1.7 ..
X
.
10.2 "
Q 8
I. ."
8."
x
lriled veipted avera,... 11.6% lrieed veiCbte.d averase .. 7.9% (riled weighted .. verq,e .. 8.91
(riled wipted .uian .. 7.41 lriled _i,hted _dian .. 9.7% (riled weighted -.;Iiaa .. 9.1%
2
where s. var{dx.) }
1 1
LA. 1
j J
with kriging variance: (12)
L Ly. + 11 - Yv , v
. J IV
J
2
[ z (X.) - z* (x . ) J
1 1 2 2
if [z(X.) -Z*(X.)J >kOK'
1 1 1
o otherwise
Steps (ii) to (v) involve only marginally more computation
than that entailed in the standard cross-validation kriging
checks of variogram models. Solution time for the kriging
equations in (12) is exactly the same as for the usual kriging
system.
Note that a resistant estimator, such as the median, must
be used in step (ii) otherwise all values in the neighbourhood
of an outlier would themselves be classed as outliers.
As with traditional kriging this technique is an exact
interpolator.
Note also that the three assumptions relating to the error
are all automatically satisfied by kriging.
The value of k and c will depend on the data and the
required amount of weighting of outliers. Values of k=6 to 9
and c = 0.5 to 1 have given good results in several practical
tests.
As a simple example of the effect of sl on the weight, Ai,
attached to sample no i, consider the example shown in figure 9
where samples 1, 2 and 3 have fixed values and the value of
sample 4 is allowed to vary.
,
7.4%
.to,t
I
0.3
!. r 1_~~2,% 0..
0.1
....l, 1 2 3 .. 5 10
VALUE OF SAMPLE No 4
Figure 9 Figure 10
The variogram model assumed is the same as that given Ln
figure 8; the value of k is 6 and c = 0.5. Figure 10 shows the
weight attached to sample no. 4 as a function of its value. Note
that the weight increases with the sample value until the latter
THE VARIOGRAM AND KRIGING 105
v - CONCLUSIONS
ACKNOWLEDGEMENT
REFERENCES
Henning Omre
ABSTRACT
E{Z(x)} = m; all x E A
2
Var{Z(x)} = a; all x £ A
Figure 1.
Isograph representation
of fxh (zl,z2) ·
o
The coefficients are included to ensure scale invariance. The
corresponding pdf, g(v,w), is graphically represented by rotating
the axis in the (zl,z2)-system by n/4 in figure 1.
00
2
y(h o ) E{W 2 } J E{W IV=u}.Prob{V=u} du
-00
00
Dh o {(i,j)IZ(x.),Z(x.)
1 J
£ S; xl' - X
j
= ho }
Nho no. of elements in ~o
In reality deviations from the Ideal Case will occur. Three main
deviation types can be identified. They are in declining order
of importance:
• Distributional Deviations.
The (Z(x),Z(x+h o ») cannot be properly represented by a bi-
normal distribution. The deviations involve characteristics
of the phenomena itself, like Z(x) is non-negative for all x
or the univariate distribution of Z(x) is skewed and heavy-
tailed. From probability theory it is known that univariate
normality does not imply bivariate normality. Hence a uni-
variate transformation of the regionalized variable to
univariate normality does not necessarily remove distribu-
tional deviations. Hence in most cases the conditional
variogram values will be functions of the variable values.
• Sampling Deviations.
The observations {Z(x.);i=l,N} are not sampled in a ran-
domized way. Biased ~ampling has taken place, and the loca-
tions of the samples tend to be clustered in high or low
value areas. These deviations do not include the charac-
teristics of the phenomena, but only the sampling procedure.
In practice, clustered sampling can be seen as an indicator
of biased sampling. This, in turn, will make the traditional
estimator for the mean biased. The influence on the estima-
tes of the variogram values is not so obvious. In biased
sampling the observation pairs used in the estimation are
not representative of the values of the regionalized variable.
If the conditional variogram value is a function of v o ' as
generally is the case, then the traditional estimator based
on biased sampling will probably be biased. The degree of
biasedness in the traditional estimator depends on how biased
the sampling is, and on the functional dependence on vo.
• Outlier Deviations.
The regionalized variable cannot be expected to be precisely
observed always. Usually the observations are correct, but
once in a while an erroneous observation is made. These
deviations only include the sampling procedure. In practice,
112 H.OMRE
N
[F(z)js L 8i ' I (x i ;z) all z
i=1
where
the subscript outside the paranthesis refers to the set on
which the estimator is based.
I(x'z) = { 1 if Z(x) ( z
, 0 else
8.;i=1,N are known weights assigned to each element in S.
l
THE VARIOGRAM AND ITS ESTIMATION 113
Dh : {(i,J·)IZ(x.),Z(x.)
o 1 J
E S; X. - x.
1 J
= ho or x.1 - x.J -h o }
Nh
o
= no. of elements in Dh 0
In the following discussion Z(X j ) and Z(x.) are interchangeable,
hence the symmetric representatIon of paits in ~. The corres-
o
ponding set of observation-pairs is:
Sh : {(Z(x.),Z(x.»)I(i,j) E Dh }
o 1 J 0
.. . .
between figure 1 and figure
2. The bivariate observations ,
in the ho-seattergram is • • ""
."....
....,/..
~/
obtained from sh while the "
marginal cdf's a~e defined by ,,"
. ...
the estimate [F(z)]. It is •
worth noting that tfte margi- .~,
nal cdf can be more precisely
estimated than the bivariate
".'"
characteristics since the
elements in Sh is constructed
1 --------_
from a subset 8f S. [tIz"l,
Figure 2.
An example of a ho-scattergram
which represent Fh o (zl,z2).
114 H.OMRE
L n. ex.IJ"II(x.,x.;zl,z2)
. .) E_ J
; all zl,z2
( I,J no
1
where
I if Z(x) ~ zl and Z(y) ~ z
II(x,y;zl'z2) = { 0 else 2
Lex. . 1
(i,j)EDh IJ
o
a .. ~ 0 ; <ill (i,j) E Dh
IJ 0
Now two estimators for F(z) are available. One based on the
complete set of observations, S, and the other based on Sh •
The weights in the latter are considered unknown. Given tRe
realization of S and hence Sh , the weights a .. ; (i,j) E Dh
may be determined by minimiziRg the deviation 1 tletween the e~ti-
mates [F(z)] and [F(z)] under the constraints given on the
•
weIghts. S sh 0
i-I
L
j=1
Ci(J') i=2 , Nh-1
0
The weights, ;(.); i=l; N, are obtained by letting the cdf step-
function definea by the elements in sh and their associat~d
weights deviate as little as possible ~rom the estimate [F(z)].
There is a one to one correspondance between the elements in S'~
and Sh and also between their respective weights. Hence a solu£
tion f8r a .. ;_(i,j)£Dh is defined through ;('); i=l, Nh • Denote
the soluti6rl a .. ; (i,j~ £ Dt It is obvious 1 that this ~olution
~ld
meets t h e requIre .0.
constraInts.
L ; ..• II(x.,x.;zl,z2)
IJ 1 J
(i,j)£Dh
o
From this estimator and the definition of the variogram value,
it is reasonable to define the estimator of the variogram value
as:
2
t L
~
a .. ' (Z(x. )-Z(x.))
(i,j)£Dt1J 1 J
o
116 H.OMRE
1
2Nh o
I (Z(x.)-Z(x.»)2
(i,j)e:Dh 1 J
o
,),(h)
20
• • • • • •
%
• •
a1
•
6.0
•
5.0
•
••••
4.0
5 10 15 20 30 40 50 60 h
3.0
1.0
15.0 ,),(h)
a1
2.0
• • • • •
•
•
m=1.06
1.0 •
10.0
a2 =2.16
•
N=14400
•••
•
5 10 15 20 30 40 50 60 h
5.0
B. Histogram and variogram
of Deposit II.
1.0
Dh :{(i,j)IZ(x.),Z(x.) £ S ; X.-x. ho }
o l J l J
Nh = no. of elements in Dh
o 0
where
{(i,j)lz(x.),Z(x.) £ S x.-x.
l J l J
Nh no. of elements in Dh
o 0
1.0
.5
ill III
U!
.1 ill
--+- --+- -+-- --+- --+- --+- --+-- --+-- -+-- -+-- -T- ---+-- -...,.
5 10 15 20 30 40 50 60 h
Figure 4.A
The distribution of the standardized estimates from the simulation of the Ideal Case.
The traditional, the Cressie and Hawkins and the robust estimator are applied for each distance.
The standardized variogramvalues of the simulated deposit are presented as long horizontal bars.
The .17, .50 and .83 quantiles of the distributions of the estimates are presented, and the
averagevalues of the estimates are marked as tildas. ?=
o
;;:::
~
tTl
THE VARIOGRAM AND ITS ESTIMATION 121
1.0
.5
fIt
.1
!
-+-- --+- --+- --+- --+- --+-- --+- - - - --+-- --+- --+- --+
I I
-+-
5 10 15 ~ ~ ~ ~ ~ h
Figure 4.B
The distribution of the standardized estimates from.the simulation of distributional deviations.
The traditional, the eressie and Hawkins and the robust estimator are applied for each distance.
The standardized variogramvalues of the simulated deposit are presented as long horizontal bars.
The .17, .50 and .83 quantiles of the distribution of the estimates are presented, and the
averagevalue of the estimates are marked as tildas. ?=
~
i'T1
'"
123
THE VARIOGRAM AND ITS ESTIMATION
')'(h)
B.O
7.0
6.0
5.0
4.0
3.0
2.0
1.0
-+--+--+--+-~-+--+--+--+--+--+--+--+
5 10 15 20 30 40 50 60 h
Figure 4.C
The distribution of the estimates from the simulation of distri-
butional and sampling deviations.
The traditional" the Cressie and Hawkins and the robust estimator
are applied for each distance.
The variogramvalues of the simulated deposit are presented as
long horizontal bars.
The .17, .50 and .83 quantiles of the distribution of the esti-
mates are presented, and the averagevalue of the estimates are
marked as tildas.
124 H.OMRE
V. CONCLUSIONS
VI. ACKNOWLEDGEMENTS
VII. REFERENCES
Paul Switzer
Stanford University
1. INTRODUCTION
We consider first the case where the variogram model has been
specified completely up to a scaling constant. This would be the
case, for example, for the family of linear variograms. For a
general stationary spatial process Z(x) let
° AoZ
2
cov(AoZ) - T 0 L:
where
r.. y(x.-x.).
1J 1 J
~ = ToD ;
2
cov(~) = T 01
mxm
2
E (~ .) = 0, E (~ ~) T , for i=l, .. ,' m .
1 1
2
L:~./m,
1
-~
ToO
11
[2ix.-x·iJ
1
.
1 1
2
y(h) = 'z°E[Z(x) - Z(x+h)J
1 for (8h) 2. 1
2
cov(O) = -T oAr(8)A' = T 2~(8)
where
r .. (S) = y(x.-x.; S) .
1J 1 J
and, in particular,
2
Let r i (SO;8) be the rank order (fixed) of ES{6 i (8 0 )}
among the set of m such values for i=l, ... , m. Similarly,
let R. (8 0 ) be the rank order (stochastic) of the observed value
21
of 6. (8) among the set of m such values for i=l, ... , m.
1
Then both rank vectors r(So' S) = [r l (So,8), ... , rm(SO,e)J and
132 P. SWITZER
If the rank vectors are independent then the upper a/2 point of
the random sampling distribution of K is given approximately by
Further, suppose that the data pairs are disjoint and well sepa-
rated from one another so that, for any reasonable value of 8,
the data pairs are uncorrelated. For testing that 8=8 0 , a
transformation to uncorrelated equal variance linear contrasts is
given by
... ,
!.,:
.6.(8 0) = [Z(x.) - Z(x!)]/[Y(x.-x!; 8 0)J2 for i=l, m.
1 1 1 1 1
2
The rank vector, R(8 0), is given by the ranks of the 6 i (8 0)
computed using the observed Z data. The fixed rank vector,
r(8 0 ), does not depend on the data and is given by the ranks of
y(x.-x!;
1 1
4. TESTING RANDOMNESS
2
cov(D) T 02: (n)
where
and
r 1J
.. y(x.-x.)
1 J
~en) T(n)oD
with
2
cov (~en)) T °I
mXm
n
136 P. SWITZER
6. ILLUSTRATION
Then
2 2
6. (n) ~[Z(x.) - Z(x!)J ![n+[x.-x! [J ,
1 1 1 1 1
and the ranking criterion for the fixed rank vector r(n) is a
monotone function of the interpoint distances within pairs.
...
N
1.0
a: 0 . 8
5
4
l
w ~
~
w ~
J
iC
"
..J
"
U
/I)
0 .2 > 1
_I ---.L
0 0 .2 0.4 0 .6 0 .8 1.0 1.2 1.4 0 2 3 4 5
NUGGET EFFECT ,2'1 DISTANCE h
7. CONCLUSION
REFERENCES
Roland Froidevaux
ABSTRACT
INTRODUCTION
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 141-164.
© 1984 by D. Reidel Publishing Company.
142 R. FROIDEV AUX
Freq.
:~ /
(zl
:
: {EmPirical sample grade distribution F
;
~V \ ~
; .•.... ~"'-
:
:
".....
~------~-T--------~~--~~----~~z
zc, Truncated Distribution
Figure 1
T (z )
v c
144 R. FROIDEVAUX
m (z )
V c
i) T"'(z )
V c
ii) m*(z )
V c
CT~ : : : :. F( L, L) (1 )
N
where:
(2 )
n
\)2
EV
(z )
c
E[[mv(Zc)-m'ftv(zc)] 2}:::.:: 1L [m (z )-m ll ( i) (z )] 2 (3)
n i=l V c V c
TABLE I
statistical and structural
Characteristics of the Simulated Deposits.
..•a
U>
I
•s
.3
IV
\\ i
.3
V
" "3--~-~-,-
, ---~-~2 .I !-3---"'2--.~1-z.-."---~---"2 1.-±3--"'---.~I-..-~.--~--'
Figure 2
Experimental \t(zc) normal cases
;;;:
00 00
••
30
'0 .0
1//'/
~'
/,'/
10
I
10 1 '1/
,• X
•• ~/ ~ '(/
•
/ // /
.'
/ :' ,,
, "
/ / /
/
· \jV
1, ........,' ..
,/ '~,
10 _ ... ..... / "
¥
•• •
\ 'I •
\\
, ,,'
~'
...
____ 51011. 5101
__ ..... 10101 . lOu
Figure 3
Experimental \t (zc) lognormal cases
PRECISION OF ESTIMATION OF RECOVERABLE RESERVES 149
and
O(h,z )
c
PRECISION OF ESTIMATION OF RECOVERABLE RESERVES 151
CASE No.3
Normal distribution
x .0.0
a2 ·1.0
0.5
0.0 .......--r-r-...,.........,...--r-...--Zc
+-...,....~---.-.,-...,....
-4 -3 -2 -1 o 2
CASE No.6
Lognormal distribution
X = .51
3.0 a2 •• 22
2.0
1.0
0.0 Zc
o .1 .2 .3 .4 .5 .6 .7 .8 .9 1.0
Figure 4
Conditional dispersion variances
.!>. .:...~ ~
J!> .23
~o
..
.
-10 I• U
-,.
'" '-0
.2 S3
.,
---- .,. •
- so 10 13
>S •1
. 00 1
.,"
0 0 U)
_>s
,.
SO
..
3J
'"
I.
0
0 ... 'Ou
". 20v
P
50 .,. .... 200
P
Figure 5
Condi tional semi-variograms
a(z ) e
-4[1-T(z
c
)J +1
c (4b)
a 2
2 -
o2 (V/O) o (0/0) - O(V,V)
K (z ) A • A . A ·A [~(zc) Fl (L,L1zc) 1
M (z ) (5)
V c 1 2 3 4
(it Fl (L,L) T (z ) ' V c
V c
where:
+--,----.-.,..--,----.- Ze .1 -1---.----.-.,..--.----._ Ie
K'¥(Zc)
50 CASE NO.5 ..
1''''''
/"
10
,/ '
/
/
.1 Ze .1 Ze
K... (Zc l l(ylZel
ClSE NO.3 50
l
5
f
10
~ . .
.5
.1 +-...,...._-~~ _ _ 2, •I +---.----.-.,..---r----.- Ze
3 2 o o 2 6 8 10
Figure 6
Predicted vs . Actual K (z )
Point support (u x Vu ) c
156 R. FROIDEVAUX
K.,(ZC} K~(Z()
.5
.1
\/
'.;<.
I
Z, .1 z,
K... lZd K¥(ZC)
.5
-J ~ . ..
.1 Zc .1 Z,
K..,IZC' K..,IZcJ
CASE No.3 50
5
.5
.J .~ ,~
10
FIgure 7
Predicted vs. Actual KV(zc)
Su x Su SMU .
PRECISION OF ESTIMATION OF RECOVERABLE RESERVES 157
K'I(Zc:l
\.
.5 ., \
i
.1 L _ _ _ _ _ __
5
.,\j.
~" .,.
. 1 '---_ _ _ _ _ _ __ __
1+-------~~~-- lc
". '
50 CASE No. 6
CUE No.3
5
10
.5 1- • .10
., Ze: .1 +-_-_---,.-
3 2 1 0 1 2 o 2 • 6 8 10
Figure 8
Predicted vs . Actual KV(zc)
lOu xlOu SMU .
158 R. FROIDEVAUX
E (z )
v c
Table II
F (z )
V c
prob{Z (x)< z }
V c
prob[ Z(x).c:'. z~} F(Z/)
c
where
(f (z -m) +m
- C
Uv
T (z ) T(Z I )::::: N( z I ) IN
V C C C
(f2(z )
V c
Fl (L,L/z c )
Fl (L,L)
F (L,LI Zl )
1 c
Fl (L,L)
and
F(L,Llz / ) N
c M (z I)
V C (6 )
F(L,Llz' )
c M (z') (7)
V c
N (z' )
c
Z I Z
C C
and
M (Zi) 1
V c
(8)
CONCLUSION
may seem (especially the first one and the third one) and,
al though the approximations proposed here are sufficient to
obtain an order of magnitude of the conditional estimation
variance, more research will be required to make them fully
operational.
APPENDIX
m (z )
v c
and
<JV
m+-(Z-m)
(j
and
with
(j
Zl m+(z-m)-
\fv
E[ a; (z(x)-m)+m/ Z(x) ~ z~ J
ACKNOWLEDGMENT
REFERENCES
8ruce E. 8uxton
ABSTRACT
INTRODUCTION
To date, not a great deal of work has been done in the area of
assessing the precision of recoverable reserve estimates,
although some important advances are beginning to be made. For
example, see Dagbert and Myers, 1982; Harrison, 1983; and
Fro i devaux, 1984.
G. Verly et at. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 165-183.
© 1984 by D. Reidel Publishing Company.
166 B.E.BUXTON
Since each datum is used for the estimation of only one grid
cell, the elementary estimation errors for the N cells can be
considered to be approximately independent. The estimation
variance of estimating the average value of z(x) over a grid
cell by a randomly located datum is SImply the dispersion variance
of a point (core support) within a square of dimension C,
D~(o/C). And then, if the average value of z(x) over the entire
domain D is estimated by the average of the N elementary estimates
(i.e. the average of the N data), then the global estimation
variance is approximately
(1)
D
in situ point dispersion variance,
D
The point recovery functions can be defined in terms of the
previous random functions. as follows:
D
recovered quantity of metal,
o
recovered ore grade.
mCz) = QCz)/T(z). (\ 0)
168 B. E.BUXTON
~~(C.C.z) ( 11)
N
~ ZJ:' t ' ( z)
}
= s ( z) m( z) I A( z) ( 14)
Vex)
in situ average ore grade (as before),
m = ~~ZV(U)dU = E(Z(x» ( 16 )
D
in situ block dispersion variance,
2 2-
Cf" = Dz.(V/D). ( 17)
recovered tonnage,
(20)
Q,,(z') = TV(z')mV(z,)
= mt- ~~1TCZ) + ;:Q(z).
(22)
The two simulated deposits used for this study are called ST2B
and ST2C. Both deposits are two-dimensional, 220 m in the north-
south (NS) direction, and llurn in the east-west (EW) direction
(see Figure I). Core-support ore grades have been unconditionally
simulated (see Journel and Huijbregts, 1978, pp.491-554) at 1m
spacing for ST2B; so ST2B is exhaustively informed by 24200
core-support grades. ST2C was then created by simply rearranging
fhe 24200 core-support grades in space. Therefore, ST2B and
ST2C have exactly the same univariate core-support grade distribu-
tion, as shown in Figure 1, but different variograms and multi-
variate grade distributions.
220m
1 .."'-
' • . M 1_
1• • 4! I "
'
:::!! ::
H • • 749%
< .. 453%Z.
N ~ 24200
''' , 191-
U.n.l"
II . Ml_ SpiKe of 22.2"
tJ-
11 . _1 I.'
II .•' ,_
(Grades: 0.)
i-llm x 11m BLOCK U.nl_
" ," l_
It .,. I-
I ,'" , _
to ."
•. U ,-
1_
STRATIFYING 1' . .. 1_
.. . . , 1 -
~ ........... .
I ,U 1M
.... ·"1-
) ,n
. ,- ----
-----
l .. _ _ _ ~ _ _ _ a _ _ _
: J .l '
I." ••
J ............... u
_ ••_ _ ......_ . - -
. . . ._ •• _ •• _ _
::~ 1
~:-~:::::::::::::::=...:::.:===-
• •• t - -........_ . - _ ......· .._·-
: I ."
. . ... 1
I
o 2 " ORE
Because ST2B and ST2C are simulated deposits, all of the point
grades, and hence, all of the block grades are known. Therefore,
every point or block grade above cutoff is known for all cutoffs;
and 50, the true recovery functions can be calculated exactly.
Eleven point-grade cutoffs have been selected, and the true point
recovery functions for those cutoffs have been calculated and are
presented in Table 1. (Note that the first cutoff corresponds
to complete recovery.) Of course, true point recovery is the
same for ST2B and ST2C. These eleven point-grade cutoffs corre-
spond to eleven block-grade cutoffs for ST2B and for ST2C through
the affine correction given by expression (18). Notice that the
affine corrections for ST2B and ST2C are different since the in
situ block variances are different. The true block recovery func-
tions for the block-grade cutoffs have been calculated for ST2B
and ST'2C and are presented in Table 1. (Again the first cutoff
corresponds to complete recovery.)
To obtain the true GEV's, the ST2B and ST2C deposits must be
repeatedly sampled. As discussed earlier. the predictive formulas
for GEV assume a random stratified sam~ling plan. Therefore,
for this stUdy, 100 repeated samplings have been taken from 5T2B
and ST2C, where each sampling selects SO random stratified data
ESTIMATION VARIANCE OF RECOVERABLE RESERVE ESTIMATES 173
-1. 1.00 .749 .749 -1. 1.00 .749 .749 -1. 1.00 .749 .749
O. .778 .749 .96] .134 .840 .740 .881 .04] .810 .749 .925
.175 .714 .744 1.04] .278 .750 .722 .962 .208 .725 .7]9 1.020
.325 .655 .729 1.114 .401 .665 .692 1.040 .]49 .670 .72] 1.080
.510 .575 .696 1.211 .55] .590 .655 1.111 .524 .560 .674 1.204
.740 .466 .628 1.]47 .742 .4]0 .552 1.285 .741 .440 .598 1.]59
.905 .]90 .565 1.450 .877 .370 .505 1.]65 .B96 .]75 .544 1.451
1.0] ,]22 .500 1.551 .9BO .]10 .449 1.449 1.014 .]25 .497 1.529
1.14 .Z71 .445 1.6]9 1.070 .255 .]9Z 1.536 1.1IB .305 .476 1.560
I.ZB .215 .377 1.75] 1.185 .230 .]6] 1.578 1.250 .Z]5 .393 1.671
1. 73 .093 .196 Z.I11 1.555 .100 .179 1.792 1.674 .100 .197 1.969
--
Table 1. True Recovery Functions.
174 B.E.BUXTON
Figure 2 shows that for these two deposits the affine correction
provides a very satisfactory approximation of the block recovery
functions. The affine correction performs well in part due to
the fact that the change in dispersion variance from point
support to block support is not too large - 33X for 5T2B and
llX for 5T2C. However the affine correction also performs well
in spite of the fact that the multivariate distributions of
5T28 and 5T2C are not simple.· 5T28 has multigaussian tendencies
ESTIMATION VARIANCE OF RECOVERABLE RESERVE ESTIMATES 175
1.0
ro
~ 0.8
(j)
w
~ 0.6
z
Recovered z
lSI
Tonnage f- 0.4
>
u
W
IX
~ 0.2
--'
CD
0.0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GRADE
1.0
~ 0.8
f-
(j)
:E 0.6
Recovered lJ..
lSI
Quantity 0
Of Metal > 0.4
u
w
IX
~ 0.2
--'
ro
0.0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GRADE
2.5
ro
N
til 2.0
w
0
~ 1.5
l')
Recovered w
Ore Grade ~ 1.0
>
u
W
IX 0.5
~
--'
ro
0.0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUTI2lFF GRADE
1.0
u
~ 0.8
Ul
w
~ 0.6
Recovered z
z
lSI
Tonnage f- 0.4
>
u
w
~
':L 0.2
---'
rn
0.0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUTeJFF GRADE
1.0
~ 0.8
f-
Ul
Recovered I: 0.5
LL
Quantity lSI
0
Of Metal > 0.4
u
w
~
':L 0.2
---'
rn
0.0
0 0.75 1 1.25 1.5 1.75 2
CUTeJFF GRADE
2.5
u
tn 2.0
(\J
w
0
& 1.5
L'J
Recovered w
Ore Grade ~ 1.0
>
u
w
~ 0.5
':L
---'
rn
0.0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GRADE
Global estimation variance has been predicted for ST28 and ST2C
with expressions (11),(12),(13),(23),(24),(25). In all
cases the true recoveries, in situ mean ore grade, point disper-
sion variance, block dispersion variance, and correlation coeffi-
cients have been used. Also, the exhaustive semi-variograms
(from all 24200 point grades) have been used to calculate the
average semi-variogram values needed in expressions (11),(12). The
exhaustive semi-variograms have been modeled by an isotropic linear
structure. This type of model provides a satisfactory approxima-
tion, as shown in Figure 3, for separation distances up to 31 m
(the maximum separation in the stratifying grid cell). (Only 2
of the 44 sets of semi-variograms are given in Figure 3
as an example.)
The results of this study for the ST2B and ST2C deposits are
presented in Figures 4 to 7. The predicted GEV's are plotted
against their corresponding true values. The GEV's are all given
in terms of relative estimation standard deviation as defined
in expression (26).
NS
~ 0 .8
~
U
~ 0 .6
I-
(f)
(f) ~ ___ -
: --
___ • -~L.INEAR
HODEL
~ 0.4
t!J ~--------------------EW
0::
<
> 0 .2
......
L
W
(f)
~ O . O ~~~~~~~~~~~~~~~~~~~~
N 0 10 20 30 40 50 60 70 80
SEPARAT10N DISTANCE (METERS)
1.0
0
"<t NS
i':
~ 0 .8 .SE
::J
U
U
C\J 0.6 -NE
.....
(f)
(f)
LINE .... R
_ --- HODEL
~ 0 .4
t!J
0:: EW
<
>, 0.2
......
L
W
(f)
..... 0.0
N 0 10 20 30 40 50 60 70 80
SEPARATI0N DIS TANCE (ME ERS)
40
N
30
ID
,-
(\J
(/)
Rgcovergd 0
(/)
20
Tonnage w
Ct:
z
10
"""
"""
0-
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUHlFF GRADE
40
~
CD
30
C"\J
I-
Recovered (/)
1
Quantity 0
(/) 20
Of Metal w
Ct:
:c
I.J..
lSI 10
0
I-
0-
0
0 0.75 I 1.25 1.5 1.75 2
CUT0FF GRADE
10
:::::
8
ID
C"\J
l-
(/)
6
0
Recovered (/)
w
Org Gradg Ct:
4
0
Ct:
l!l
IJ.J
Ct: 2
lSI
I-
0-
DD 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GRADE
40
'" 30
co
(\J
l-
V>
Recovered 0
V>
20
Tonnage w
cr
Z
I- 10
y
-'
co
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GRADE
40
co 30
(\J
I-
V>
Recove r ed 0
Quant i ty V>
w
20
0::
Of Metal
I:
LL
IS) 10
0
y
-'
CD
00 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUT0FF GR"DE
10
CD
8
(\J
.....
V>
0 6
V>
Recollerl!d w
0::
Ore Grade 0 4
0::
l!l
W
0::
IS)
2
y
-'
CD 0
0 0.25 0.5 0.75 I 1.25 1.5 1.75 2
CUT0FF GRADE
40
R 30
u
N
>--
V)
RQcove.-ed D 20
11l
Tonnage w
~
Z
tSl
l- 10
I-
n..
0
0 0.75 1 1. 25 1.5 1.75 ;:>
CUT0FF GRADE
40
'"
u
30
N
I-
RQcovered V)
D
Quant ity V) 20
w
Of Metal ~
1:
U.
tSl 10
0
l-
n..
00 0.25 0.5 0. 75 1 1.25 1.5 1.75 2
CUT0F F GRADE
10
'" 8
u
N
l-
V)
6
D
Recovered V)
w
O.-e G.-ade ~
4
D
~
<!)
w
~ 2
tSl
l-
n..
0
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2
CUl0FF GRADE
40
;::;
30
u
N
>-
(/)
Recovered 0 20
(/)
To nnage w
cr
Z
lSl
>- ]0
:,.::
...J
CD
a
0 1.5 1.75 2
40
'" 30
u
(\J
>-
Recovered (/)
Quantity a
w 20
(/)
Of Metal (Y-
1:
lL.
lSl 10
"
:,.::
...J
CD
0
0 0.75 ] 1.25 1.5 1. 75 2
CUT0FF GRfIOE
]0
'"
u
8
N
>-
(J')
0 6
Recove red (J')
w
Ore Grade (Y-
o 4
(Y-
l!)
W
t:r
lSl
2
:>l
...J
CD
0
0 0 .75 I 1.25 1.5 1.75 2
CUHl FF GR" OE
CONCLUSION
REFERENCES
ABSTRACT
G. Verly et al. (eds.). Geostatistics for Natural Resources Characterization. Part 1, 185-199.
© 1984 by D. Reidel Publishing Company.
186 Y. C. KIM AND E. Y. BAAFI
l
lOX Kriged Value
* * 429* 450* 465* 467* 456* 451*
* * 14* 17* 25* 38* 37* 38*
* * 0* 0* 0* 0* 0* 0*
+----+----+----+----+----+----+----+
* 447* 455* 461* 464* 458* 4SS* 455*
* 65* 40* 17* 13* 47* 60* 64*
+----:---~:---~:---~!_--~!__-~:---2!_--~:
465* 472* 478* 482* 482* 472* 469* 464*
I * U06* 79* 55* 37* 39* 6.3* 77* 78*
* 0* 0* 0* 0* 0* 0* 0* 0*
lOX Kriging
+----
Variance * 542 492 501 ----+----+----+----+----~
500* 505* 502* 486* ~84*
* 101 L78 m52 39* 51* 72* 85* 85*
* a ~ 0 VQ 0 0* 0* 0* 0* 0*
+----+----+---- ----+----+----+----+----+
* 645* b06* 568* 535 513 522* 524* 521* 50S* 507*
* 103* 102* 69* 65 e32 13* 44* 70* 85* 83*
* 0* 0* 0* a a 0* 0* 0* 0* 0*
:-691~24 579:-54"5* 526* 537:-;;2:-53;:-531:-529:
* 80 ~84 71* 49* 28* 19* 42* 61* 71. 72*
• 0 ~ 0 0* 0* 0* 0* 0* O. 0* 0*
:-&87* 628*-;78@[40-523537553-557:-543:-532:
* 63* 69* 51 822 26 rr37 ~41 40. 46* 51*
* 0* 0* 0 0 0 ~ 0 ~ 0 0* O. 0*
+----+----+----
• 666* 606*
--- ----+----+
* 522* 521* 528 544 562 540* 510.
* 65* 67* * 13* 25. 41 'TD35 9 19* 20*
* 0* 0* * 0* 0* 0 ~ 0 1 0 0* 0*
+----+----+----+---+----+----
* 614* 566*
----+----+
* 508* 506* 514* 526* 534* 515* 486*
*. 77* 73* * 39* 33* 38* 38* 25* 23* 10*
* 0* 0* * 0* 0* 0* 0* 0* 0* 0*
+----+----+----+----+----+---+----+---+----+----+.
* 562* 530* 507* 498* 495* 503* 513* 512* 491* 473*
* 81* 13* 62* 43* 11* 23* 41* 41* 36* 28*
• 0* 0* 0* 0* 0* 0* 0* 0* 0* 0*
+----+----+----+---+---.--+--...;.-+---+--+---- ...
(1)
where
(2)
m
Z;
1
=2.:
k-1
AiJc Z1l<. 1-1,2, ... ,n (3)
where Zil, Zi2' ..• , Zim are sample values used to krige block
V. and Aik are the respective weights. In equation (3) above,
i£ is assumed that the same number of points are used to krige
each mining block. If it is difficult to satisfy the above
assumption of using the same m samples to krige each block, it
can still be satisfied by letting the weights of excluded assays
to be zero, at least in theory.
188 Y. C. KIM AND E. Y. BAAFI
(4)
and
1 2 * 1
-;1 V1 Var [ Zv - Zv
1 1
LOCAL KRIGING VARIANCES FOR SHORT-TERM MINE PLANNING 189
(6)
(7)
f
m
2:
Aik
Cov[ (z" - z"* ). (ZV - Zv* )1 ~ r(x - Zik) dx
i i j j Vi
k-l Vi
f
m
~
+L k-l
Vj
Vj
y(y - Zjk) dy
1
ViV j
ff y(x - y) dxdy
Vi Vj
m m
k-l r- 1
L
n
I1(D) - ~ Vi )leVi) (10)
i-1
Equations (9) and (10) imply that if the same sample points
are used to krige a group of blocks one at a time, it is possible
to obtain the solution of combined kriging system of equations
from the solutions of individual (local) systems of equations.
Once the combined kriging weights and the Lagrange parameters
are known, the combined kriging variance can be computed using
the standard kriging variance relationship given in equation (11).
oK
2 (D) - -C (D,D) - " L..J (11)
..
where C is an average covariance term.
192 Y.C. KIM AND E. Y. BAAFI
The program does not group working areas into disjoint sets
which is needed when the working areas are scattered in the mine.
Such a grouping must be done manually by the user. However, the
program has an option to krige blocks with sample points inside
a polygon whose coordinates are input by the user. This option
ensures that the same sample points are used in the kriging pro-
cess. The program additionally checks that identical sample
points are used in the kriging of individual mining blocks. A
run is terminated when different sample points are found to be
used in the kriging of any block.
+. Sample Locatio~
+
+ + +
+
+
+
'+
+ +
+ + + +
r+;\ ~Polygon A
+
+ +~
t-:;::J+ +
+
~ +
0.+ +
~f\n : \.\
\!-J ~+
+ +~:+10
+ +& +
Polygon B ~
+
5
0 + ++0 + +
-I"
+ +
+
:<
o
~
Fipure 3. Nine ~vorking Areas Enclosed by Two Polygons ~
~
1"'1
"
:<
CD
[:
:!l
LOCAL KRIGING VARIANCES FOR SHORT-TERM MINE PLANNING 195
Yes
CONCLUSION
The reliability of an estimate of a combination of kriged
values is conveniently obtainable provided that the correct
assumptions are made. If the same data points are used to krige
a set of mining blocks, t~e combined kriging variance can be
computed using the individual kriged rE:sults. Hhen mine workings
are scattered in a mine such that it is not possible to krige all
the mine working with the same data points, some grouping of work
ing areas may be needed prior to the use of the local kriging
results.
LOCAL KRIGING VARIANCES FOR SHORT -TERM MINE PLANNING 197
* *
- E[(;'. - z".)(ZV. - z".)]- [E(ZV - ; , ) ] [E(ZV - Zv )] * * (12)
1 lJ J i 1 j j
But,
Hence.
Cov [(ZV *
- Z ) , (ZV - ; , ) ] *
1 V1 j j
- E [(ZV *
- ; ' ) ( L - Z )] * (13)
i i "j Vj
Cov[(;,
i
-z;) ,(;, -z;)]
i j j
-tt k -1 r-1
(14)
But.
(Z
Vi
- zik) (Zv. - Zjr)
J
• v1 v
ijvV
II [Z(x) - Zik1 [Z(y) - zJ·r1 dxdy
i j
JJ
Hence,
.. v \ . ff
i J V V
E[Z(x) Z(y) - Z(x) Zjr -Z(y) Zik
i j
.-1....
ViV j fi
V V
K(x - y) dxdy - !Vi f
Vi
!C(x - Zjr) dx
i j
r
Vj
J
V
K(y - Zik) dy + K(Zik - Zjr) (15)
J
198 Y. C. KIM AND E. Y. BAAFI
(16)
K(n) = K(O) - y(h)
Cav[Cl.; -
i
Z;) ,(ZV
i j
- z; )]-
j
If y(x - y) dxdy
m m
-L L AikAjr Y(Zik - Zjr) (17)
kal r- 1
ACKNOWLEDGEMENT
REFERENCES
Raymond L. Sabourin
1.0 INTRODUCTION
G. Verly et af. (eds.). Geostatistics for Natural Resources Characterization. Part 1. 201-215.
© 1984 by D. Reidel Publishing Company.
202 R. L. SABOURIN
UG( NO
( " ..... ~ . S. \JI .. IEII,I'-ot!
G.. !t<t,UI c:..! .,,,,,11
i. '~.-. i..
border line. This coal field covers an area of almost one thou-
sand square kilometers. However, only a small fraction the total
geological resources can be economically exploited by the actual
mining methods. The classification methodology will be applied
on geological resources of the Short Creek zone without taking
into account economical or mining conditions which would define
the reserves. However, the method is not limited to classify
geological resources and can equally well be applied to the
reserves.
~ r-------------------------------~
~;;' /'
,;
.
Reasonabl y Assured
+ +
,0
.......-----L- - - -- = - - -- - - - - - - - - - - - - - 4".oo
+ +
=0 for
h
0
The method will be applied as proposed and ratios 0.2, 0.4 and
0.8 between a~ and a~ will be considered to arbitrarily define
three KV limits: a~M = 0.2a6, a~I = 0.4a6 and a~F = 0.8 a6for
the measured, indicated and inferred tonnage resources respective-
ly. The values shown below on the contour maps represent twice
the standard deviation of the block KV. Consequently, the above
criteria should be transformed to: CokM = 2/.Cob, CokI = 21.lbb
and 2akF = 2/.8a b• A new set of criteria will be calculated for a
different size of block. The size of the chosen blocks are re-
spectively: 1000 m x 1000 m, 500 m x 500 m a~d 250 m x 250 m. 2he
values of the distribution variances are: ab(lOOO) = 0.91957 m ,
a~(500) = 0.97467 m and a~(250) = 1.00237 m. Table 2 gives the
criteria used to classify the resources of the Short Creek coal
zone.
Blocks Blocks Blocks
(1000 m x 1000 m) (500 m x 500 m) (250 m x 250 m)
Measured (M)
( CokM) 0.8577 0.8830 0.8955
(Tonnage) 105 x 10 6T 89 x 10 6T 81 x 10 6T
M + Indicated (I)
(Co kI) 1.2130 1.2488 1.2664
(Tonnage) 163 x 106T ,154 x 106T 152 x 106T
M + I + Inferred (F)
( CokF) 1. 7154 1. 7661 1.7910
(Tonnage) 213 x 106T 212 x 106T 211 x 106T
+ +
+
Figure 5. Isocurve map of twice the standard deviation (m x 100)
of kriging variances for 1000 m x 1000 m blocks.
These criteria vary slowly with the size of the blocks. This
can be explained by the large value of the range of the variogram
compared to the dimensions of the blocks. Indeed a~ should vary
more rapidly when the range of the variogram is small.
7.0 DISCUSSION
+ +
+ +
Figure 8. Isocurve map of twice the relative preclslon (xlOO) of
kriging estimates for 1000 m x 1000 m blocks.
~~-------,~~~-~~~==~--------,+
.,
""
0-
~ a -
'"
0
I> -
-D
0
M
~
N
X 0
n
~,
'"
0
0
~
0
f::::
- 0
~
~ 0
'E n
~
~ 0
f:::: 0
'" / 1"
'"0 I
'"
0
...,~ / '
0
Q
0.7500 .875 I.oao 1 . • 25 •. 250 •. 375 !.50C 1.625 •. 750 ;.0 ' 5 2.00:
KRiGI G STA O~RD OC fATICN LI ITS (Ml
(4)
The corrected KV limits for the three sizes of block and f',)['
M + Indicated (I)
ah/a6 0.3460 0.3830 0.4
( ;bkI) 1.1281 1.2219 1.2664
Tonnage 146 x 106T 150 x 106T 152 x 106T
M + I + Inferred (F)
a~F/a6 0.7820 0.7943 0.8
(;b kF) 1.6960 1.7598 1. 7910
Tonna~e 211 x 106T 211 x 106T 211 x 106T
The initial ratios of 0.2, 0.4 and 0.8 for a~/ag were arbi-
trarily chosen by the author to show some properties of the
method. The tonnage (mainly in the measured and indicated
resources categories) would be considerably reduced if the
initial ratios would have been applied on points instead of
blocks. In fact, no corrective factors could have been cal-
culated if the nugget effect would be larger. To eliminate
this problem, the initial ratio should never be applied to
points because it involves using the nugget effect. Instead,
the initial value ofat should be equal to the sill value less
the nugget effect, i.e., a~ = C. The above initial ratios
could be legitimately used to define the three categories of
resources.
APPLICATIONS OF A GEOSTATISTICAL METHOD 215
References
I - INTRODUCTION
G. Verly et ai. (eds.) , Geostatistics for Natural Resources Characterization, Part 1,217-230,
© 1984 by D. Reidel Publishing Company.
218 M. DAVID ET AL.
paper will try to do is to cope with the bias and get rid of it
as much as possible. Using simulated deposits, a number of
situations will be reviewed, covering normal, lognormal and
highly skewed distributions. Simple, ordinary and lognormal
kriging will be considered using variograms, relative variograms
or lognormal variograms. The impact of the bias on reserve
estimation will be reviewed and a solution proposed.
The main remark used in this paper will be that ordinary
kriging is equivalent to simple kriging performed with the opti-
mally estimated mean (Matheron, 1971). As for computer reasons
ordinary kriging is performed with a small number of points, it
means that it is equivalent to having estimated the local mean
with the same small number of points and then having done the
simple kriging. We will see what improvement can be gained
simply by getting a better estimate of the mean with a larger
number of samples. It will be seen that this surprisingly
simple tool does get rid of most of the bias. Other possible
avenues to correct the bias will also be discussed at the end.
(4) 2
KS
D2 I Yi cov(Z , Z. )
U Z v l
V i
CONDITIONAL BIAS IN KRIGING 219
(6) i y.J = 1
j
Again Zv and Z*KO have a binormal joint distribution. It
follows, considering ~5) and (6) that
D~tKO _ w
(7) E [ZVIZ~KO] = m + (D2 ) (Z~KO - m)
Z~KO
(9) E[~ -Z*KSK)2] = 0K2S + a20~ + 2a[ cov(L: A.Z.,iii) - Cov(Z ,iii)]
v v m. 11 V
1
The last term vanishes if the same data set is used to
estimate the mean and the block value (case of KO).
If iii is a good estimate (N large) the the last two terms
are very small compared to 0Rs and we may write the relation:
o~S<0RSK<ORo; using properties of the normal distribution (app.3)
cov(Z ,z* )
(10) E[ Z IZ*K ] = m + 2V vKSK (Z* - m)
v v SK DZ* vKSK
2 vKSK 2 2
where cov(Z ,z* ) = DZ'~ - a 0_ - a[2cov(D.Z.,iii)-Cov(Z ,iii)]
v vKSK vKSK mil 1 V
220 M. DAVID ET AL.
The last two terms being small relatively to D~* , the bias
is small and decreases as the precision of the ¥~efmate of the
mean increases with special cases being ordinary kriging and
simple kriging. Comparing equations (7) and (10) leads to
(11) ).l = ao~
This r~ation was derived differently by Matheron (1971).
It expresses the Lagrange multiplier ).l as a fonction of the
imprecision on the estimation of the mean with again, as limiting
cases, ordinary kriging and simple kriging.
Table 1
KO KS KSK
Error 0.0077 0.0040 0.0068
Absolute error 1.705 1.665 1. 678
Squared error 4.659 4.437 4.502
bO global* 3.279 -0.2712 0.1285
b 1 global 0.9346 1.0055 0.9976
bO average* 2.983 -0.477 -0.16
b~ average 0.941 1.010 1.003
r (correlation) 0.7905 0.7999 0.7965
Table 2
KO KS KSK I
£ -0.002 -0.003 -0.002
£2 0.0106 0.0096 0.0098
Td 0.0587 0.0562 0.0571
bn 0.018 (,.005 0.006
b, 0.857 0.965 0.964
Table 3: Global comparison on 100 data sets, of KO, KS and KSK
222 M. DA VlD ET AL.
'( 08"
S
Co: 050
/ C '60
A I 000
.1,) Co 12 1 • h
() Co "0 I I, ) Co.IO
C ' 1. e • 3S C : .2
A 114 10 A = 1;)100 A - 1170
H. h '06 h , .. h
Table 5
KO KSK
t-<:>. T t P % T t P % potentiaLI
10 28 s 'mulat ons: Norm dist ibuti n (50 25)
(140 )0 blo ks of 50 m
70553 -0.1 I
45 13354 50.28 70482 0.3 13472 50.24 71763
47 11843 50.78 44624 1.5 12018 50.73 44779 0.0 47367
48 10497 51.17 33286 2.5 10635 51.14 33405 0.0 36669
I
I
I
49 8818 51.66 23474 4.0 8921 51.65 23596 0.0 27383
50 6956 52.23 15491 6.7 6982 52.24 15626 O. :~ 19621
2 0 20 s 'mulat 'ons; ~ogno nal di ~tribu ion m = O. 2,
(] 1m = 2.0 10000 bloc )
.02 9920 .117 962 2.6 9984 .116 962 2.6 965
.03 9524 .120 861 3.4 9909 .117 862 2.9 875
.05 7995 .134 673 5.8 8855 .126 676 3.2 720
.07 6149 .155 521 9.8 6802 .147 527 2.5 599
.10 4059 .188 356 18.3 4149 .187 362 4.5] 461
Table 6: Recovery for the normal and lognormal case
KO KSK
t", T t P % T t P % Potential
1 0 Abso ~ute v ~riogr ~m
.10 363 .89 288* 3.0 399 .83 290 5.8 294
,20 310 .99 246* 7.0 386 .85 250 6.8 266
.40 231 1.13 169* 24.3 337 .94 183 6.4 225
.60 187 1. 33 136 23.8 248 1.13 131 3.8 193
2 0 Rela ive v ~riogr am
.10 325 .97 283* 2.1 376 .87 289* 3.6 294
.20 275 1.09 244* 6.0 336 .94 249* 5.9 266
.40 192 1. 32 176 20.8 231 1.20 186 11.2 225
.60 159 1.53 147 21.1 156 1.55 149 14.2 193
30 r- Loga ithmi vari ograrn
.10 337 .95 286* 20.7 344 .93 287* 32.8 294
.20 251 1. 20 250 26.5 267 1.14 251 39.8 266
.40 190 1.45 200 36.7 191 1.44 199 53.6 225
.60 143 1.71 159 52.5 157 1.60 157 72.7 193
* The "no selection" option gives better profit than the above.
Table 7: Recovery for the skewed distribution
226 M. DAVID ET AL.
VI - CONCLUSION
ACKNOWLEDGEMENTS
APPENDIX
1 - Preambule
N
Z 'Y. 1
. 1
~
and we have
2 (J2
0
K
=
m ).10
2 - Derivation of formula 9
2 D2 2 Cov(Z v,Z~KSK) + D2
(JKSK Z
v Z~KSK
D2 - 2 Z A. Cov (Z , Z . ) + Z Z A. A. Cov(Z.,Z.)
Z
i j 1 J J
1 V ~
v i 1
3 - Derivation of formula 10
4 - Lognormal KSK
i.e. m = ~ .
1
y. y.
1 1
with ~ Yi =1
i
o~S - a20~
exp[Y* + 2 -aCov(LA.Y.,ill)]
vKSK ill
is unbiased.
This result should be compared with the expressions for
simple lognormal kriging and ordinary lognormal kriging given in
Journel (1978). Again KSK can be viewed as the general expres-
sion for which KS ~nd KO are special cases.
Now, the variance of the kriging error is
with
D~ M2 [exp(D~) 1]
v v
- 2exp(aCov(Y ,m))]}
v
Again, if we compare this expression with the expressions
(35) and (45) presented in Dowd (1982), we see that KSK leads to
a general expression for which KO and KS are particular cases.
This expression also shows that the better the estimate of
the mean (o~ small) the smaller the kriging error.
m
NOTE: (if we don't use Rendu's approximation, there is no
simple expression either for 0 2 * or 0 2 or 0 2 and
zvKO z~KS z~KSK
therefore comparison is impossible to make) .
230 M. DA VlD ET AL.
REFERENCES
ABSTRACT
231
G. Verly et al. reds.;. Geostatistic-I> for Natural Resources Characterization. Part 1, 231 ·244.
© 1984 by D. Reidel Publishing Company.
232 R. 1. BARNES AND T. B. JOHNSON
LX + 1"- = C5 {Z.Zal
l' x = 1 {Z.Zb)
These condi t ions may be written more conveniently in a
partitioned matrix form [4J as
{z.3l
r
sample covariance matrix, L, is positive definitive, the
inverse of ~ necessarily exists and thus,
[f· *
{Z.4)
is guaranteed to exist.
The Positive Kriging problem can be stated as finding the
set of non-negative weights which minimizes the estimation
variance (a quadratic form), subject to the unbiased condition
(a linear constraint). Thus, Positive Kriging is again nothing
more than a quadratic programming problem (it is, however, more
comp 1i cated than the OK problem). The Pos it i ve Kri gi ng problem
can be written in a matrix equation fonn as
234 R. J. BARNES AND T. B. JOHNSON
1:'~ =0 {2.6c}
x ) 0 and 1:) 0 {2.6dt
where
J [~ ]
1
-x
G, 0
-1
01 A
= {2. 7a}
1:
with side constraints
£I~ =0 {2.7b}
xi = 0 or ~i = 0 for all i.
This is known as a complementary slackness condition, and
variables
xi and ~i
are said to be a complementary pair.
where the subscri pts "b" and "n" i ndi cate basic and non-basi c.
Carrying out this rearrangement of rows (and symmetric columns).
- ,l· . -
the Positive Kriging system may be written as
[:]
[iJ.' iJ.' l
- b - n 0 O.le}
x
[:]
and {3. Id}
- >0 1:: > 0
Where
~b is the square matrix of covariances between the samples
associated with the basic x's.
~n is the squa re mat ri x of cova ri ances between the samples
associated with the non-basic x's.
~n is the rectangular matrix of covariances between the
samples associated with basic x's and samples associated
with non-basic x's.
~b is the rectangular matrix of covariances between the
samples associated the the non-basic x's and samples
associ ated with the basic x' s; !t,n and ~b are transposes
of one another.
POSITIVE KRIGING 237
(J
-{l
is the vector of covariances between samples associated
with the non-basic x's and the volume to be estimated.
Multi plyi ng out the expanded system {3.1} and remembering that,
by definition, the non-basic x's and the basi c fl'S are equal to
zero, yields:
~b ~ + ~
" =~ D.Za}
~b ~ + In " - -1n~ = (J
-fl
D.Zb}
I' ~
-b 1 {3. Zc}
l-~b b] [~'I' l~ -]
Rewriting conditions {3.Za} and {3.Zc} in a part it i oned
matrix form yields
D.3}
l'b 0 ," J 1
This system of equations has N+l equations and N+l unknowns, and
may thus be solved directly. In fact, system {3.3} is identical
to the Ordinary Kriging system {Z.3} except it is solely in
terms of the basic x's and the unbiased condition's Lagrange
multiplier. System {3.3} will be called a Reduced Ordinary
Kriging system in the basic variables ~.
4. SIMPLEX-BASED ALGORITHMS
6. COMPUTATIONAL DETAILS
The following proposition presents a useful identity for
the inverse of a symmetric matrix partitioned on its rightmost
column and bottommost row.
POSITIVE KRIGING 241
I·
then the following i dent ity hol ds:
[!
A d
[ =f b' -s
then the following i dent ity hold s :
A-I C + (1/B).9i'
j
PROPOSITION 4 - Given that
- ~ 11 ••• bI .•••
. J
~1mr [<W"
·· d.1'" ~1m . J
b . I ••• b ..••• b. d· l ••• d ..••• c .
. Jm
.J .JJ
·.J ..J J .Jm
amI··· bmj ••• amm em 1 • •• dmj ••• cmm
then the inverse of the matrix after row j and column j have
been permuted to the far right and bottom of the matrix is
242 R. J. BARNES AND T. B. JOHNSON
~ 11 .• , a.lm
bmj
-' -1
.1 J
r :11" ,
c ml '"
:1m
c mm
dl
• J· 1
bj 1'" b. b .. dj 1'" d.
InJ
dd .. I
Jffi JJ Jm JJ
-< ...
i
...
i s given by
A d .
b'
7. CONCLUSION
8. REFERENCES
Kateri Guertin
Applied Earth Sciences Department
Stanford University
ABSTRACT
I. INTRODUCTION
G. Veri), et al. (eds.). Geostatistics .fvr Natural Resources Characterization. Part 1. 245-260.
o 1984 by D. Reidel Publishing Company.
246 K.GUERTIN
conditional unbiasedness
E{(Z(x)-Z(x)*)/Z(x)* ~ zo} = 0 for all Zo,
conditional estimation variance
E{(Z(x)-Z(x)*)z/Z(x)* ~ zo} to be minimized, for all ZOo
K(Z(x)*) = Z(x)**,
to be applied to any initial linear estimator Z(x)* to turn it
into a conditionally unbiased estimator Z(X)**. Such a corrected
estimator would insure accurate prediction of recoverable
reserves from an ore deposit [2): Z(X)** would be such that
Zv
Regression curve:
E{Zv/Zv*=zo}=Kv(zo)lz o
m+--------./..-~~
Zv
Regression curve:
E{Zv/Zv**=zoJ=Lv(zo)=zo
m+------,~---,,,
~-==------I-~~------~
m" Zo
Zv"
CII
CII
E{Z(x)/Zv} = Zv.
Consider also an isofactorial binormal distribution model for the
corresponding' couple (Y(x),B), with a coeffic.ient of correlation
'r'. Then, a model 'of the distribution of the block grade Zv
CORRECTING CONDITIONAL BIAS 251
00
Zv = ~ Tn rn 7In(B). ( 3)
n=O
02(V/O) =~ Tn z r 2n , (4)
n=1
00
00 00
00
where B* = ~·'(Zv*).
IV
V.
IV
2.4 1.0
~
2 .2
l 0.9
2 .0
::> 0.8
u
,.., 1.8 Zv ' - '" '3
Z. " • 0
u
S 0.7 r lv' - •
1.6 00 0 ~
~ -,.., Zv" - 0
~ 0° ~ + +
0
1.4
0
J .
,. 00 •+ >
• • •
N 1.2 N 0 .5
'- o • '-
>
N " 0 •
~ ~ ,'",. + •
1.0 o .+
LL
tS>
o ..
'1-+ + l N
>
°t
0.4 ~
~
0
0 .8 Ot N
w • 0
:J 0 . 3r- . .
--' 1 LL + 0
<: 0 .6 tS> +
> o 0 0 0 0 0
w o 0 0 0
w 0.4 :J D.?
19 --'
r+ · · o a
<: <: o 0
Q:: > 0.1 0
L'.l 0.2
".> W
<: 19
0.0 L , .I " , I. " I " I., , ,I , " I " , I , , , I , " I , <: 0.0
Q::
0
J~
Cr .::> 0.4 U.5 0.8 1 1.;> 1. 4 1.6 1.8
, ,\ LRA"F:' VAL l iE' I'F {lv ' fly ' • 1~) In %ru
2
,~
2 .2 2.4 w
>
0 0 .2 0.4 0 . 5 O. B I I.;> 1.4 1.6
CUT0FF GRADE ZD APPLIED 0N Zv' (%Cu)
1.8 2
<:
STAt. -~·, r ~ ~ : 1" ,1 ~1',r'- Cll.ISTrRF() DATA Pf1INTC:: ST ANHIRD 2 : 135 N01\·CLUSTERED DATA F01t\TS
Fig.3 Cumulntive conditional binS given a Fig.4 Cumulative conditional squared error
cutoff Zo; Zv" holds for the initial given a cutoff Zo; Zv" holds for the
~
estimator Zv* and its corrected estima- initial estimator Zv* and its corrected C'l
tor Zv**. Each symbol '+' or '0' holds estimator Zv**. c:
!Tj
Two sets of data points are taken from the deposit (with and
without clusters). From each data set. 500 mining blocks of equal
size v are estimated by the polygonal method: each block receives
the grade of the data point closest to its center. Then the two
sets of linear block grade estimates are corrected for their con-
ditional bias using the model described above.
2.25 +
'3 ... '3 •
u u ..
2.00 t '" ;>.00
'" •
... + .
~ >
> ~ ... ... + N .
N \.75 ...... ... ++
+ , 1.75
W W +
t :t.. ++ f 0
0
... I
\.50 ..... .......t .... n:: \.50
..,:t:-4+ ... "'t-""
.. ++ ... +
+r+ . . .
n::
"" .... ""
+ ... + ++ ++ +++ ++
1.25 .. + ++
++
'""-u
1.25 ~ \lR- 31
'"u"-
IS) >'. + ... CSl
++ -'
.+. ......
-' 1.00 m
m I +t 1.00 ...
~ + '+
w w
:J :J
0.75 n:: 0.75
t-
"" 0
o 0.50 m -.5815: 0' -.261
0 0 0.50 m -.5816: ,,' -.261
If) If)
m·•• 5589 : a2--. 32~ m···.5677: 0"" -. 237
Co" .CooF .-.834
Fig.7 Scattergram of Zv vs Zv*; notice the fig.S Scattergram of Zv vs Zv**; notice the
strong conditional bias and the unbal - reduced conditional bias and the more
anced areas of misclassified units: WR = balanced areas of misclassified units:
25 and WS = 64 for 20 = 0.6 %Cu. WR = 31 and WS = 41 for 20 = 0.6 %Cu. ~
C'l
C
to1
'...,Z"
CORRECTING CONDITIONAL BIAS 257
fig.9 Correction of Zv* by Zv**. Notice that fig.10 Correction of Zv* by Zv**r. The '+'
the unique correction function Zv** = holds for the blocks located within the
Kv(Zv*) is strictly increasing. non-clustered areas (nca), the '0' for
those located within the clustered
areas (ca). Notice that both regional :><:
correction functions Kvrnca and Kvrca C'l
C
are strictly increasing. m
::<>
-l
Z
CORRECTING CONDITIONAL BIAS 259
Using again the 150 clustered data points and the initial
poly estimator Zv*, such a regional corrected estimator Zv**r is
obtained. In terms of conditional bias and squared error, Zv**r
is practically equivalent to Zv** (previously defined given a
unique overall average covariance, figures 5 and 6); the plots of
the global correction function Kv(Zv*) (figure 9) and of the two
regional correction functions Kvr(Zv*) (figure 10) are also very
much alike.
VI. CONCLUSIONS
REfERENCES
A.G. Journel
Abstract
Keywords
Introduction
G. Verly et al. reds.;, Geostatistics for Natural Resources Characterization, Part 1,261-270.
© 1984 by D. Reidel Publishing Company.
262 A. G. JOURNEL
Characteristics of deviation
(3)
Estimation of F(z/(N»
Prob*{Z(x) IE 1z , z ,J/Z(xa~";Z
q q ..:
, "'€l"')~ ~
•
C(-q (4)
with: q' > q Eo [0, 1J , an~ z~ being the q-quantile of the esti-
mated conditional cdf Fx (z./on) .
~
If the estimated conditional cdf's F (z/(N» are
asymmetric and account for reversal of sk~wness, so will the
prediction intervals ]z , Zq') .
q
Expressions for the mAD
with:
JU z.dF(z)
o
= zF(z)l u0 - J,uF(z)dz
0
= uF(u)-J~F(z)dz
0
The integral J:
F(z)dz has no limit when u
to introduce the following "cdf integral":
+ +00, so it is better
(7)
A(u) = J~[l-F(z)]dz u[l-F(u)] + m(u) .F(u)
Generalization
L(z.u..)
Remarks
Conclusions
References
• u
Prob {Y $ y} .. A(y)/A(-) • 1.
m
• A(y) (3)
E lz 2} .. 2E {Z} • E {y }
A(y ) .. pm 05)
p
270 A. G. JOURNEL
A (u) •
n
f0 z n-l [I-F(z)] dz, with: A1(u) = A(u) (16)
EIF-AQUITAINE (PRODOCTION)
ABSTRACT
Hhen processing seismic times, it is of upmost impor-
tance to take into consideration the action of
faulting, to allow both a correct estimation and
display of the data. Kriging with a fault throw as an
imposed drift allow to include a very precise model of
faults to the standard drift functions. After
examining in detail an academic example in I-D,the
routine 2-D practical situation is studied and 2 real
examples are shown. KE~mRDS : Seismic data, kriging
with faults, contouring.
Introduction
Conventional seismic data appear at a given stage and
after convenient processing, as two-ways times to a
given horizon, q~signed to a series of shot points
along.a profile. Next step in data processinq will
consist in plotting the isochrone curves relative to a
given area, based on the time values of a certain
number of profile (fig.!).
271
G. Verly et al. (eds.). Geustatistics fur Natural Resuurces Characterizatiun, Part 1, 271-294.
o 1984 by D. Reidel Publishing Cumpany.
272 A. MARECfAL
I
I X.CI
\~>< {~ = ~t )::tI,...{,.-~
,: " J. f; e
¥
~ ~
Ao~3o AA;: ~
~2- ::. 5f't4 --~':Z.f- - f- <;n-~-A.
~
-- (A-l) 3" -3rt'-
~-p_A
A3~
These results are in agreement with the usual results
found with a a linear variogram and linear trend
1""(
o 2. 4 , g
r--. "-
4. The two sample values ;5~,;. and~", when
corrected of the throw cOlncide with the extrapo-
lation from the last undisturbed point zp.
Indeed, the throw function PCk) for each sample
k, k=p+l,n has the following value :
~k-;:C 5~ - D(J.t)
D(k}~ 3~-3p" {~-f)~i3. - (k_!) ~:1:1t;-
3" ;. 3-r- 4- ("-f' 3 ?!?. ft
3p", ~ >ft--f ¥
Practical case of a succession of faults
\
'\,
I
I
25011
~---r-/
I Actual Profile I
I
I
Fault
f-_
2~~--~----4-----~-- __________ ~ __ ~4- __ ~~ __ ~ ____ ~_
5000
Fig: 5 EJW Profile no 4 of data of F 19 11.
-
10.
~4---------------~----------------~---------------------
282 A.MARECHAL
Theoretical solution
{(~'}1c: qo~d~,})-t~.{)(,~"("t~J44,~~()((~)
2. by digitizing the shape of a throw function, of
unit vertical shift and given lenght of influence
along a given fault segment L. The unknown throw
will then be considered as proportional to this
unit throw function 1L()(t~) :
+(~l d-):::' a- TL ('Xl~)
Practical solution
As each individual fault in 2-D is described with a
much larger number of parameters than in the 1-0 case,
it is easy to realize that the strict theoretical
approach is usually difficult to implement.One can
th~nk, however, that such solut~on could be
implemented in a graphical way similar to the one
applied in C.A.D. when defining the fault traces in
a given area, the interpreter either would draw a
model of the fault movement later digitized by the
system, or he could pick up one particular model in an
already existing catalogue. At the end of this step,
the computer would resolve the above kriging system.
In order to test the practicality of the use of a
throw function in kriging, the strict theoretical
solution given above was abandoned and replaced by an
approximation which can be operated on the existing
krig ing sol ftwares (KRIGEPACK or BLUEPACK).
o
Fig. 7 Radian Data : Standard krigepack processing considering
faults as barriers
Fig. 13 Sattlegger data: Final grid with fault throw function contoure.-j
directly by a contouring with fault program.
110·
In.
ACKNOWLEDGEMENTS
BIBLIOGRAPHY
Donald E. Myers
Department of Mathematics
University of Arizona
Tucson, Arizona 85721
Abstract
INTRODUCTION
295
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 295-305.
© 1984 by D. Reidel Publishing Company.
296 D. E.MYERS
The "under sampled" case is the one that has received the
most attention in mining applications but only for a few
variables and a few sample locations. Section three describes
a simple algorithm for obtaining the system of equations ir.
matrix form for any "under-sampled" problem from the general
form of Co-Kriging.
NOTATION
-L LY(X.-X.)AiA.;> 0
~ J J
(1)
n n T- T-
TrE L
L i=l fi Z(x.) Z(X.)f j ]
j=l ~ J
n n T
Tr[ L Lfi K(x. -x. )f j ] ) 0 (3)
i=l j=l ~ J
for all f l' .•• , f with
n
Al A2 ••• ,A n
ss' ss' ss
are diagonal entries in f 1, "', fn respectively.
A!t' "', A;t are also diagonal entries. Let
298 D. E.MYERS
A
s
[A 1
ss'
... , Anss 1 (8)
Bt
1
[A tt , ... , Antt 1 (9)
K = (10)
uv
k' (x 1-x ) ... k (x'-x)
uv n uv n n
then (16) can be written
A KAT + A K B T B K B T 0 (11)
s ss s s st t + t tt t ~
and (17) as
A KAT ~ 0
s ss s
BtKttB t T~ a (14)
BLOCK CO-KRIGING
Stationary/Covariance case
C(h)
C l(h) •.. C (h)
m mm
then
E[Z(x 1. )]T[ZV] = = l J C(x.-x)dx = C(x. ,V) (16)
V V 1 1
and
- T- 1-
E[ZV] (Z(x.)] = -V J C(x-x.)dx = C(v,x.) (17)
J v J J
- T-
E[ZV) [ZV] =
1
-z
J J -C(x-y)dxdy = C(V,V) (18)
V VV
The Kriging equations become
n n
L C(x.x.)r.+~=C(xi'V), L r. I (19)
j=l 1 J J j=l J
with Kriging variance
(20)
(21)
and becomes
(22)
2
0CK = Tr Im y(x.,V) f.- Tr]J - Tr
1 -
L f f y(x-y)dxdy (26)
1 1 1 V VV
-4 f f y(x-y)dxdy (27)
V V V
2
are not required otherwise the last term in o CK is simply
n 1
I 2" f f Y i i (x-y)dxdy (28)
i=1 V· VV
1
r:GJ :
V is the square with vertices A,B,C,D.
Data Z1(0),
Z2(0), Z2(A), Z2(B), Z2(C)
Z2(D)
It is desired to
estimate Zv where
CO-KRIGING - NEW DEVELOPMENTS 301
y *(so ,SR) I rr l
!
-*
y *(SR,SR) 1* I r2
1* 0 L~
where
y12(80,SR)1
~
Y*(SO,SR)
(30)
Y22(SO,SR)J
~ y,:(8 ,8J
Y ;(SR,SR)
(31)
R
~ y,:(v.J
y*(SO,SR)
(32)
[: :] [: :]
*
1* I (33)
,
I .. I 0 ~ I
302 D.E.MYERS
jo jo
A. l' •.. , >... are all zeros.
10 lOrn
A
_
Let r.O
J
be the modified
_
r J.0 Similarly let
y*(x.O-x) be the same as y(x.O-x) except that in the
J p J _*p
iO row all entries are zeros and y (x -X. O) is the same as
_ p J
y(xp-x. O) except that all entries in the iO column are all
J _ - T
zeros, that is, y*(x -x. O) = y*(x.O-x ) Let I* be an
P J J p
identity matrix except that the iO row is all zeros and
(I*)T = I* •
j* j 0 J J
+ Tr [ I I r. Ty (x -x.) r j
1* j 0 j* j 0 1 n J
T-
I
A
+ r i y(x.-x· o ) r·o
i* jO 1 J J
T-
+ I r J· o y (xJ·O-x i ) r J.
j* jO
(35)
=I
~
L I r _ + I r J. O (36)
j*jO J
Consider the following identities
(i) (37)
(ii)
(iii)
(iv)
(v)
(vi)
K = (38)
.r (xn -xn )
E (39)
(40)
(41)
i.e.
h' :] n [:J -
=
-
K r + E 1.1 = K0 ET r =I
(42)
( 43)
This not only reduces the size of the system but also avoids
the possibility that (43) is ill-conditioned.
AX = B (47)
f. (X)
1
X- ~
ai
[A.1 < X, A.>
1
- A.B.]
1 1
(48)
(49)
By writing
B = [B ••• B ]
1 P
CO-KRIGING - NEW DEVELOPMENTS 305
References
A.G. Journel
Abstract
Keywords
Introduction
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1,307-335.
© 1984 by D. Reidel Publishing Company.
308 A. G. JOURNEL
with:
1
z (x) Iv(x) z(u) du, being the mean grade of the smu v(x),
v v
- quantity of metal
recovery factor
(2)
(3)
Fv(h- l(z)/(N»
I (x;z) = 1 - J (Xjz)
_ {I,
if Z/x) ~ z
, then:
(9)
v v
0, if not
00)
F (z/(N»
v
LN = ~ + 1: 111 Aa Z(x a) + \ • 1 J C HN
THE PLACE OF NON-PARAMETRIC GEOSTATISTICS 313
~: Space of all
measurable
functions of
Z(x~), ~ e (N)
Fv(z/(N»:
projection of
IvCx;z) onto HN
[FvCz/CN»Jt ;
projection ~f
IvCx;z) onto
LN C HN
0,
.
It thus appears that the determination of the projection:
..
~(Z(xa), at!::(N» = LIv(x;Z)JtK :: [Fv(Z/(N»]LK ,of Iv(x;z)
Some examples
r 1, if Z(x) ~ z
~k(Z(x» = I(x;z) = to, if otherwise
t! ! ( 13)
SCUn, I) + 1: Va(z) F(Z(x a » + Aol (z). I(xrv;z)
u.
+ ). (z)}
0
<=1 <'/"1
""
the ¢n being known coefficients. The estimators [hn(Y(x»]SK
are recombined into the DK estimator of Z(x):
00
It thus appears that all three, IK, PK, DK, are service-
variable approaches, considering different transforms t(za) of
the initial N data. The question is then to determine which
transform is the most adequate for the particular goal at h· .,1,
which is the estimation of the conditional cdf Fv(z/(N) thr·.·agh
that of the indicator R.V. Iv(x;z), cf. relations (9) and (lO)
and Figure 1.
(16)
Z (x) =
v v
I vex) Z(u)du.
with: Y(h)
l..e. ro
t( z) . dF (z) <
v-
r
o
t(z)dF(z)
VadZ(x) / Z(xw = z 0( ,
aC£(N)} - Y(v,v/(N».
322 A. G. JOURNEL
(19)
Remarks
i(x . ~ =
e(J)
i I, if ZCt. < z,
0, othe(w{.'>e
- • a{ fiil (N)
Unbiasedness conditions
It comes:
N
F(z). I \JCt.(Z) + "o(z)
<r1
THE PLACE OF NON-PARAMETRIC GEOSTATISTICS 325
1 " o (z) o
t v (z)
a
" (z) =
o
= 0,
(1- L
for all at::(N), and all z
N
0.=1
"a(Z)].F(Z)
(25)
The PK system
~
~l
AS(z) Cr(xa-xS;z) + !
~l
'VS(z) Cur(xa-xS;z) = Cr(x-xa;z)
for all a E (N)
- rndicator covariance:
(28)
Cr(h;z) = EtI(x+h;z) I(x;z)3 - [F(z)r = S2(z) - Yr(h;z)
(30)
E(Un(x+h). r(x;z)3 - ~F(z) Cur(O;z)-Yur(h;z) 20
with:
2
-~ s (z)
2
z F (z)
since: Jo F(u) .dF(u) = - -
2
A :(Z(".~) >z
hex) , z
Z(lt+h) B: (Z(x+h) ~ z
tz(x) > z
c: {Z(X+h) oS- Z
Z(x) ~ Z
D: lZ(x+h) > Z
tz(x) > Z
E: rZ(x+h) " z
\Z(x) ~ z,
F: Z(x+h) '"... Z
G: Z(x+h) "i'" Z I
Z(x)
o
Z f-+--+----l-
Remarks
-
213 • Cur(h;z) ~ 0, yielding:
S(u)
(34)
Pz(h) • F(min (z,z'» + [l-Pz(h)] • F(z) • F(z'),
with:
YZ(h)
1 - ~,.------ being the correlogram of z(x) •
Var {Z(x)}
Conclusions
References
Bruce M. Davis
The deposit and sampling present problems which may cause diffi-
culties in other linear and nonlinear estimation procedures.
These problems are either mitigated or avoided by using I.K.
INTRODUCTION
G. Verly et·al. (eds.) , Geostatistics for Natural Resources Characterization, Part 1,337-348.
© 1984 by D. Reidel Publishing Company.
338 B.M. DAVIS
SUPPORT DIFFERENCES
*
*
*
V 27+ ** * *
0 * **
1 *
u
m
-*
e 18+**
- *
--* *
9+
- *3"** 2 3 *
0+**23* * 2x ________ 3x ________ 4x ________ 5x
222 1x ________
D________
Grade Au
Figure 1 - Scatterplot of Sample Volume Versus Grade
30+
F I
r 20+ *
.
e I
q I
U I
e
n
_I
10+ 1
*
I
C _1*1
Y -*** I
-**.»* I
*
_1****** * *
_*******ff*** 1* ** I
0+---------1x--------2x--------3x-------~4x--------5x
Grade Au
Figure 2 - Histogram of Au Grade
+0* 00000000111112222222222233344
+0. 55566788889
1* 01
1. 6799
2* 0
2.
3*
3.
4*
4. 6
Figure 3 shows the normal score plot (6) of the natural log-
arithm of grade. Also shown is the correlation between the nor-
mal scores and the log-grade values which is a measure of the
INDICATOR KRIGING AS APPLIED TO AN ALLUVIAL GOLD DEPOSIT 341
+ *
« ** •
*
III
*
L +
n *2**
* 22
G **
r 12*
a + *2222
d **2**
I
e
A 1*1*
u + *
**
+
*
+---------+---------+---------+---------+---------+
-2.5 -1.5 -0.5 0.5 1.5 2.5
Normal ScoreS
Figure 3 - Normal Frequency Plot of the Natural
Logarithm of Grade
One method proposed to the author for dealing with the ano-
malous values was rejection of the data above an arbitrary cut-
off. Another method involved the adjustment of all values above
an arbitrary cutoff (e.g., 2.5 standard deviations above the mean)
to the mean value of the histogram class containing the grades
greater than the cutoff value. After either rejection or trim-
ming it was suggested that a natural log transformation be made
to deal with the skewed distribution of the data.
GRADE ESTIMATION
~'
if z(~) > z
c
i c (x)
- =
1, if z (~) < z
c (1)
18+
*
•
• *•
15+
*. * *
* *
*** *
* *
T 12+
*
**
* *
h *** *
2* * 2
c
k -*
2 * *
n 9+**
e
s *
s - 2
6+* *
*
**
3+
+--------1x--------2x--------3x--------4x--------5x
Grade Au
Figure 4 - Scatterplot of Gravel Thickness Versus Grade
Indicator Variograms
.40
.35
.30
"'~.
.25
~-'~~~------~~'
~
~ .20
,
.,,10'
'\(
.• 5 ~ .
.• 0
.05
Y(h) 0.10 + 0.14Sph(1100)
, ~.
oO''--..L--,.20"::0--'----40,I:0--'---::6:':00:-'---::8:':00:-'~10=0:"0-'--"12=0:-0--'---"1:':40::::0-'--;;:;1600
LAG [FEETl
. ~o
.35
,-~ .
. 30
,""
"'~ . ,~'1-'
• 25
'"
~
~ .20
i3
.15
.10
.05
Y(h) = 0.10 +0.14Sph(250)
o ,L[
o
--'~"::[
200
--'-.,.L---'---,-L--'__L---'---,L:-"c"..-,L!---"'--.• ..L!__' - - '
400 600 800 1000 1200 1400 1600
LAG IFEEjl
sistent with the sill values of the fitted models for each of
the three indicator variograms.
¢*(x
-0
;Z c ) = -A'i
- + (1 - -A'e)F*(z
- c)
(2)
100+
90+
80+
F 70+
r
e 60+
q
u 50+
e
n 40+
C
y 30+
20+
10+
0+-------ZC1------------ZC2-----------------ZC3-----
Grade Au
Figure 6 - Estimated Cumulative Frequency Histogram
of Grade Inside the Pit
346 B.M.DAVIS
¢**(~o;zcl) Max{¢*(~o;zcl)'O}
¢**(_x o ;zc2) Max{¢**(x'z
-0' c l) ' ¢*(x'z
-0' c 2 )}
¢**(x;z ) Max{¢**(x;z ),¢*(x;z )}
-0 c3 -0 c2 -0 c3
¢***(x'z ) = Min{¢**(x'z ) l}
-0' c3 -0' c3 '
(3)
x {¢***(p'z
, c3 ) - ¢**(p'z
' c2 )}
(4)
where zci, i = 1-3 are the three cutoff values and Zm is the
average value of the last class of the raw grade histogram. The
value Zm is used instead of the usual (zmax + zc3)/2 to resist
the anomalous value zmax'
INDICATOR KRIGING AS APPLIED TO AN ALLUVIAL GOLD DEPOSIT 347
CROSS-VALIDATION
*
60+
o
50+
*
F 40+
r o
e
q
u 30+
e
n
c
y 20+
*
0
0+-------ZC1------------ZC2-----------------ZC3-----
Grade Au
Figure 7 - Comparison of True and Estimated Proportions
348 B.M.DAVIS
CONCLUSIONS
REFERENCES
I.C. Lemmer
ABSTRACT
I. PRELIMINARIES
G. Verly et al. (eds.), Geostatistics jar Natural Resources Characterization. Part 1. 349-364.
© 1984 by D. Reidel Publishing Company.
350 I.e. LEMMER
1, if z(~) < z
i(z(~);z) i (~; z) (1. 1)
0, if z(~) > z
and their mean value over panel area A (within deposit D),
¢(A; z) = 1/ A Ji (~; z)
A
dx (1. 2)
t (A; z) 1 - ¢(A;z)
(1. 3)
q (A; z) d¢(A;u) du
dU
(1 .5)
a a
simple kriging (SK) (J. 7a)
¢*(A;z)
I k (z).i(x ;z),
a --a
ordinary kriging (OK) ( J. 7b)
a
J=
vS(z).PI(has;z)
S=I
I/A f. PI(hao;z)dxo = PI(a;z)
N A
OK: I AS(z).PI(haS;z) a = I, ... , N (J. 8)
B=I
(J. 9)
2. INDICATOR COVARIANCE
faS(zo'zS) = f(za)·f(zS) I
Tn(a,S) Xn(Za) Xn(zS) (2.1)
n=o
to the grades. We recall here that the orthogonal polynomials Xn ,
having the marginal density f as weight function over their
orthogonality interval, are uniquely determined by this weight
function, provided that certain rather general integrability
conditions are met (3). From this fact it follows that the first
two polynomials are
The actual values of the a n 's (some or all of which may be zero as
ESTIMATING LOCAL RECOVERABLE RESERVES VIA INDICATOR KRIGING 353
One can hope that breaking off 01 at the first term will be
good when p(h) is small. Anticipating such an approximation
comes from practical experience. One knows empirically that, due
to the nugget effect, p(h) exhibits a very sharp drop from unity
as h increases. Thus, excepting h = 0, this empirical circumstance
suggests the linear form (2.8) as a realistic expression for
0I(h;z). As with the full series, this approximation still
contains the destructurization effect in the vanishing of the
cofactor of p(h) at large z.
'.•r~'----------------,
flr hO, 9
o o fllI . OT qI
o
o !'tl hO,S 0; o FU hQ)
o
. RtI . 0.5
' .0 1- - - - a Rr t- O,'
o 0.' o
• I
O,t
••
.
0
q. o ~,
0.' o~ 'p
, V(h)
· OJ' 0 0
0
0.1
10
'. '00 '\0
-0,4
3. INDICATOR KRIGING
Cutoff z Weig~t
where 9 SK weights of F
F(z)=0.3 7.3 8.1 7.3 4.8 13.7 13.8 14.1 13.7 8.1 9.2
F(z)=0.5 7.1 8.0 7.1 4.3 14.5 14.6 15.0 14.5 8.0 6.9
F(z)=0.7 7.2 8.1 7.2 4.5 14.3 14.3 14.7 14.3 8.1 7.5
F(z)=0.9 5.6 6.0 5.6 1.7 16.9 15.7 19.3 16.9 6.0 6.3
F(z)=0.98 -0.1 -0.0 -0.1 0.0 3.5 2.1 5.8 3.5 -0.0 85.2
9 OK weights
F(z)=0.3 8.2 9.3 8.2 6.0 14.6 15.0 14.7 14.6 9.3
F(z)=0.5 7.8 9.0 7.8 5.3 15.2 15.5 15.4 15.2 9.0
F(z)=0.7 7.9 9. I 7.9 5.5 15.0 15.4 15.2 15.0 9.1
F(z)",0.9 6.2 7.0 6.2 2.7 17.5 16.7 19.5 17.5 7.0
F(z)=0.98 9.4 9.5 9.4 9.6 13.0 11.6 15.1 13.0 9.5
Grade(OK) 4.9 5.5 4.9 1.0 19.0 17.4 22.7 19.0 5.5
Note the constancy of the weights, and, since they are a hundred
356 I.e. LEMMER
Let us suppose next that the A's in eq. (3.1) have been
determined by OK. The next question is: with respect to estimating
what quantity is the estimator ie*(~;z) unbiased? To answer this
question, consider ie* and ie as particular realizations of
random functions Ic* and Ie in the usual way. Then if F(z) is
the distribution of actual grades,
= J F(z(~»).fG(z(~) - z) dz(~)
or
E{IG(~;z)} = f E{I(~;z~)}.fG(z~ - z) dz~ (3.6)
(3.7)
(3.8)
(3 . 9)
Then
(3. 10)
>-
(J
z
W
:l
ow
0:
U.
% CU
0,'
1,0 _ ____ _ ___ _ o-
f
~ "
Cutoff z
where 9 OK weights
F(z)=0.3 5.9 6.9 5.9 2.3 17.7 17. I 19.8 17.7 6.9
F(z)=0.5 5.5 6.4 5.5 1.7 18.2 17.2 20.9 18.2 6.4 I,
F(z)=0.7 5.3 6. I 5.3 1.4 18.5 17.3 21.5 18.5 6.1
F(z)=0.9 5.8 6.7 5.8 2.1 17.9 17. I 20.? 17.9 6.7
F(z)=0.98 6.9 8. I 6.9 3.9 16.3 16.3 17.2 16.3 8.1
Grade(OK) 4.9 5.5 4.9 1.0 19.0 17.4 22.7 19.0 5.5
_ "TRUE _ TRUE
04 04
---- IK ---- IK
03 03
02 02
"
1:
01 :'·'--1 01
00
CUTOFF GRADE , OS 10 15 CUTOFF GRAOE ,
_ TRUE _TRUE
04 ---- GIK 04 ---- GIK
03 03
02 02
01 0'
10 15 CUTOFF GR ... OE , OS 10 15
_ TRUE _ TRUE
---- MG 04 - --- MG
03 03
02 02
"-
01 01
00
05 10 1.5 CUTOFF GRADE , 05 10 15 (UTOF~ GR ... OE ,
Figure 4a
Figures 4a and 4b (next page) show the true and estimated incre-
ments in CP(A; z) c:orresponding to 18 chosen cutoff grades, comparing
lK, GlK and MG for four panels or 'blocks' of Stanford II.
362 I. C. LEMM ER
b .(It,ll 1t11t,1I
,LOCK 3/16 6
BLOCK 5121
OS 05
TRUE
_TRUE IK
04 ---- IK 04
OJ OJ
02
'"
01
00
CUTOfF GRADE • 1S
CUTOFf GRADE Z
_TRUE _ TRUE
04 ---- GIK 04 ---- GIK
OJ OJ
02 02
01 01
00 00
OS CUTOff GRAOE Z CUTOFF GRADE.
03 03
02 02
01 01
00 00
1S [U TOFF GRAOE •
Figure 4b
In each estimation only the 9 closest drillholes were used. A
'spike' for zero cutoff grade (representing waste in the panel)
is shown just to the left of zero cutoff grade.
ESTIMATING LOCAL RECOVERABLE RESERVES VIA INDICATOR KRIGING 363
:! r'
0,1
0.1
0,1
,~-----~-----~----__4
',' C,' I.e
ACKNOWLEDGMENTS
REFERENCES
APPENDIX
(A. I)
JEFF SULLIVAN
ABSTRACT
ACKNOWLEDGEMENTS
The author would like to thank the freeport Gold Company for
graciously supplying data from their Jerritt Canyon deposit.
365
G. Verly et al. (eds.), Geustatistics for Natural Resuurces Characterization, Part ],365-384.
© 1984 by D. Reidel Publishing Cumpany.
366 J. SULLIVAN
INTRODUCTION
i (x; z) = fH Z (x)
o if Z (x) ) z
~ z
i(x,z)
o z (x) z
ri(h,z)= .5 E[(I(x+h,z)-I-(x,z»21.
n
~>"o.ri(xo. - XB'Z) = 'li(XB'V,Z) J3 =1 to n
1
n
a 2 ik = ~>"o.ri(Xo.,V,z) - ~i(V,V,Z)
1
n n 1\
.p'1r(V,z) = l~o.ri i(xo.,z) + (l-lAo.)f(Z)
1 1
1\
where Fez) is the mean of all indicator data for cutoff z. The
recovered tonnage at cutoff z is simply:
(1-~*(V,z»(tonnage factor).
....)
-:f
~ (V , z)
c
z
Fi +,
subject to
> Fi for = " n-l
F1 >0
Fn (1.
where Fi is the final solution (dash~d line)
Ii iS'the initial solution given by IK
(11«V,z»
and Wi is a weight given to each cutoff. if
Wi = l/n for all i the the objective
amounts to minimizing the average squared
deviation.
o 2
Z (X)
The shaded area is the value o. the indicator variogram for
cutoff z This value, although v~ry useful does not give any
indication of how the mass is distributed within the shaded area.
For instance, a bivariate distribution which is closely packed
about the line Z(x)= Z(x+h) and a distribution which is loosely
packed could easily have identical indicator variogram values.
Compare distributions A and B below.
o !
o " - - - . . L2- - - - - - : :Z-;( x)
A B
n n
~*(v,z) = ~A~(X~ z) + LV~U(x~)
1 1
n n
~~~Pi(X~ - x~,z) + rV~pui(X« - x~,z) + ~1 = Pi(X~'V,z)
1 1
n n
r~~pui(X« - x~,z) + rv~pu(x« - x~,z) + ~z = Pui(X~,V,z)
1 1
n
r~o. =
1
n
LVo. = 0
1
372 J. SULLIVAN
n
o 2 pk =p;ev,V,z) + ~, - ~~«puieV'Xa,Z)
1
COMPARISON OF IK AND PK
varCfCV,z» - varCf*Cv,z».
02ik (appendix).
02pk (appendix).
~:t
0 .3
mean =.5
variance=.25
0.2 200 data
0 .1
20 40 60 80 100
Indicator Variogram Stanford 2c
median cutoff
0.4 ,,~~~-r-r~~~~~~~~~~-r~
~I
0.3
mean =.5
variance=.25
0.2 200 data
0 .1
20 40 60 80 100
1.0
0 .8 true
c: 0 .6 IK
ro _.- PK
Q)
e 0 .4
0 .2
0 .0
0 0 .25 0.5 0.75 1.25 1.5 1. 75 Ceu toff)
Mean recovered tonnage as a function of cutoff.
This plot shows that both IK and PK
are nearly unbiased over the entire
range of cutoffs.
~ 0.10 , true
c:
ro I ...-
...- .- '. IK
.~
0.9
<lJ
0'10 .8
C\l
c:::
C::: OJ
o
E-<
'0 0 .5
<lJ
~ u.5
:>
8 0.4
<lJ
0:; 0 .3
Q)
::l
1-1 0.2
E-<
0 .1
0.9
Q) 0 .8
0'1
g 0. 7
C\l
o
E-i 0.5
'0
Q)
1-1 0.5
Q)
:>
0 0.4
U
Q)
0:; 0 .3
Q)
::I 0.2
1-1
E-i
0.1
After examining the data base and checking that there were
no obvious data entry errors, variography was performed. The
variogram of grade (fig 4) shows no structure beyond 100ft.,
however, the twin hole information indicates a fairly low nugget
effect. Given this variogram model; the spacing of the data; and
the fact that the smu's are much smaller than the size of the
panel, it would be impossible to accurately krige the grade of
the s'elec~ive mining units (smu's) due to the smoothing of
ordinary kriging. Hence. the necessity of a conditional
distribution approach such as PK which will ultimately determine
the distribution of smu's within each panel. At this stage.
however~ only point recovery functions will be estimated since
the true smu recovery functions are not available for comparison.
Z Variogram
't z 0.06
mean =.105
(587) variance=.044
0.04 193 data
0.02
"'t1 0.20
(307)
(587)
0.15
mean =.71
variance=.206
193 data
0.10
0.05
0.00
0 100 200 300 400 500 600
Conditional Variogram(70th percentile)
- 0.15 ~-r-~-r-~-r-~-r-"'--r-,....,-'--"""'-'--"""'--'--'
(587)
-0.10
mean =.254
C . (0,z)=-.103
Ul
-0.05
193 data
0.00 0
600
CONCLUSION
0.9
.'
0,8
<1.1
C'
100.7
c
c
~ O,6
"
'0
~ O.5
ClI
>
8 0.4
(])
0::
<1.1 0.3
...
~
E-< 0.2
0 .1
0 .0 It! I! I,
I t I 1
Summary Statistics
mean variance
true .791 .0573
estimated .811 .0604
/'
1.0 I' iii i i 'l"'T'
,
, .
,
0.9
/.
0.8
ill
cr>
to 0.7 ' ,
c
t:
0
E-< 0.0
'0
.... 0.5
ill
ill
>
0
u
0.4
Q)
c:t:
Q)
0.3
+,'
:J
....
E-< 0.2
0.1
Summary Statistics
mean variance
true .499 .0890
estimated .512 .0778
:> 0.5
0 "
.)
l>
rc:
0.4 '.
1)
0.3
...
:l
c....
0.2
• 1
0.1
. ..
0.0 L.Li......J....J....l....l....1....I....I.....J....~: , ' I' ,,I ! • , I, ' , l..w..L-Li..-J..-.L..LL...J
o ll.l c.? 1.1 0.4 o.t:', 0.", 0.7 0.8 0.9
PK Estimated Recovered Tonnage
PK Estimated Recovered Tonnage
Summary Statistics
mean variance
true .201 .0454
estimated .219 .0513
APPEMDIX
Smoothing of IK and PK
n
SZ(z) = 1/k l(~(Vk>z) - "l(V,Z»2
1
n
"lev,z) = l/k l~(V,z) = ~eD,z)
1
n
SZ(z) = 1/k l?2(V K,Z) - ?z(O,z)
1
n
E(Sz(z» = l/k lE(~2(VK'Z) - E(?Z(O,z»
1
E(~2(V,Z» = l/Vz ffE(I(x,z)I(x',z» dx dx'
=C(V,V,z) + f2(Z)
E(?2eO,z» =Ci(O,D,z) + f2(z)
ECSkCz» = Ci(V,V,Z) - CiCD,D,z)
n
~*(V,Z) = l~«(z)[i(x«,z) - fez») + F(Z)
1
the sk system is:
n
l~«Ci(X« - XIPZ) =Ci(A,X,8'Z) B =1 to n
1
n
a2iK = Ci(V,V,Z) - l~«(Z)Ci(V,X«'z)
1
n
S*2(z) = l/k l(~*(V,z) - ~*(D,z»2
1
384 J. SULLIVAN
n
= l/k r(?*(V,z) - f(zj)2
1
E(S*z(z» = l/k ~t~«(z)~p(z) Cj(xa - xp,z)
= l/k ~~~aCZ)Cj(V,x«,z)
= l/k tCCj(V,V,z) - aZjk(V,Z»
= C(V~V,z) - azikeV,Z)
REFERENCES
Alain MARECHAL
ELF-AQUITAINE (PRODUCTION)
ABSTRACT
1. INTRODUCTION
1.1. Recall on definitions
In the natural resources industries, there is a clear
distinction between resources and recoverable reserves. The
resources present in a given orebody cannot be extracted
completely for two kinds of reasons :
Technical reasons : equipment size, strength of materials,
organizational problems.
Economical reasons : the cost of extraction is uneconomical.
The term "recovery estimation" adopted in geostatistics refers
to a very particular situation :
the technical environment of the exploitation can be summe-
rized in the definition of a minimum selection volume,
the economical environment can be summerized in the defini-
tion of a selection cut-off applied to the selection volume
average grade.
385
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1,385-420.
o 1984 by D. Reidel Publishing Company.
386 A.MARECHAL
T(zc) = TM *
1 - F (zc)
(2)
I * * *
Q(zc) = TMI'~z(z ) F (dz )
The recover~ tonnage will depend directly on the histogramm
of the estimated block Z*v : whatever the mathematical esti-
mator involved in Z*v, less sample information will result in
a smoother distribution for the estimated grades Z V.
388 A.MARECHAL
2. EARLY SOLUTIONS
2.1. Recovery estimation in the early 70's
As soon as the practice of 3-D block kriging was made possible
on computer around 1970, a block by block estimation of big
open-pit mines was performed : selecting the blocks which kriged
grade was above the cut-off was supposed to provide an estimator
of the recovered tonnage in the open-pit. MARECHAL (1972), DAVID
(1972) showed by comparing these results with the actual recovery
that the method was biased. The smoothing effect of kriging
(Information effect) reduced the variability of estimated grade,
hence changing the proportion of blocks above the cut-off grade.
Estimating the grade of each small block was unrealistics and
moreover without real importance. Which block was to be selected
as ore was a matter of mining practice, when the real estimation
problem was : which proportion of all the blocks to be mined in
the pit will be selected as ore ? This is equivallent to
estimate, locally or globally the histogramm of blocks grades.
Only large panels beeing correctly estimated, is it possible
to derive from the data an estimator of the distribution of
block grades within such panel ? Two kinds of solutions were
proposed, according to whether actual recovery figures were
available or not.
P (zv)
o
= a 0 o + a 0 1 zv + a 0 2 z 2v
For new pa~els, zv is not known, but can be estimated by
kriging z • An unbiased estimator Po* of Po (zv) is then
p. a 1 2 2 2 2
o = a o + E\a ZVI + a. ( ZV ~) +a c (J" IV
Iz(x) =
1 i f z (x) >z (3)
o if z (x) <z
Considering now Z(x), a stationary R.F., Iz (x) defines a new
stationary R.F. with the following characteristics
(8)
(9)
defining ~(z' ) =/
>..~( z. z' ) dz it comes:
t ( X) =~ f i. ( z' ) Iz' ( x. ) dz'
Each term Jp.: ( z' ) l ( x. ) d. is the definition of a function
h (Z )' through the indicator functions. We are back to the
us~al p~esentation (7) of the D.K estimator.
RECOVERY ESTIMATION 393
(13)
5. INDICATOR KRIGING
5.1.• Simplifying the complete indicator cokriging
JOURNEL (1982) proposes to simplify the complete indicator
cokriging shown in (10), replacing it by ordinary kriging with
indicator values at the same cut-off z. In other term, he
396 A.MARECHAL
~ (dz, dz' ) =p( hH)z( dz' )F( dz) +( l-p( h) )F( dz)F( dz' )
In this model, Z(x + h) and Z(x) are equal with the probabi-
lity ~(h), ~r independant with same p.d.f F(z) with the proba-
bility -1- f(lt) .
This is the only model for which the indicator functions are
a set of "autokrigeable" functions, that is functions for which
cokriging reduces theoreticaly to kriging. It can be seen that
such a model is in contradiction with most of the experimental
observations, whereas the "destructurization of high grade" ty-
pical of a bigaussian distribution seems to correspond much
better to the reality. Then, to quote MATHERON (1982), "in the
bigaussian case, substituting indicator kriging to dis- junctive
kriging must result in a very substantial loss of information".
Using indicator functions as service-variable is not recom-
398 A. MARECHAL
Estimate T( Yx ) =R[ .
by TI><= L f'~ ( Y~ )
<I
f.( Y. ) =~-t:. Xn ( Y. )
By linearity of the D.K. estimator of T (Yx), this one can be
expressed as :
TOI( ')'co T
=1Iob n
XDI<n
XOl(
n being the D.K. estimator of the function Xn(Yx)
(22)
Hence
(25)
(26)
a!.=l-),o.p,.(a xl
(27)
}.:p (y.
" n''4
-x• )=Pn (x a. -x) 'v'a.=l.n
Then, whatever be the function T(Yx) known by an expansion
given above, we have :
(28)
(29)
IIlK (
• X
)::l-G(
Y ~ n! ~ " " •
>:
) -~ g( ylH,) y) ()' H (Y ) )(30)
written slightly differently, this formula reads :
400 A.MARECHAL
(
I ~ ( X) = l°Of~ y) dy
(31)
f~( y) = ~_H~( X~~", y) g( y)
I"" ( X) is obtained by integration from a density function
as is the D.K. estimator of any function
(33)
i=O,l. .. N
Z.=z. if Zo~Z.<Z",
~ =2:W·CP( i )Hn( i) N
1-1
u= y- X:x '>j3
y= ¢-~ Z)
(J~
(43)
Yi =)(lxYB +Y~
Then
(46)
(47)
1J
V v E( Z x I Z v =z) =z (48)
(48 bis)
(49)
f(YxIXv)=~Tn(v,v)Xn(Xv)Xn(Yx)g(Yx)
(53)
(54)
(55)
Z =m+o-Y (57)
Z v =¢.(
v Xv )'
Z = ¢ (X ) =~ 'i'n d n ( V) H n( Xv )
v v W n!
no()
The basic hypothese now is that, any pair Yo. • Yj3 are
bigaussian and also any pair \1. • Xv.,
I
with correlation coeffi-
cients fJafj and lai respectively. For such an assumption to be
valid, one has to consider each point sample value Za as located
at random within a small block v., so that the dependance
between any point transformed gr~de Yo. and its "companion block"
transformed grade XVi will be fixed by a unique parameter
r = cov (Yx, Xv).
Applying condition (48) on the above model leads directly to
the D.G. change-of-support formula:
7.6.
Comparison between the practice of different change-of-
support methods
The two major change-of-support formula currently in practice
are the Affine Correction and the D.G. (otherwise termed genera-
lized lognormal permanence). Various checks have been done with
production data (MARECHAL 1975, GUIBAL, REMACRE 1983, DAVID,
1977) but few direct comparisons exist between both methods
(~ruGE 1982), and only one exists comparing theoreticaly exact
results and change-of-support methods (MATHERON 1981). In gene-
ral, for global recovery estimation purposes, for small change-
of-support and not heavy-tailed distributions, the results of
recovered tonnage obtained by the two methods do not differ
considerably, but difference may appear in recovered metal.
The major advantage of Affine Correction is simplicity : the
method can be used without any mathematical modeling of the
point grade histogram. Global recovery estimation can be
obtained directly form the ordered sequence of sample values.
As Zi, -m ~ Zx -m
(j"v (j"
Hence (61)
p( Z v >z) =p( Zx >(j"v
~ ( z-m) +m)
I*(v)=~
So Z N' ~ rank of Zot such that
Z 0. =.Q:..(
(Tv
z-m) +m
The Affine Correction honors the condition of decreasing
benefits for selection on large blocks.
The Affine Correction is exact for gaussian distribution and
good for quasi-gaussian distributions. It preserves lognorma-
lity, but in a rather strange way. Let Zx be any 3 parameters
lognormal.
ZX =Z N +E exp( m+aYx )
RECOVERY ESTIMATION 413
I
The expectation is E( Zv Z, .•• ZN ) = ~ E( Zx Z, f I
The mean of this distribution is the aver~ge in the block of
• ZN ) dx
(68)
(69)
7.9.
Estimation variance of recovery estimators
The situation is very different for global and local recovery
estimation.
The tonnage and quantity of metal in point recovery can be
considered as the average, for the whole ore-cody, of the
service variables Iz(x) and Z I z (x) : after fittir:.g the
x
experimen-
tal variograms of these two variables, tlte estimation variances
can be computed without difficulties using the standard approxi-
mation for global estimation variances. In the case of block
recovery, we are faced with the need for adopting a change of
support formula, the preCision of which cannot be Judged theori-
tically. With an affine correction, Iz(v) and ZvIz(v) can be
lin-
early related to the point variables, so that the estimation
variances will be derived from the variances computed in the
point recovery case. With a D.G. model, Un =E( Zv I Za; can be
computed for each sample point, so that the service-variable
giving directly the tonnage and metal recovered on block can be
studied, their variogram fitted and estimation variances
computed.
Experimental checks have shown that usually a local recovery
estimation is fairly imprecise at the level of each panel,
especially for the quantity of metal at high cut-off values.
The precision is better expressed in term of conditional
variance : this quantity is theoritically accessible only when
using the conditional distribution approach in the D.G. model.
Other estimators (D.K., I.K., etc.) provide only as uncondi-
RECOVERY ESTIMATION 417
8. CONCLUSION
Let us recall two major observations :
Estimating recoverable reserves is of uppermost economical
importance, because it is the only way to account for the two
major factors influencing the selectivity : the support effect
and the information effect. It is only by acknowledging these
effects that geostatistics was able to resolve correctly the
"problem of vanishing tons" (to quote M. DAVID).
In many orebodies, the average grade estimation of each
individual panel is not very precise : what about estimating
the local point grade histogram !
It results from the two above points that performing an
effective recovery estimation deserves a special effort
necessity to use somewhat sophisticated methods based on well
defined probability models. Such methods exist and have been
reviewed in this paper : they are not difficult to understand
if correctly taught to the right people, and the softwares ne-
cessary for applying them exist or can be easily developped.
necessity to devote more human time and thinking to fit and
check the models. It has been noticed that, the more complete
the model, the more satisfying the estimator : but it is known
too, that the more precise the model, the more likely the
model is wrong !
Consider the example of the Petroleum Industry : for such
important problems as seismic data processing or reservoir
production simulations, physically complex models are used,
involving high level mathematics and resulting in large
softwares of high C.P.U. consumption. There are no reasons why
one of the major challenge of the modern Mining Industry,
namely est~mating, efficiently planning and operating big
-selective mine, would not be considered with the same
attention as oil men consider the acquisition and processing
of exploration and production data. Consequently, the future
for recovery estimation methods lies in developping more or
less complete probability models adapted to various particular
cases, together with developping the methods for fitting and
validating these models. Example of such cases are the
negative binomial model adapted to diamonds estimation, the 3
parameters lognormal adapted to the south-african gold
orebodies, or the gaussian transformed model largely used for
estimating disseminated massive orebodies. The estimation me-
thods corresponding to theses models exist : they are the
conditional distribution and disjunctive kriging. It would be
very important, in addition, to develop the simulation methods
corresponding to theses new non-gaussian models.
ACKNOWLEDGMENTS
REFERENCES
Georges MATHERON
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1, 421-433.
© 1984 by D. Reidel Publishing Company.
422 G.MATHERON
T
O(T) = jz(t) dt meT) = 0(T) /T veT) 0(T) - T z(T)
0
In other words, the convex function V(z) and the concave func-
tion OCT) are in duality, in the sense of convex analysis, see
[2J, pp. 104 & Sq. In particular, they satisfy the reciprocal
relations:
V( z) = Sup { Q( T) - z T}
TC[O,1]
(a)
and thus
z
¢(z) =a + bz + 1(z - t) ll(dt)
o
00
v (t)
2
(c)
426 G.MATHERON
(b)
(d)
T[. =
1
u{s.:
J
j EJ, S.nT[.
J 1
fl!}
FS =L p(S.) FS F =
T[
L p(T[.) F
T[i
J j 1
pectively).
Then:
Comments.
h ( X) = E(Z!X )
But it is clear that this comes back to working with the ran-
dom variable H = h (X) = E(Z/X). In other words, the true
grade/tonnage curves depend on the distribution FH of H and not
on the distribution Fn of Z. Even though it might seem paradoxi-
cal at first, we may say that we are mining conditional
expectations ~ rather than grades ~.
5. A TYPICAL EXAMPLE.
than FS ' but more selective than FH (see above), that is:
or
Var(i.nH)
TABLE I
-------------- ------------------------------------------------
CUTOFF GRADE VARIABLE T(z) Q(z) m(z) V(z) x 100
z
-----------------------------------------------
H .872 .949 1.09 51. 3
z = 0.5 ')l;
S
.735
.577
.909
.884
1.24
1.53
54.2
59.5
S
.362
.309
.638
.691
1. 76
2.24
27.6
38.3
S
.177
.183
.413
.538
2.33
2.94
14.8
2(,.4
-------------- -------------------------------------------------
S = f F (z ) [1 - F ( z ) ] dz
a
00
1
S = E(V(Z)) =fV(T) dT
o
!2 S
1S = Cov (z(T), ~ - T)
S ~ 0//3
+00
~ 5 = f
-00
¢ (y ) [2 G(y) - 1] dy
00 (_1)n C
S =.:..1. L; 2n+1
;rr n=O 22n(2n+1) n!
THE SELECTIVITY OF THE DISTRIBUTIONS 433
If F is N (m,02), we find
5 = ojfi
(this value is very near the limit oj 13, so that a normal dis-
tribution is not very far from a uniform one •.• ).
REFERENCES
ABSTRACT
INTRODUCTION
The general problem comes from the fact that the charac-
teristics of many deposits are incompatible with the requirements
of non-linear methods (i.e. strictly stationary hypothesis).
For instance, in the case of porphyry copper, there is a high
grade zone which often leads to preferential sampling at the de-
triment of the border which is poorer. We are thus confronted
with a double problem: presence of a large scale drift and an
irregular preferentially sampled grid. We know that these cir-
cumstances do not have a great effect on kriging with a univer-
435
G. Verly et al. (eds.). Geostatistics for Natural Res(Jurces Characterization, Part 1, 435-448.
© 1984 by D. Reidel Publishing Company.
436 D. GUIBAL AND A. REMACRE
PRELIMINARY CONDITIONS
The data come from one level of a mine where the analysed
copper grades of blast-holes were available. They are considered
to represent the reality. More precisely, the data correspond to
a bench 13 m thick in a porphyry copper type deposit, which is
thus reduced to 2 dimensions. The presence of different zones
can be observed (in our case, slightly higher grades towards th.e
south of the bench), which confirms the type of deposit mentioned
above.
We then have two files: one contains the basic data, the
second one the "reality". The basic data are used to compute the
recovery functions of blocks 10 x 10 x 13 m using the method
mentioned below and the results are then compared to the contents
of the second file. We have obtained 173 basic data, 583 blocks
and 50 panels with 281 blocks.
fn
(1) fey) I --. H (y)
n. n
L
Our objective is to estimatei~l Hn (YVi)/L using three different
sets of assumptions.
a) Multigaussian.
[] (Y )J
n Vi :1G
+"
L
f f 1 fn H (yk /s )
o 1 L L..J n! L n vi Vi
n~2
LOCAL ESTIMATION OF THE RECOVERABLE RESERVES 439
b) Uniform Conditioning.
[H
n
(Y .>] UC
Vl
= E [H (Y
n
.> IY]
Vl
= pn H (ykjs)
vi n
c) Disjunctive Kriging.
A
a,n
pn
as
=.1
L i=1
L: pn
aVi
Then, the DK estimator can be expressed in the form:
f
f
o
- f1 L:a AY +
a a
L:
n~2
~
n!
L: Aa,n H (Y )
n a
We note that for n = 1 all the methods give the same re-
sults, which is very interesting in view of the importance of the
left-hand term of the development of (1).
Zv = I/L L: ¢ (Y ) = L
r vi n
which comes back to the same thing.
HISTOGRAMS
TABLE 1
r-----------------------------------------------------------
t-----------------------------------------------------------
ACTUAL BLOCS Z PANELS BLOCS*
IS. IS.
II. II.
S. S.
I. I. y
I. 1. 2. S. ~. I. 8. 7. I. I. 1. 2. S. ,. I. I. 7. e.
",,'
.... JIl.'~"'.
t ••
...,.
------------------------------------;~-------
",,_.r~' _ _- - - -
VARIOGRAMS
11
2 5 7
10 3 1 8 12
4 6 9
13
2. 2. TH - T IX.HI
T - T IXI
"
1. 1.
B. B.
-1. -I.
-2. Y -2. T
-2. -I. B. 1. 2. -2. -1. B. t. 2.
ECTH/T:J- f'Y lOCH., ITHI /T:J- p3H., ITI
tt2 H~
2. 2.
,
I 1.
1. I
8.
8.
-1-
-1. T -2. T
-2. -1. B. L. 2. -2. -1. B. 1. 2.
TON ~ O~
100. OK 125.
100.
75.
75 .
5iL
51'!.
25.
25.
0. ZC :'; I'!. ZC ;;
TON ~
100.
75.
50.
25.
0. zq
0.7S 1. e I .25 1 .50
i <!~.
7S
IS.
lil
~5.
5.
IL..G "3.d '3 -;"" ~ .'is.Idr"",- oC..S ).d ').~ 4.('~ L.&.C; S.1d ~ c:
TON \
11'10
411.
75.
'lB.
50.
21'1.
25
lB.
II. e. zc(,
1.52.1'12.5 'l.a 'l,S 4.1'14.5 1.5 2.e 2.5 '9.1'1 '9.5 4.1'1 4.5
CONCLUSION
block grades, whereas one would expect this only if the sample
grades were the same.
REFERENCES
Georges MATHERON
ABSTRACT
449
G. Verly et af. (eds.), Geostatistics for Natural Resources Characterization, Part 1,449-467.
© 1984 by D. Reidel Publishing Company.
450 G.MATHERON
W..
1J
W~ W~
1 J L:
n~O
T (v a' vb) /(i) Xnv(')
n n J
(S)
V V
w..
1J
W. H.
1 L:
J n;;:O
T (V a'V b ) / ( i ) XnV(')
n n J
CD) w..
1J
(CS)
(2)
2 - A GUIDELINE.
-A t
P .. (t)
IJ
L e
n
W.
J
n>-O
A. . - c. 8 .. + C. II ..
IJ 1 1J 1 1J
(5) C. 1oJ. II .. c. W. II ..
1 1 IJ J J Jl
ISO FACTORIAL MODELS AND CHANGE OF SUPPORT 453
A Xn
In other words, the condition (5) means that the generator A must
be a self-adjoint operator on L2(iN,F), and from the condition (6),
the spectrum of this operator must be discrete, with ei~en values
- An~ 0 (A o = 0). Moreover, the eigen functions, i.e. the factors
X , must form a Hilbertian basis of L2 (IN,W).
n
(7) H. a. = H. J b.1.+ J
1. 1. 1.+
a.
1.
= p(v+i) b. = i
1.
(v > 0, ° < p < I)
a
-G.(s,t) A G. -[v p + i(p+I)]G. + p(v+i)G. I + i G. I
at 1.
1. 1. 1.+ 1.-
and the initial condition Ci(s,O) = Sl. From the intuitive inter-
pretation of the process, we expect G. to be of the form
1.
Gi(s,t) = H(s,t)[y(s,t)]i
1 dH - v p(I - y) H(s,O) 1
H at =
n (S)
n
= LW,
. J
H (j)sj
n
v
q (I - s )
n
I (I -p s )
n+v
J
\!
H (i) 1- F(-n v+i v'q)
n n '"
p
L
n?O
AV pn H (i) H
n n n
(J') W.
J
<5 ••
lJ
and the stationary distribution of NCO) and N(t) has the bivari-
ate generating function:
(9) G(s,s';t)
456 G.MATHERON
Change of Support.
v n qV(I_s,)n
AV n q (I-s)
(10) L
n~O
n p
(I -ps)
n+v (I-ps') n+v
v
qq' ) q,v(l_s)n qV(I_s,)n
(
q'-(p-p')s-p'qss'
L (I_p's)n+v (I_ps,)n+v
n:?.O
v n v+v" n
q (I-s) q (I-s)
(12)
n+v v+v"
(I-ps) (I-ps)
v v v v Vv
(13) Wi V( Pv ) Wj v( Pv ) "~ Anv( Pv )n HnV(1",
~ Pv ) Hn ("],PV )
n;:.O
Taking into account that mv/my = v/V, the first of the two
relations (14), written with n=l, clearly implies E(Nv/NV)=v/V NV.
In other words, the distribution (13) satisfies Cartier's condi-
tion (2), so that it is an admissible model of change of support.
4 - A DISCRETE/CONTINUOUS MODEL.
r(v)n!
r(v+n)
which make the calculations easy, by taking into account the 1n1-
tial condition L (x) : 1.
o
Change of Support.
(15)
ISO FACTORIAL MODELS AND CHANGE OF SUPPORT 459
Fi rs t Mode 1.
Second Model.
·this implies:
(18)
g(x,x') "
g (x) g (x') LJ r n Av L (x) L (x')
v v n~O n n n
Let us put Q =
tOT
t
i v
X dT and cP (A,\l;X) = E rexp(-AQ - \l X )/X = x].
t L t to
Clearly, the function CPt is determined by the evolution equation:
H
o
9ft = - $2 - $ + A $ = \1
at t t o
( C exp(t/2) )\!
C Ch(Ct/2) + (1+2\1) Sh(Ct/2)
v
~ (A)=( C exp(t/2) \
t [C ch(Ct/4) + Sh(Ct/4)] [Ch(Ct/4) + C Sh(Ct/4)]/
(20)
C = II + 4A
(22)
6(t + t exp(-t) - 2 + 2 exp(-t)) -l
1:2 (exp(-t) - + t)3/2 IV
On the other hand, the model (21) leads to the same value as
the gamma distribution itself, i.e. 2/1V, which coincides with
the value of (22) at t=O. But it turns out that x 3 /(cr t ) 3/2 remains
almost constant if t is not too large:
t x31V/cr~ t x31V/cr~
0 2 1.5 1.9384
0.1 1.99967 2 1.8985
0.5 1.99215 3 1.8056
1 1.97060 4 1.7071
TABLE 1
0.95399 0.99447
x = 0.16 { 0.97950 { 0.99690
0.97893
x = 0.25 { 0.88194
0.87873
{ 0.97698
x = 0.5
0.66911
{ 0.65896 { 0.89919
0.89472
0.68562
x = 1
0.36995
{ 0.36788 { 0.68343
0.32697
x = 2 { 0.11466
0.11460
{ 0.32766
0.13678
x = 3 { 0.03550
0.03573
{ 0.13785
0.00341 0.01994
x = 5 { 0.00347 { 0.02033
x = 7 { 0.00033
0.00034
{ 0.00257
0.00265
ISOF ACTORIAL MODELS AND CHANGE OF S UPPOR T 465
w2 + t 2 /16
exp(t/ll II _n__ ~ ____
f (w)
I - Ft(x) = E (B /b ) exp(-b x)
n n n
(23)
K (x) E (B /b ) (x + I/b ) exp (- b x)
t n n n n
+
G (S,t)
q C exp(qt/2)
(24)
Ch(Ct/4) + C Sh(Ct/4)] [C Ch(Ct/4) + q Sh(Ct /4:J"
C = l(l+p) 2 - 4ps = /q2 + 4p(l-s)
REFERENCES.
Ga ry F. Raymond
Senior Mine Project Engineer, Cominco Ltd.
1. INTRODUCTION
Tabular style lead-zinc ore at Pine Point exists as a number
of small, flat lying, very erratically mineralized orebodies sur-
rounded by barren waste rock. Both exploration and production
ore estimates are a problem. To examine the potential for apply-
ing geostatistics in tabular ore, a detailed back analysis was
made of a 3.2 million tonne (3.5 million ton) mined out orebody
based on exploration drillhole (DH) grades, production blasthole
(BH) grades, and actual mined ore. Ore/waste boundaries estima-
ted using ordinary kriging of BH grades were compared to actual
mined limits based on visual estimation during mining. Explora-
tion conditional probability estimates were calculated using or-
dinary kriging, relative variograms, and a 3 parameter lognormal
approximation to conditional distributions. Comparisons were
469
G. Verly.et al. reds.), Geostatistics for Natural Resources Characterization, Part 1, 469-483.
© 1984 by D. Reidel Publishing Company.
470 G. F. RAYMOND
r
E
0
........ ....
....0
..,
0
:>
...
:>
u u
QI ~
0
> 0
0 .rJ (lJ
<: CD (/)
Q)
X X
m CD +-'
ro
• + .+-'....E
(/)
W
..-
ro
::;
(/)
>
>.
CO
"Cl
Q)
.....C
::;:
V1
+-'
E
.....l
Q)
S-
O
Q)
s...
::;
.....
CJ)
w...
"":'.>
GEOSTATISTICAL APPLICATION IN TABULAR STYLE LEAD-ZINC ORE 473
represent high grade solution channel way ore. The central area
represents lower grade bedding-replacement ore, disseminated ore
and a few local collapses. There is an obvious influence of the
crescent shape of the mining face on these ore/waste boundaries.
Geostatistical study on BH grades involved calculating rela-
tive variograms in various directions, jackknife analysis, and
contouring ore/waste boundaries using ordinary kriging. Relative
variograms were calculated by dividing variogram values at indi-
vidual distance lags by (grade + constant) squared. This approach
will be further discussed under conditional probability. Jackknife
analysis involved removing BH's one-by-one, estimating a grade at
their locations using ordinary kriging of nearby BH's, and comparing
estimates to actual. Since the mining cutoff is determined by
combined Pb+Zn grades, and since initial studies showed similar
variogram shapes for Pb, Zn and Fe, detailed studies as described
in this paper, were based on combined Pb+Zn grades. However,
final methods have since been extended to include individual Pb,
Zn and Fe estimates. Because of grade-density dependency, all
statistics and estimates were density-weighted. Sample density
is calculated from Pb, Zn and Fe grades assuming pure minerals
and constant porosity. Because of the limited vertical extent of
tabular ore, only 2 dimensional (horizontal) estimates were used
in this study.
F~gure 2 shows L37 BH variograms relative to (grade + 0.5%
Pb+Zn) calculated in along-trend and cross-trend directions
(~ 22.5 degrees). The first distance lag represents only a few
sample pairs. Succeeding lags represent 5,000 - 15,000 pairs. A
spherical variogram model with a 2:1 geometric anisotropy has been
fitted to the plot as shown. Relative to grade squared (no additive
constant), the nugget effect of the variogram would be 1.4.
Despite the high variability among BH assays, the BH jackknife
estimate was nearly perfectly conditionally unbiased. In practical
terms this means that, if we ~ to mine to kriged ore/waste
boundaries from BH assays, we would obtain the estimated tons and
grade on average. This conclusion, proven elsewhere by long-term
comparisons to mill head grades (3), is valid only if BH assays
are unbiased and if kriging variance is not too large. Using the
BH variogram and the nearest (corrected for anisotropies) 10 BH
samples per estimate, ore/waste boundaries shown in Figure 3 were
calculated using ordinary kriging. Kriged boundaries are very
similar to those obtained by visual estimates (Figure 1). Compar-
isons between BH kriging and milled ore from the portion of L37
pit studied are shown in Table 1.
Note that milled ore itself is not an exact figure, being
back-calculated through several stockpiles based on mill head
grades and truck counts. Reconciled milled ore tonnage from L37
is thought to be an over-estimate of actual production. From
474 G. F. RAYMOND
Z.O
~
'" 1.0
...,
o
;:
>
>
'"''" Fitted l10del
<i
« 0.5
20 40 60 80 100
DI stance (m)
[
E
0
......2 ......2
OJ OJ
U U
Q)
> ~
0
n Qj
ex:
x
""
x
a:> a:>
• + VI
>,
rei
VI
VI
<::(
:r:
co
0'>
C
-,...-
0'>
-,...-
S-
::.<:
iB
VI
-I-'
',...-
E
',...-
-l
ClJ
S-
o
(V')
ClJ
S-
::l
0'>
-,...-
l..L..
476 G. F. RAYMOND
of ore, and waste plus overburden within the pit actually mined,
and within the pit as redesigned based on conditional probabilities.
20 ~
10 ~
~
.&J
"-
..
: 5 ~
o
."
~
U
o BH Assays
l~ o A hploration OH's
30 50 70 90 98 99_5
'.
'I~n
..
o
.,g
«
:::
~
'"
:c
o 0
• +
<-
ttl
C
o
'...-
~
'...-
"o
C
U
C
o
.
A
. .
482 G. F. RAYMOND
REFERENCES
1. David, M.: 1977, "Geostatistical Ore Reserve Estimation".
Elsevier, 364 pp.
2. David, M.: 1969, The Notion Of Extension Variance And Its
Application To The Grade Estimation Of Stratiform Deposits.
In: "A Decade Of Digital Computing In The Mineral Industry."
AIME, pp. 63-81.
3. Raymond, G.: 1979, Ore Estimation Problems In An Erratically
Mineralized Orebody. CIM Bulletin, Vol. 72, No. 806, pp. 90-98.
4. Raymond, G.: 1982, Geostatistical Production Grade Estimation
Using Kriging In Mount Isa's Copper Orebodies.
Proc. Aust. IMM, No. 284, pp. 17-39.
5. Skall, H.: 1975, The Paleoenvironment Of Pine Point Lead-Zinc
District. Economic Geology, Vol. 70, pp. 22-47.
LOCAL ESTIMATION OF THE BLOCK RECOVERY FUNCTIONS
Amilcar Soares
INTRODUCTION
G. Verly et al. (eds.) , Geostatistics for Natural Resources Characterization, Part 1,485-493.
© 1984 by D. Reidel Publishing Company.
486 A. SOARES
BLOCK-RECOVERY FUNCTIONS
Defining:
1 iff X. ~ X
1 0
(X. )
1
o iff X. < X
1 0
z ¢ (x.), with:
v. v 1
1
MULTIGAUSSIAN APPROACH
N
Hence E { Xi/ya} L Aa(i) YN =
a=l "
The weights Aa are given by a simple kriging system.
1
where S J r x,i dx
v
v
Y. 1 J P.
E Yv ' } dx
J 1 V·1 v·1 J,X
E Y. Y is denoted by P
1 X i'X
S =
1
V.
I (j, ,/ Px,y dx dy) (4)
1 1 1
~ "~real"
~amples r small block - v 5x5 m
l
1O;>------~+-----9---~
:J
0
0
PANEL 1 PANEL 2
'"
I -SO m -
'1
Structural analysis
Results
For each class the relative deviation for the 3 methods are
shown in Fig. 2, 3 and 4. These relative deviation are then
average over the 10 panels considered, c.f. Table I.
From Fig. 2 to 4 and Table I it can be withdrawn the follow-
ing non trivial conclusions:
·eo -8
-, -'00:
I
I
l I
, 4 •j ) 4 Ij 6 7 e 9 HI F'AN£lS
2 l S 6 1 PANELS
• ••
Fig. 2 : Relative deviation values of the Fig. 3: R~Jative deviation values of the
1 st class along the 10 panels 2 class along the 10 panels ....
'oQ
492 A. SOARES
10 10 10
12: IFi-Fi*1 (t) 1.L. Fi -Fi * ('!;) 1.L. (Fi-Fi *) 2
N i=l F1. N i=l Fi N i=l
Z {37.5% 42.9 37.9 37.8 -21.2 -10.3 - 8.6 48.4 53.1 55.1
37 .5,!; (z (45'1; 12.4 13. 15.4 - .8 - 9.4 -12.4 69.6 59.5 92.4
m<z 39.5 38.7 39.1 -20.6 -14.8 - 4.5 51. 50.1 59.0
Average of
31.6 29.9 30.1 -14.2 -11.5 - 8.5 56.3 54.2 68.8
the 3 classes
CONCLUSIONS
ACKNOWLEDGEMENTS
REFERENCES
Georges Verly
ABSTRACT
I - INTRODUCTION
495
G. Verly et al. (eds.), Geostatistics for Natural Resources Characterization, Part 1,495-515.
o 1984 by D. Reidel Publishing Company.
496 G.VERLY
= { o if not
(1)
To L
t(zc) = -- L
iclzyCl») (2)
L 1=1
To L
q(zc) = -- L zyCl) -i clzy( 1)) C3)
L 1=1
To L
r(zc) = L hlzyCl») (4)
L 1=1
where hlzy(l») is some function of the smu grades. From now on.
this general expression will be used.
C5)
To L
EnR(zc) = -- ~ En{hlzy(l)]}
L 1=1
To L +co (6)
= -- ~ f h(z)ef y( ll(zI(H»edz
L 1 =1 -co
= En{(RCzc)]Z} - [EnR(zc)]Z
T02 L L +co +co (7)
= - ~ ~ f f h(z)eh(z')·f y( lly(1')(z.z'I(N»edzdz'
J2 1=1 1'=1 -co -co
-[E nR(zc)]2
1 M 1 M
Zy(l) ~ - ~ Z(i) =- ~ ~(Y(i») (10)
Mi=1 Mi=1
( 13)
where
Note that the conditional mean vector Vila (12) is none other
than the SK estimate of the Y(i)'s.
One possible model for the smu grades is the Discrete Gaus-
sian model [11].
where [Y(1), •.. ,YCM»)' is a M-normal vector with a zero mean vec-
tor and a covariance matrix fully characterized by the variogram
of YCX). The cdf (19) is then obtained by using a Monte Carlo
method.
502 G.VERLY
The conditional cdf of Zy(l) used in eqs. (5) and (9), can
be expressed as
N
L Ai-r(i,a) = R(a,1) (l =1 to N
i =1
where r(i,a) = E{Y(i).Y(a)} and ( 21>
R(a,1) = E{Y(cr.)-BCI)}
(22)
"OESPIKE" PR0CEOURE
RAND0MLY
ACCI!IRDING T0
LilCAL AVERAGES BJ ESTIMA TED PDF'S
~y(hl
Al VARl0GRAHS
h v
k =1. NB (23)
A POINT MULTIVARIATE NORMAL DISTRIBUTION 505
+w +~
Freq.7.
SAMPLES t
SPIKE 31.6%
AVERAGE .59
STATISTICS VARIANCE .384
10
2 3 Z
IV - CASE STUDY
~y(h)
):(
):(~------------~~):(~----
°O~------~------L-------~------~--~h~
100 200
cutoff percentil e
I value
-
1 0 0
2 0.001 31. 5
3 0.170 41.1
4 0.360 51.2
5 0.579 60.1
6 0.989 78.0
7 1. 423 88.7
2.2
2.0
1. 8
1.6
w
o
~ 1.4
C,!)
~ 1.2
<
ffi> 1.0 ......
<
w 0.8
::>
Ct:
.... 0.5
......
0. 4
C0R. C0EF .• 0.94
0.2 REl. BIllS · -3.5%
1.0 2.0
O.g 1.8
0.8 ...J
1.6
<
I-
0.7 W 1.4
L
L1J U.
C9 0.5 (SJ 1.2
z< >-
z
(SJ
r-
~ 1.0 c-
0.5 I-
.
I- 0 •• / •
Z
L1J <
:::J 0.4 :::J 0.8
0:: o
I-
W
0.3 :::J 0.0
r:r
I-
0.2 0.4
ceR. C0EF. 0 0.9.
C~R .C0EF . 0 0.86
0
REL. BI AS 2.42
0
0.1 l·/ REL. BIAS '5 . 5~ 0.2
;~~
0,0 I I , ,
, , , I , , , I I , , ,, , I I , , ' I ,I,w ,
::.. Yf ~ , I
O.Ow; .. ".;' ." ,I " , "'-1, .,., • • •I l
a 0 .2 0 .4 0.6 0 ,8 1 1.2 1.4 \.6 1.8 2 2 ,2
o 0.1 0.2 0 .3 0.4 0.5 0.5 0.7 0.8 O.g ~IG QUANTITY 0F MET AL (/To) ; "3 (0 ,170;0
MG T0 t\NAGE liTo); CUT0FF "3 (0.170:<)
PANEL (20X20l / SMU (5X5)
PANEL (20X20l / SMU (5X5)
o
<:
ITl
Fi g. S. Recovered tonnages and quantities of metal, true versus estimated, given the 3rd cutoff ::0
t'"'
(0.170), This cutoff corresponds to the 41th percentile of the distribution. <:
;!>
:3
z
o-j
s::
c:
2.2 ti
1.0 =2
;!>
2 .0
0.9 ;;:
'"
o-j
1. 8 t'1
0.8 Z
o
...J 1.6
< s::
0.7 I- ;!>
'r"
w 1.4
w ~ 1:1
~ 0 .15 LL Vi
o-j
is)
z 1.2
z >-
~ 0.5 I- '0;"
..... c:
l- 1.0
. ::l
w Z f ,• o
~ 0.4 < /.' z
I-
::> O.B
CI
0.3 w
::> 0.5
0::
I-
0.2
0.4
C0R. C0EF. - 0.92
CHR. C0EF. - 0.93
0.1 REL. BIAS - ·3.q~
0.2 l -' , REL . BIAS - '1.3X
./ ' . -
0.0 bl'ir:;1*w.1~ •• I • I III •• I. I I II I • L.u I I ! ! I I I IZLLI.II ! I I! I. II !! I II '
0 .0 1•.• 1 •.• 1 ••• 1. .1 .n
o 0 .1 0.2 0 .3 0 .4 0.5 0.5 0.7 0.8 O.P o 0 .2 0 .4 0.5 0.8 1.2 1.4 1.6 1.81 2 2.2
MG T0NNAGE liTo); CUHlFF ~5 (0.579 4) MG QUANTITY 0F METAL (ITo) ; a5 (0 .5794)
PANEL (20X20) I S'1U (5X5) PANEL l20X20l I SMU (5X5)
fi g. 6. Recovered tonnages and quantities of metal. true versus estimated. given the 5th cutoff <J>
(0.579). This cutoff corresponds to the 60th percentile of the distribution.
N
-
• .. 1... 1 ... 1 ' , , 1 . , ,...1 '1 ' , . 1 ' I·', ,, ',I ' , , 2.2
1.0
2.0
0.9
1.8
0.8
-'
1.6
0.7 <
I-
UJ w 1.4
L
~ O.B
LL
z IS) 1.2
z
~ 0.5 >-
I-
~ 1.0
UJ I-
~ 0.4 z
I-
< 0.8
:::>
a
0.3 w
:::> O.B
0:
0.21- / I-
C0R. C0EF. - 0.89 0.4
0 .1 I- / REL. 81AS - O . 9~ CilR . CilEf . - 0 .90
REL. BIAS - 6 . 3~
J 0.2l . /
0.0 V!,T'. J ,',', i i ,. 1 , " " , ,I "I", ",I, 1 ' 0
O.OlZl'"j"'Li!7!!",-.. I II , I " , I . . I" I, I, ,n
o 0.\ 0.2 0.3 0.4 0.5 O.B 0.7 0.8 0.9
MG T0NNAGE liTo); CUHlFF 11.7 (l .423;() o 0 .2 0. 4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 2 .2
MG QUANTITY 0F METAL (/Tol: "7 (1.423, )
PANEL (20X2OJ I SMU (5X5)
PANEL (20X20l I SMU (5X5)
o
<:
tTl
fig. 7 . Reoovered tonnages and quantities of metal, true versus estimated. given the 7th cutoff ::<'
l""
(1. 423), This cutoff oorresponds to the 89th percentile of the distribution. -<
A POINT MULTIVARIATE NORMAL DISTRIBUTION 513
o
t-
--ACTUAL
/'
-' - - ESTIMATED
<: /'
.....
w /'
1: LU
/' D
LL
C>l
<
0::
>- /. /' GRADE (!)
t-
.....
Z
<:
::l
CI
0.5 1.0
~
D
Z
-<
W
(!)
<
z ~
z
C>l
t-
the relative bias never exceed 6Y. in absolute value, and most of
the time is less than 3.5%. The higher relative bias on the
recovered tonnages and grades at low cutoffs (respectively -13.8Y.
and 12Y. given the .001 cutoff) is due partly to some panels for
which the available information does not reflect the content.
However, the grade of these panels is generally low and the
amplitude of this bias quickly decreases as the cutoff increases.
v- CONCLUSIONS
REFERENCES
15. Switzer, P., and Parker, H., (1975), "The problem of ore ver-
sus waste discremination for individual blocks: the log-
normal model". NATO Advanced Study Institute, D. Reidel
Publishing Company, Dordrecht (Holland), p. 203-218.
16. Verly, G., (1983). "The multigaussian approach and its appli-
cations to the estimation of local reserves". Journal of
Mathematical Geology, Vol. 15, No.2, pp. 263-290.
THREE-DIMENSIONAL, FRECUENCY-DOMAIN
SIMULATIONS OF GEOLOGICAL VARIABLES
Leon Borgman
University of Wyoming, Laramie, Wyoming
M. Taheri
INTEVEP, Caracas, Venezuela
R. Hagan
Anaconda Minerals, Denver, Colorado
517
G. Verly et af. (eds.), Geostatistics for Natural Resources Characterization, Part 1,517-541.
© 1984 by D. Reidel Publishing Company.
518 L. BORGMAN ET AL.
Two-Dimensions:
= L: L: A(fl,f2,f3)exp[i211(flxl+f2x2)]dfldf2
( 1.4b)
Analogous formulas may be written for three or more
dimensions.
-1
M = (Nt,x ) (1 .5)
N-l
t,x L W
n
exp[-i2'TTmn/N] (1.7a)
n=O
N-l
W = t,f L A exp[i2rrmn/N] (1.7b)
n m
m=O
This is called the discrete or finite Fourier
transform pair and is exact for any finite sequences
of complex numbers [4].
N-l
Sm = ~x L Ckexp[-i2nmk/N] (l .8)
k=O
where
(1. 9)
N,-' Nr '
i:':,X,i:':,X Z 1: 1:\J eXPL -i 2rri!!!1~1 + !!!2.!!.2.)] (1.10)
n, ,nz N, N2
nl=O n2=O
N1-l N2-1
i:':,f 1i:':,f 2 1: 1: A exp[i21T(.~J~1+ .!!12!12)] (1.11 )
ml,m2 Nl N2·
ml=O m2=O
In a similar- way. the discr-ete Four-ier transform can be
generalized to three or more dimensions.
CB = BL (2.1)
!
W= Il + BL2 Z (2.2)
(2.3)
(3.1 )
(3.2)
(3.4)
COy
m1,Nz-m z
S Ol
COy (3.5)
o Sm"N 2-m2
N,-'l: NZ-l
llX Z l:
kl=O kZ=O
A -m N -m = complex conjugate of A
n1 I l' 2 2 ml ' m2
(3.7)
An -m ,m = CO[i'P lex conj ugate of Am ,N -m
l l 2 l 2 2
t:
analogies to the corresponding integral relations
Let
(4.3)
(4.4)
(4.6)
and let
p = h/ho (4.8)
g(p) (4.9)
(4.12 )
Ir I~ IP
21Tab
S*(r) = -
(1 - ~+~)
2 . 2
JO(rp) m dp (4.14)
(4.16 )
(4.18)
(4.21 )
where
(4.22)
~~ijREE-DIMENSIONAL SPECrRA
Cov
[ 6c -q~ -~ ~l
S 0
(6.1)
q cos
(7.2)
(8.1)
(8.2)
(8.3)
(8.4)
(8.5)
(E) ....
lJ exp( i 21T!!li~) (8.6)
~x EW (8.7)
- [!l )
-W-- !!2
- (8.8)
,...,
where ~2 is the conditional simulation of .\!2. given .\!l
= ~l' Let.A be the Fourier coefficients for an
unconditional simulation, .W. Then
~
A (8.9)
3-D FREQUENCY-DOMAIN SIMULATIONS OF GEOLOGICAL VARIABLES 535
o
o A B
NI -I
Figure 3_1 Geometric arrangement of 2-0
Frequency domain structure
u.
...,
0-,
TARLE 7.1 Real part of Fourier coefficients for firs,t horizon of Example R.
Ml= 0 2 4 6 8 9 10 11 12 13 14 15
M2
0 -19.6 -12.4 -59.1i 27.6 1.2 6.7 -.7 -1.6 -.1 1.6 .7 6.7 1.2 27 .6 -59.6 -12.4
1 -22.3 -10.5 -49.0 -20.9 5.9 .7 -3.5 -.5 0.0 .2 2.4 .7 2.4 1.3 -68.1 -12./j
9.4 7.6 -16.1 I,. 1 3.7 1.0 -.1 0.0 0.0 -.1 .7 .2 -.9 7.0 8.6 -3.3
-.9 - .1 0.0 -.2 .2 0.0 0.0 0.0 0.0 0.0 0.0 .2 -.1 1.1 •1 .3
4 0.0 -.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 C.O 0.0 0.0 0.0 0.0 0.0 0.0 0.0 -.1
13 .9 .3 -.1 1.1 -.1 -.2 0.0 0.0 0.0 0.0 0.0 0.0 .2 -.2 0.0 -.1
14 9.4 -3.3 8.6 7.0 -.9 .2 .7 -.1 0.0 0.0 -.1 1.0 3.7 4.1 -16.1 -7.6
15 -22.3 -32.6 63.1 1.3 2.6 .7 2.4 .2 0.0 -.5 -3.5 .7 5.9 -20.9 -49.0 -10.';
TABLE 7.2 Imaginary part of Fourier coefficients for first horizon of Example fl.
Ml= 0 2 4 6 8 9 10 11 12 13 14 15
M2
0 0.0 -39.3 25.5 -1.9 -21.0 -7.9 -14.2 -.1 0.0 .1 14.2 7.9 21.0 1.<) -25.0; 39.3
.2 -52.2 -4R.9 7.9 -R.3 -13.5 2.3 -.6 .2 -.3 -2.8 6.0 19.5 -';7.3 26.2 -68.9
2 -4.7 1.4 7.7 -9.3 3.3 -1.2 -.3 -.1 -.1 .2 -2.1 -[.4 -.9 -12.2 -4.1 2.4
3 .4 -1.1 .2 -.7 .4 0.0 .1 0.0 0.0 0.0 -.1 •1 0.0 0.0 0.0 -.6
4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 t""
7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ",
8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 1).1) 0.0 0.0 0.0 0.0 0.0 0.0 0
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ::0
9 C)
10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ;;::
0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 ;.-
11 0.0 z
12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
13 -.4 .1) 0.0 0.0 0.0 -. I .1 0.0 0.0 0.0 -.1 0.0 -.4 .7 -.2 1.1 ...,m
4.7 -2.4 4.1 12.2 .9 1.4 2.1 -.2 .1 .1 .3 1.2 -3.3 9.3 -7.7 -1.4 ;.-
14 t""
2.3 .3 -.2 .f, ~.3 ll.5 8.3 -7.9 48.9 52.2
15 -.2 611.9 -21\.2 57.3 -IQ.5 -6.0
w
TABLE 7.3 Simulated values for first horizon of Example B 6
'Tl
m
I:)
Ml= 0 2 4 5 6 R 9 10 H 12 13 14 15
'c:::"
m
M2 Z
-.8 -.2 .8 .CJ 0.0 -.6 -.4 -.7 -1.2 ("l
0 -.8 -.5 -1.0 0.0 1.2 .7 -.4
.3 -.7 -1.4 -.8 .f> 1.2 .5 -.5 -.8 -1.0 -1.2 ><:
1 -.1 .4 .3 .3 .9 6
2 .5 1.2 .4 .5 .5 -.1 -1.0 -1.CJ -1.4 .3 1.2 .7 -.4 -.9 -1.2 -1.0 0
3 .9 1.7 .8 .5 .2 -.5 -1.0 -2.1 -1.8 -.2 .R .7 -.3 -.CJ -1.0 -.8 =::
4 1.0 1.9 1.0 .5 0.0 -.5 -.7 -1.9 -2.0 -.7 .2 .5 0.0 -.5 -.7 -.5 >
5 1.0 1.7 .9 .4 -.1 -.2 -.1 -1.4 -2.0 -1.3 -.5 .3 .3 -.1 -.2 -.1 Z
6 .8 1.2 .5 .3 .1 .4 .7 -.7 -1.9 -1.7 -1.0 .2 .6 .2 .1 .1 '"
0.0 -1.6 -1.9 -1.3 .2 .9 .5 .3 .2
:i
7 .5 .6 .1 .2 .3 1.1 1.5 c:::
1.6 2.2 .5 -1.2 -1.8 -1.3 .3 .9 .5 .2 .1 t-
8 .1 -.1 -.3 .2 .6
9 -.4 -.7 -.7 .1 .9 2.0 2.5 .9 -.7 -1.5 -1.2 .2 .8 .5 .1 -.2 ~
10 -.9 -1.1 -1.1 0.0 1.0 2.1 2.5 1.1 -.3 -1.1 -1.1 .1 .6 .5 .1 -.6 0
.2 -.7 -1.1 -.5 .2 .5 .1 -.9 Z
11 -1.5 -1.5 -1.3 -.2 1.1 2.0 2.1 1. I
12 -1.9 -1.8 -1.6 -.3 1.2 1.7 1.6 .9 .6 -.2 -1.0 -1.0 -.1 .5 .2 -1.2 '"0
'Tl
13 -2.1 -1.9 -1.3 -.4 1.3 1.5 1.0 .7 .8 .3 -.8 -1.2 -.4 .5 .2 -1.3
-.6 .4 .1 -1.3 P
14 - I. 9 -1.8 -1.8 -.4 1.4 1.3 .5 .3 .8 .6 -.3-1.1 m
15 -1.4 -1.3 -1.5 -.2 1.4 1.1 .1 -.2 .4 .8 .3 -.6 -.6 .1 -.2 -1.1 0
t-
O
P
TABLE 7.4 Real part of Fourier coefficients for second horizon of Example 8.
n>
t-
Ml= 0 2 3 4 5 6 8 9 10 11 12 13 14 15 <:
M2
>
o -20.2 24.0 -63.7 26.3 -2.9 5.1 -.4 -.7 .4 -.7 -.4 5.1 -2.CJ 26.3 -63.7 24.0 ';;"
2.4 -3.3 -.2 -.2 0.0 1.3 1.3 -2.7 -2.7 6CJ.6 -17.1 t=
1 -49.5 -15.2 -54.2 -25.3 1.3 t-
2 5.3 7.1 -14.0 4.5 1.7 .9 0.0 0.0 0.0 .1 .6 .3 -1.0 4.5 6.5 .1 m
3 .8 .3 .3 -.2 .1 0.0 -.1 0.0 0.0 0.0 0.0 -.2 -.1 1.1 .3 .6 '"
4 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
5 0.0 0.0 0.0 0.0 0.0 0.0 . 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 D.!]
9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
13 .8 .6 .3 1.1 -.1 -.2 0.0 0.0 0.0 0.0 -.1 0.0 .1 -.2 .3 .3
14 5.3 .1 6.5 4.5 -1.0 .1 .6 -.1 0.0 0.0 0.0 .9 1.7 w
'"
-2.7 1.3 1.3 0.0 -., -.2 -3.1 2.1-+ 1.3 -A:j :~t9 l~:; ->
15-49.5 -17.1 69.1i -2.7
<A
W
TABLE 7. ') Imaginary part of Fourier coefficients for second horizon of Example B .. 00
Ml= o 4 6 8 9 10 II 12 11 II. 15
M2
o 0.0 -58.9 2.3 -11.7 -22.3 -10.5 -13.2 -.4 0.0 .4 13.2 10.5 22.3 13.7 -2.3 53.9
1 1.9 -40.9 -57.5 9.8 -4.7 -13.~ 4.2 -.3 .1 -.3 -2.3 3.5 22.4 -53.5 33.2 -71.6
2 -8.9 5.4 7.0 -8.9 5.4 -1.5 -.1 -.2 -.1 .1 -2.3 -1.5 0.0 -11.1, -3.1 2.0
3 .5 -1.4 .2 -.8 .3 0.0 0.0 0.0 0.0 0.0 -.1 .1 .1 0.0 - .1 -.8
4 0.0 0.0 -.1 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
5 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
6 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
7 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
8 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
9 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
10 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
11 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
12 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0 .1 0.0
13 -.5 .8 .1 0.0 -.1 -.1 .1 0.0 0.0 0.0 0.0 0.0 -.3 .3 -.2 1.4
1', 8. I) -2.0 8.<) 11.4 0.0 1.5 2.3 -.1 .1 .2 .1 1.5 -5.4 8.9 -7.0 -5.1,
15 -1.9 71.6 -38.2 53.5 -2~.4 -3.5 2.3 .3 -.1 .3 -4.2 13. Ii 4.7 -9.8 57.5 40.9
Ml= o 2 4 8 9 10 lJ 12 13 14 15
M2
o -.9 -.1 -.3 .2 .3 .1 -1.1 -1.5 -.9 .3 .7 -.2 -.3 -. ') -.8 -1.5
1 -.2 .8 .5 .6 .6 -.2 -1.3 -2.6 -1.5 .2 1.1 .2 -.9 -1.0 -1.2 -1.5
2 J~ 1.6 1.2 .9 .4 -.5 -1.3 -2.1 -1.9 0.0 1.1 .4 -.8 -1.2 -1.4 -1.3
3 .3 2.2 1. Ii 1.1 .3 -.6 -1.2 -2.1 -2.2 -.4 .7 .5 -.n -1.2 -1.1+ -1.1
4 1.1 2.1, 1.9 1.1 .3 -.4 -.8 -2.0 -2.3 -.9 .1 .3 -.2 -.8 -1.0 -.7
5 1.2 2.1 1.6 1.0 .3 0.0 -.1 -1.4 -2.2 -1.5 -.6 .. 2 .2 -.4 -.5 -.3
6 1. 1 1.8 1.2 • I) .5 .6 .Il -.7 -2.0 -2.0 -1.2 .2 .6 .1 -1 .1
7 .9 1.2 .7 .7 .1 1.4 1.7 0.0 -1. 7 -2.2 -1.5 .2 .9 .4 .3 .3
8 .6 . ') .1 . ') 1.0 2.0 2.4 .7 -1.2 -2.0 -1.5 .1 1. 0 .6 .4 .2 r
9 .1 0.0 -.3 .1, 1.2 2.1, 2.7 1.1 -.7 -1.7 -1.4 .2 • Q .. (1 .4 0.0 to
o
10 -.4 -. ') -.6 .2 1.3 2.4 2.6 1.2 -.2 -1.3 -1.3 -. I .6 .6 .3 -.1
11 -.9 -.9 -.9 0.0 1.3 2. I 2.0 1.0 .1 -.8 -1.2 -.4 .3 .6 .2 -.7 ""
12 -1. 5 -1.2 -1.2 -.1 1.2 1.6 1.1 .6 .3 -. ') -1.2 -.8 -.1 .6 .2 -1.1
~
;..-
13 -1.8 -1.4 -1.3 -.2 1.1 1.1 .4 .2 .3 -.1 -.9 -1.1 -.3 . 'j .1 -1.4 Z
14 -1.'8 -1.3 -1.3 -.2 1.1 .7 -.2 _.1+ .1 .1 -1.0 -.5 ..1 -1.5 tTl
-.5 -. I >-l
15 -1.5 -.9 -.9 -.1 1.0 .4 -.7 -.9 -.3 .3 .1 -.7 -.7 0.0 -.4 -1.5 ;..-
r
TABLE 7.7 Data covariance computed from the simulation shown in Table 8.7. mean = -.076, variance .985 6'T1
::<:I
tT1
,0
~1l= 0 2 3 4 5 6 7 8 9 10 11 12 13 14 15 c::
tT1
1'12
-.b,
Z
0 1.0 .7 .1 -.3 -.4 -.3 -.1 0.0 0.0 0.0 -.1 -.3 -.3 .1 .7 (")
1 .9 .7 .2 -.3 -.4 -.2 -.1 0.0 0.0 0.0 -.2 -.1 -.4 -.3 .1 .F; -<
2 .7 .6 .2 -.2 -.3 -.2 0.0 .1 0.0 -.1 -.2 -.3 -.3 -.3 0.0 .4 6
0
3 .5 .5 .2 -.1 -.2 - .1 0,( .1 .1 -.1 -.2 -.2 -.2 -.2 - .1 .2 ;;::
t. .2 .3 .2 0.0 -.1 0.0 .1 .1 .1 -.1 -.1 -.1 -.1 - .1 - .1 0.0 ;.-
5 -.1 .1 .1 .1 0.0 0.0 .1 .1 .1 0.0 -.1 -.1 0.0 0.0 -.2 -.2 Z
Vl
6 -.3 -.1 0.0 .1 .1 .1 .1 .1 .1 0.0 0.0 0.0 .1 0.0 -.2 -.3
7 .-4 -.3 0.0 .1 .1 .1 .1 .1 .1 .1 0.0 .1 .1 .1 - .1 -.4 a::
c::
8 -.5 -.4 -.1 .1 .2 .1 .1 .1 .1 .1 .1 .1 .2 .1 - .1 -.4 t""
;.-
9 -.4 -.4 -.1 .1 .1 •1 0.0 .1 .1 .1 .1 .1 •1 .1 0.0 -.3 o-j
10 -.3 -.3 -.2 0.0 .1 0.0 0.0 0.0 .1 .1 .1 .1 .1 .1 0.0 -.1 0Z
11 -.1 -.2 -.2 0.0 0.0 -.1 -.1 0.0 .1 .1 .1 0.0 0.0 .1 .1 .1 Vl
12 .2 0.0 -.1 -.1 -.1 - .1 -.1 -.1 .1 .1 .1 0.0 -.1 0.0 .2 .3 0
13 .5 .2 -.1 .-2 -.2 -.2 -.2 -.1 .1 .1 0.0 -.1 -.2 -.1 .2 .5 'T1
14 .7 .4 0.0 -.3 -.3 -.3 -.2 -.1 0.0 .1 0.0 -.2 -.3 -.2 .2 .6 C'l
tT1
15 .9 .6 .1 -.3 -.4 -.3 -.2 0.0 0.0 0.0 -.1 -.2 -.4 -.3 .2 .7 0
t""
0
C'l
TABLE 7.8 Data covariance computed from the simulation shown in Tables 8.10 mean = -.079, variance = 1.146 n
;.-
t""
1'11= 0 2 3 4 6 8 9 10 11 12 13 14 15 <:
~12
;.-
::<:I
0 1.1 .8 .2 -.3 -.4 -.3 -.1 .1 .1 .1 -.1 -.3 -.4 -.3 .2 .8 :;
1 1.1 .8 .2 -.2 -.4 -.2 0.0 .1 .1 0.0 -.1 -.1 -.4 -.3 .1 .7 t=
t""
2 .9 .7 .2 -.2 -.3 -.2 0.0 .1 .1 0.0 -.2 -.3 -.3 -.2 .I .6 tT1
Vl
3 .6 .5 .2 -.1 -.2 -.1 0.0 .J 0.0 -.1 -.2 -.3 -.2 -.2 0.0 .3
4 .3 .3 .2 0.0 -.1 0.0 0.0 0.0 0.0 -.1 -.2 -.2 - .1 -.1 -.1 .1
5 0.0 .1 .1 .1 0.0 0.0 0.0 0.0 -.1 -.1 -.2 -.1 0.0 0.0 - .1 -.1
6 -.3 -.1 .1 .1 .1 -.1 0.0 0.0 -.1 -.2 -.1 0.0 .1 0.0 -. I -.3
7 -.4 -.3 0.0 .2 .2 .1 0.0 -.1 -.1 -.1 -.1 0.0 .1 .1 -.1 -.4
8 -.5 -.4 -.1 .1 .2 .1 -.1 -.1 -.1 -.1 -.1 -I .2 .1 -.1 -..II-
9 -.4 -.4 -.1 .1 .1 0.0 -.1 -.1 -.1 -.1 0.0 -1 .2 .2 0.0 -.3
10 -.3 -.3 -.1 0.0 .1 0.0 -. J -.2 -.1 0.0 0.0 -1 .1 .1 .1 -.1
11 0.0 -.1 -.1 0.0 0.0 -.1 -.2 -.1 -.1 0.0 0.0 0.0 0.0 .1 .1 .1
12 .3 .1 -.1 -.1 -.1 -.2 .2 -.1 0.0 0.0 0.0 0.0 -.1 0.0 .2 .3
13 .6 .3 0.0 -.2 -.2 -.1 -.2 -.1 0.0 .1 0.0 -.1 -.2 -.1 .2 .5
14 .9 .6 .1 -.2 -.3 -.3 -.2 0.0 .1 .1 0.0 -.2 -.3 -.2 .2 .7
w
'"
15 1. I .7 .1 -.1 -.4 -.1 -.1 0.0 .1 •1 0.0 -.2 -.'.. -.2 .2 .R
'"
540 L. BORGMAN ET AL.
Alain GALLI
ABSTRACT
1. INTRODUCTION
G. Verly et al. (eds.) , Geostatistics for Natural Resources Characterization, Paty 1, 543-557.
© 1984 by D. Reidel Publishing Company.
544 A. GALLI ET AL.
2. IMPLEMENT AnON
z(x) is the field measured at points Z (Xj), the Vi (x) are mutually
orthogonal random functions representative of differents scales of
variability.
.It is obvious that relation (1) does not specify the expectation of
the VI. It can be easily seen, for instance by writing the non bias
condition, that if :
In case 1), the expression (3) is still .yalid, but the cokriging is
performed without conditions on the ~ IJ.
It can then. easily be shown that, if Z * (x) i.s the kri~ing estimate
of Z (x), VI * (x) the cokriging estimate of VI (x), m the kriging
estimate of m, we have:
On the*other hand, we know from the dual formalism (cf. [2J, [3J)
that Z (x) can be written as :
Unique neigbourhood
Moving neighbourhood
3. MAGNETIC DATA
Figure 1
Magnetic values
(The rectangle shows the test zone for stability study)
Figure 2
There are two possibilities for the term m(x) : to work either on
anamorphosed data with a zero mean, or on raw data with the
arithmetic mean removed. In this case we have no universality
condition and m(x) = 0, or else we can work with a non-zero mean. In
this case m(x) is the constant term of the drift, similar to a local mean,
and is estimated by kriging on the same neighbourhood as the 2
components.
Raw data
3a 3b
VI * estimated with m * estimated with
8 neighbourhood points 8 neighbourhood points
550 A. GALLI ET AL.
Raw data
4a 4b
y 1 * estimated with m * estimated with
90 neighbourhood points 90 neighbourhood points
Anamorphosed
data
5a 5b
y1* estimated with m * estimated with
16 neighbourhood points 16 neighbourhood points
FACTORIAL KRIGING ANALYSIS 551
Anamorphosed
data
6a 6b
V1* estimated with m* estimated with
90 neighbourhood points 90 neighbourhood points
Conclusions on Stability
Anamorphosed
data
7a 7b
V1 * estimated with V2* estimated with
16 neighbourhood points 16 neighbourhood points
552 A. (jALLI P'I' AI..
Anamorphosed
data
8a 8b
yl * estimated with estimated with
90 neighbourhood points 90 neighbourhood points
gure 9b
Raw data: y2* estimated with 70 neigbourhood points
554 A. GALLI ET AL.
Figure 9c
Raw data: m * estimated with 70 neighbourhood points
FACTORIAL KRIGING ANALYSIS 555
Figure 11
Map of the vertical gradient of reduced to pole field
6. CONCLUSION
REFERENCES
Luc SANDJIVY
ABSTRACT
INTPODUCTION
"harmonics" separately.
Even though the method is always the same, the simplest way
to understand it and to see its usefulness is to begin with a uni-
variate study in the stationary cace.
A UNIVARIATE STUDY.
Theory
Z(x) = L: a Y (x)
u u
u
where the Yu(x) are weakly stationary random functions which are
mutually orthogonal:
C(h) = ~ a 2 C (h)
u u
u
L; A = 0 when m is unknown
a a ~ : Lagrange coefficient
or
I~ L;S AQ
\
f'
c(a,S) = a C (a,x)
U U
when m is known.
If x E V Z*(x) = L; a Y* (x) + m*
u u
u
where
Z* : kriged value of Z at D oint x
Yi* x) .:
co }rigj:' d value of Yu at D oint x
m : known e xpectatiQn E( Z( x» or its krig= d estimate in V.
562 L. SANDJIVY
,------------ - -C2(h,
Co ,-7'------------C1(h'
','
I,
x x
Figure 1
Theory :
In multivariate analysis, the general procedure is the same
but is not so simple to apply since a number of (cross)-variograms
have to be computed and modelled and large co-kriging programs
used. This has not been done so far, but there is a very interest-
ing case that leads back to the univariate one ; that is, when
the variables have an intrinsic correlation. Then, all the (cross)-
variograms are proportional to the same variogram function.
factor to eigenvalue v ..
l
THE FACTORIAL KRIGING ANALYSIS OF REGIONALIZED DATA 565
F.(x) =
l
6 b
u
Y (x)
u
decomposition of F.(x) into the
l
u
different Y (x) corresponding to the C (h)
u u
G.(h)
l
= 6 b 2 C (h)
u U
covariance of F.(x)
l
u
Case Study
0")
C.. (h) = e .. C(h)
lJ lJ .48
( .481 .63
e .. = 1.
lJ
.42 .63 1.
C(h) = a 2 C t- a~ CI (h) + a 22 C2 (h)
0 0
Eigenvalues v. :
l
CONCLUSION
REFERENCES
(1) Matheron, G., 1982 : "Pour une analyse krigeante des donnees
regionalisees". Note nO N-732 , Centre de Geostatistique,
Fontainebleau.
(2) Sandjivy, L., 1982 : "Etude statistique et structurale de
donnees geochimiques". Note nO N-747, Centre de Geosta-
tistique, Fontainebleau.
(3) Sandjivy, L., 1982 : "Cartographie des donnees geochimiques".
Note nO N-768, Centre de Geostatistique, Fontainebleau.
(4) Journel, A.G. , 1978 : "Mining Geostatistics" , Academic
Press, London.
THE FACTORIAL KRIGING ANALYSIS OF REGIONALIZED DATA 567
PRESENTATION OF RESULTS
The mapping of ANA and FACTl along with SPl and SP2 shows
the local and regional features of the regionalisation. See maps
2,3,4 for Cu and 6,7,8 for FACTl.
F .K.A OF COPPER
Map 2 ANA
Distribution of Cu
Map 3 SPl
Local Variations
, \ t -- \
\ \ \ .) ~,'
,, "" .... J \ ......,:;: ~ _ F
,,
Map 4 SP2
~ARIABLE CUIVRE
***************
NUM PPM ANA FEP SP 1 SP2 Rl R2
169 61.0 1.942 2.228 0.596 0.633 0.914 0. 173 0.367
201 122.0 2.536 3.317 1.586 0.136 0.946 0.320 0.055
235 55.0 1.849 1.807 0.636 0.792 0.832 0.207 0.515
264 63.0 1.959 3.211 0.908 -0.195 0.977 0'.19S -0.084
301 74.0 2.123 2.695 0.920 0.376 0.955 0.230 0.188
318 60.0 1.909 2.343 1.172 0.222 0.935 0.331 0'.126
582 117. It.\' 2.'177 3.315 1.928 -0. 11 ~ 0.924 0 . .:380 -0.043
654 70.0 2.077 2.634 1.297 0.171 0.941 0.328 0.086
679 107.~ 2.'125 2."723 1. '-l18 0.497 0.912 0.336 0.235
683 82.0 2.203 3.<172 0.55'1 0. 112 0.993 0. 112 0.045
750 56.0 1.864 3.510 0.990 -Ii:!. 1371 0.957 e. i91 -0.220
793 89.0 2.300 2.776 0.694 0.559 0.939 0.214 0.268
814 67.0 2.015 1.948 0.584 0.931 0.616 0'. L 73 0.552
834 52.0 1.635 0.723 0.539 1. ')89 0.302 0. l59 0.940
838 82.0 2.233 0.885 1.202 1.656 0.335 0.321 0'.866
839 86.0 2.265 0.829 0.956 1.660 0.292 0.238 0.926
841 51.0 1.622 0.'145 0.656 1.709 0.178 B.186 0.966
850 329.0 2.977 2.347 1.067 1.651 0.691 0.222 0.G60
852 184.0 2.806 2.6'11 0.635 1.08 0. 785 0.133 0.605
857 51.0 1.795 1.120 0.703 1.17S 0.542 0.241 0.605
925 60.0 1.925 3.923 1.269 --0.923 0.927 0'.212 -0.308
980 67.0 2.034 2.778 0.772 0.276 0.972 0'.191 0.137
987 71. 0 2.099 1.678 1.186 0.789 0.603 0.358 0.'177
1001 51.0 1.808 0.530 0.858 1. 531 0.229 0.263 0.937
1002 50.0 1.783 0.496 0.564 1.669 0.203 0.163 0.966
1177 105.0 2.379 2.405 1.049 0.B47 0.663 0.266 0.430
1182 74.0 2. ilJ8 1. 871 1. 013 0.944 0.777 0.298 0.555
1195 651.0 3.054 3.'513 1.119 0.937 0.916 0.206 0.345
1212 63.0 1.977 1.184 1.007 1. 212 0.538 0.323 0.778
1228 161. 0 2.693 2.121 1.093 1.430 0.700 0.255 0.667
1256 56.0 1.878 1.012 0.862 1.279 0.'-'68 0.282 0.837
1293 95.0 2.338 1.309 1.201 1.492 0.'199 0.324 0.804
1296 67.0 2.055 1.396 0. 767 1.283 0.593 0.231 0.771
1325 133.0 2.606 1.425 1. '121 1.647 0.'190 0.345 0.801
1334 63.0 1.995 0.972 0.683 1.547 0.398 0.198 0.896
1337 82.0 2.17S 1.125 0.904 1.560 0.440 0.250 0'.863
1338 59.0 1.894 0.922 0.886 1.349 0.'117 0.283 0.863
TABLE 5
F.K.A. of Copper - Decomposition of sample values> 50 ppm
(See presentation of results for comments).
57C L. SANDJIVY
Map 6 FACTl
Distribution of the
first ' factor issued
from the PCA of Pb ,
Zn, Cu .
Map 7 SPI
Local Variations
,, , ,
I
'. , ,, ,
I
"-' ,---~~, ~
Map 8 SP 2
Regional Variations
THE FACTORIAL KRIGING ANALYSIS OF REGIONALIZED DATA 571
~RRIRBLE FRCTEURI
**-***************
NUH FRCT! PEP SPI SP2 R0 Rl R2
201 ~.2~2 3.266 1. 6~6 0.0~2 0.667 0.~62 0.009
2~5 2.90~ 1.570 8.617 1.220 0.Hl 0.270 0.615
251 3.32~ 1. 619 1. 216 1.121 0.701 0.~69 0.519
26~ 2.761 2.755 1.062 -0.~32 0.938 0.332 -0.156
301 ~.332 2.9:2 1.363 0.867 0.677 0.386 0.265
305 3.210 2.353 1. 016 0.q62 0.910 0.364 0.199
506 2.615 1.175 0. '168 I. 1~0 0.6~0 0.386 0.663
553 2.655 1.132 1.089 1. 01 ~ 0.620 0.51~ 8.594
562 3. 1~~ 2.716 1.904 -0.046 0.639 0.545 -0.816
563 2.717 2.011 1.326 -0.015 0.65~ 0.522 -0.807
6~6 2.691 1.166 I .~~6 0.069 0.601 0.596 0.8~3
TABLE 9
Andrew R. Solow
Stanford University
ABSTRACT
INTRODUCTION
573
G. Verly et al. (eds.) , Geostatistics for Natural Resources Characterization, Part 1,573-585.
© 1984 by D. Reidel Publishing Company.
574 A.R.SOLOW
lim s'l-cq) = S~
The differences between the time domain and the frequency domain
approaches are clear. • The frequency domain approach models a
process as the sum of sinusoids. while the time domain approach
models a process as a parsimonious ARMA process. While both
models are theoretically equally general. the harmonic model
presumes a greater understanding of the physical process. It is
linked to the model of hidden periodicities. The family of ARMA
models. while equally general. is less bound to a particular
philosophy. In terms of structural analysis. the differences are
perhaps even sharper. The spectrum is designed for decomposition
of variance according to frequency. The autocorrelation
function. on the other hand. is designed for the study of the
pattern of correlation.
CONCLUSIONS
LITERATURE
REFERENCES
Jenkins, G.M. and Watts. D.G. (1968). Spectral Analysis and its
Applications. Holden-Day, New York.