Computers & Geosciences, Vol. 3, pp.
173-180. Pergamon Press, 1977. Printed in
PRACTICAL KRIGING IN THREE DIMENSIONS
ISOBEL
Department of Mineral Resources Engineering,
Imperial
(Received 10
September 1976)
Abstract Little has
been published to date on the practical difficulties of ore reserve estimation
in three-dimensional deposits. Some
authors have suggested condensing the problem into two dimensions, but this is
not always practicable or desirable. A
suggestion and a FORTRAN IV Function Segment are provided which may alleviate
some of these difficulties.
Key Words: FORTRAN,
Geostatistics, Kriging, Regionalized variables, Semivariogram.
INTRODUCTION
At the feasibility study stage
of a potentially mineable deposit, the ore-reserve valuation generally is
carried out using data from (usually vertical) boreholes drilled through the
deposit. Ideally for economic appraisal
of the site, local as well as global estimates of ore grades and tonnages must
be made. Where a deposit is of a seam or
vein type, or where it comprises a sedimentary layer, the problem is
essentially a two-dimensional one. The
width (or depth or thickness) of the deposit, and the accumulation (grade times
width) of metal may be considered as measurements or point samples in a
two-dimensional plane (Sinclair and Deraisme, 1974). There may be some 'unrolling' problems, but
essentially these are not statistical.
However, in a three-dimensional deposit new problems occur. If local estimates are to be produced, these
must be block estimates on a bench by bench basis. Many such applications are to be recognized,
this author having experience with uranium, nickel, and lead/zinc deposits of
disseminated character. In these,
two-dimensional simplifications are not always possible. It is intended to use one of these deposits
as an example throughout this paper. It
is not possible to describe the deposit in full but for the purposes of this
paper, it may be said that we have a disseminated nickel deposit in northern
BLOCK
ESTIMATION
To estimate each block by the
method of Kriging, we must first set up the Kriging system of equations. These equations require that we must evaluate
the average semivariogram between each of the samples we wish to include in the
estimation procedure and the unknown area.
Also we require the values of the average semivariogram between each
pair of samples, and between all points within the unknown area. In one or two dimensions these evaluations
present little or no problem.
In three, however, the problems can be great. Two approaches have been suggested, and this
author ventures to suggest a third.
David (1976) suggests that a deposit may be considered bench
by bench. That is, we ignore all data
above and below the bench under consideration.
Thus, both information and data span only the bench width (see Fig. 1)
and the problem may be reduced to a two-dimensional one. Figure 2 shows the bench as if from
above. The average grade of the borehole
intersection then may be considered as coming from a point sample, and the
block is reduced to a panel. This
approach has many advantages as far as computation goes. Two-dimensional auxiliary functions (cf.
There are disadvantages to this approach, as illustrated in
Figure 3. Boreholes may not make intersections with the full bench width. For example: boreholes may stop before the
bottom of the bench (a); assaying may not have started until part way into the
bench; (b) there may be core loss, or unassayed sections in the hole; (c) holes
may be inclined, so that the intersection of the borehole with the bench is
longer than the actual bench width. If
these problems do not occur or it is considered an adequate approximation, then
David's approach may be used. However,
there is one other situation which cannot be handled by this approach. Suppose, as in our example, we have a bench
width which is smaller than the range of influence of the semivariogram. Imagine a situation in which a borehole goes
through, or close to a block to be estimated (Fig. 4). It is intuitively more reasonable to expect
that the portions of borehole (i) above
and below the block should have a closer relationship to, and hence more
influence on, the shaded block than the portion of borehole (ii), say, on the
same bench. That is, we would like to
include sections of the boreholes above and below the bench in the Kriging
system, because this should improve the estimation of the block grade
considerably.
How much information above and below the block we wish to
include in the estimation would depend on the bench width, size of the block,
range of influence, and also on how the core is to be weighted in itself. Ideally, each core section which was assayed
should be treated as a separate sample, and a Kriging system produced to determine
the weight for each section. In that
situation, obviously, the more sections included the lower the estimation
variance. However, after a certain
distance the increase in accuracy must become negligible-especially compared to
the amount of work required to evaluate it.
Alternatively, the core length could be split into the 'internal' (to
the bench) portion, and the 'external' portions and each of these given an
individual weight. This may become
complicated, but might be worth doing if the range of influence is not greater
than the bench size. A third possibility
is to weight the whole length of core equally, that is average the length to be
included in the estimation for each core and accord each one weight in the
Kriging system. This approach is favored
by the author in situations such as the Norwegian nickel example, because it
forms a workable compromise between the ideal situation and the practicalities
of computing in a reasonable size 3-D deposit.
However, comparative studies remain to be carried out in any case study
to ascertain the 'best' approach for the particular situation. Having decided upon this approach, it remains
to evaluate the amount of core which yields the 'best' estimators. This could be determined by calculating the
extension variance of a core of different length to a block of the desired
size. However, this returns us to the
problem of calculating semivariograms in three dimensions.
An alternative approach to David's is necessary once we have
decided to work in three dimensions.
Because 3-D auxiliary functions are intractable, we must produce
numerical approximations to the required block-core, core-core, and block-block
semivariograms. If cores are parallel
and of the same length, core-core relationships can be evaluated by the 2-D
auxiliary function g(
; b). Because our decision to work in 3-D was
prompted partly by the fact that the cores are probably not of the same length, this simplification is unlikely to
occur. The simplest numerical solution
to calculating the average semivariograms between such pairs is to approximate
each volume or length by a close network of points. A. Marechal (1975, pers. comm.) has stated
that a block may be considered as a 4 x 4 x 4 network of points, and that this
will produce an accuracy of 1 percent in the final semivariogram. Similarly a core might be considered to be a
string of discrete points-presumably 64 points along the length would give
similar accuracy. Then the average
semivariogram between core and block would be calculated by evaluating the
semivariogram between every point in the block and each point of the core (642
combinations) and then dividing by the number of combinations. Each average semivariogram would need 4096
calls to the point semivariogram in its calculation. This approach is heavy on computer time and
(the author feels) uncertain of the final accuracy of both approximated
semivariograms and final estimates.
A third approach has been suggested by this author elsewhere (; b) may be
used here. Calculation of core-core or
core-block semivariograms where one core is not parallel to another present a
difficult problem - and are as yet unsolved.
However, for boreholes which may be considered to be vertical (sic) the
calculation of core-core relationships presents no problems at all. The author presents in Appendix I a FORTRAN
IV Function Segment which will calculate the average semivariogram between any
two cores of any lengths, in any relative positions and any distance apart so
long as they are parallel, for a spherical model.
The function is called GIMEL (the Hebrew letter G) and is used
as follows (see Fig. 5):
GAM = GIMEL (EL, B, H, D, A, C),
where A is the range of
influence and C the sill of the point semivariogram (spherical type); GAM is to
contain the average semivariogram between two parallel cores of length (EL) and b (B), a distance h (H) apart, whose ends are offset by a length d (D) as illustrated.
Restrictions on the parameters are: and b must be positive, b must not be greater than . If the user finds this restrictive, a simple insertion
of about six lines will generalize the FUNCTION. d may take any value. h
may be positive or zero. This last
enables different lengths on one borehole to be compared.
OTHER USES OF GIMEL
Although the use of the
core-core function reduces the calculation of bench/block estimations
considerably in the regular situation, this is not the only situation in which
it might be useful. Approximating one of
the three dimensions by a continuous line allows greater attention to be paid
to irregularities which may occur. For
example, in
AN EXAMPLE
Producing simple tables for
checking GIMEL is difficult because of the four arguments. However, a program example is given in
Appendix II with its output, to enable users to check the FUNCTION. Some checks also may be made against the
standard function g(; b) by
computing GIMEL (b, b, , 0, a,
c). Also GIMEL (, , 0, h, a, c)
may be compared with regularized semivariogram for length .
The example used is a program to determine the extension
variance of a vertical core section of length , at a distance h
from the center of a block along a median through the block (see Fig. 7). The disseminated nickel example has been
used, that is a = 50, c = 1.0, blocks
are 25 m by 25 m, benches are 10 m. A grid of 8 by 8 'cores' has been used to
approximate the block. Values of h range from zero, sample central to
block, to 50 m, that is sample almost out of range. Lengths of the core range from 2 m to 50 m
and are centered on the center of the bench.
Also it may be noticed, the program takes advantage of any symmetry
which might be present.
Acknowledgments - The function GIMEL was first conceived and initially
evaluated while the author was employed by the Norwegian Geological Survey in
1974. The programs described in this
paper were produced when the author was Visiting Professor of Geology at
REFERENCES
David, M., 1976, The practice of Kriging, in Advanced geostatistics in the mining
industry: D. Reidel, Dordrecht-Holland, p. 31-48.
Sinclair, A. J., and Deraisme, J., 1974, A
geostatistical study of the Eagle copper vein, northern
Please note that, due to the
inefficiency of the OCR program,
Appendices have been omitted
from this copy.
If I cannot guarantee that
it works, I won't hand it out.
If you really want a copy of
35 year old Fortran code,
please e-mail geoecosse@kriging.com