ORE RESERVE ESTIMATION USING RADIAL BASIS FUNCTIONS

This paper presents an application of multiquadric equations for ore reserve estimation. This method is compared with available traditional estimators used in the mining industry: inverse distance weighting and ordinary kriging. Radial basis function with a constant precision term has a dual form absolutely equivalent to the ordinary kriging. In this sense, radial basis functions present a potential application for ore reserve estimation. Besides, radial basis functions can be adapted easily to allow block estimates as block kriging does. It is shown that radial basis functions depend on the metrics of data coordinates and so the basis functions should be carefully chosen. This paper also discusses kernel transformation for improving radial basis function estimates. Using the exploration data from the Chapada Copper Deposit, it is shown that logarithmic transformation of basis functions provides the best results comparable to the ordinary kriging ones. Thus, radial basis functions can be applied for ore reserve estimation, becoming a reliable alternative to this task.


INTRODUCTION
Approximation of a variable at unsampled points is a common procedure used for computer aided ore reserve estimation.Usually, a mineral deposit is discretized in small blocks of dimensions compatible with the available exploration data.The discretization is made within recognized boundaries defined by the exploration data.A certain number of blocks have to be estimated and their individual ore reserves determined.The estimation or interpolation procedure is applied to determine the reserve parameters for a block using neighboring exploration data.There are basically two interpolation methods that have been used for computer aided ore reserve estimation: inverse distance weighting and ordinary kriging.The first method is indicated when available data do not justify variogram computation and so the ordinary kriging technique cannot be applied.The second one has been used extensively by the mining industry for ore reserve estimation.However, it is important to mention that interpolation for threedimensional data is not limited to these methods.There are several available methods that could be applied for ore reserve estimation.One of these methods is derived from radial basis functions, which allows approximation of any irregular surface by means of a sum of a number of kernel surfaces.Radial basis functions are a generalization of multiquadric equations, which were firstly suggested by R.L. Hardy for the analytical representation of a terrain surface from discrete data points (HARDY, 1971).Although radial basis functions have been developed for interpolation of two-dimensional data, it can be extended for three-dimensional data as well as for higher dimensions.This paper presents an application of radial basis functions for ore reserve estimation.

ORE RESERVE ESTIMATION METHODS
As introduced before the current methods for ore reserve estimation are: inverse distance weighting and ordinary kriging.
The inverse distance weighting was among the first computer-based techniques used in ore deposit calculations (PHILIP & WATSON, 1987).
Kriging is the standard name given to the collection of generalized linear regression techniques for minimizing an estimation variance defined from a prior covariance model (OLEA, 1991).There are many different types of kriging methods depending on objectives, conditions and constraints for each type.So, kriging with a known mean receives the name of simple kriging (SK), kriging for the mean value receives the name of kriging of the mean (KM), kriging with an estimated mean receives the name of ordinary kriging (OK), among other types of available kriging methods.Details for these different types of kriging can be found in (WACKERNAGEL, 1995).In this paper only the ordinary kriging will be considered because its formulation allows solving most of the ore reserve estimation problems.

Inverse distance weighting
The inverse distance weighting is an interpolation technique for estimating values of locationally dependent variables by forming a linear combination of a set of measurements (PHILIP & WATSON, 1987).The non-negative weights sum up to one and are inversely related to the distance to data points (PHILIP & WATSON, 1987).The grade or any other variable at unsampled points can be estimated by using the inverse distance weighting as: where d i is the distance from each of the n' sample locations to the point to be estimated, p is the power of distance, and F 1 , ..., F n' are the grade values.Therefore, the normalized weight for the i-th sample becomes: FIGURE 1 -Block estimation by extending the grade computed in its center.
For ore reserve estimation this technique was used for interpolating a point in the center of an evaluation block (Figure 1).The resulting grade was extended to the whole block with an associated error proportional to the block size.
This author (YAMAMOTO, 1996) has proposed a new method of block calculation using the inverse distance weighting method, according to the following equation: where i W is the average of the inverse distance weighting of the i-th sample in relation to all subblocks.This new method is an adaptation of the conventional block kriging method, where a block estimate is equal to the average of point estimates centered on sub-blocks.So, given a block it is discretized into nsb sub-blocks of equal size.For block discretization the limits given by JOURNEL & HUIJBREGTS (1978) are recommended (Table 1).
Normalized weights between samples and each sub-block are computed as illustrated in Figure 2.
The normalized weight between i-th sample and k-th sub-block is: FIGURE 2 -Computation of normalized weights for each sub-block: A) for the sub-block 1; B) for the sub-block 2; C) for the sub-block 3; D) for the sub-block 4.
After the normalized weights for all subblocks have been computed, the average weight for each sample is computed as follows: where nsb is the number of sub-blocks.Expression (5) allows computing the grade of a block directly as does ordinary kriging.Obviously, it gives the same result as the average of grades of nsb sub-blocks.However, the great advantage of expression ( 5) is the direct calculation of an estimation error of inverse distance weighting block estimates, by means of the interpolation variance as proposed by this author (YAMAMOTO, 2000) for ordinary kriging estimates.The interpolation variance of a block can be computed as: ( ) Actually, this is an extension of the same formula used for computing a measure of the reliability of ordinary kriging estimates.This extension is valid for inverse distance weighting because their weights are positive and are normalized, which assures the positiveness of the interpolation variance.

Ordinary kriging
Kriging depends firstly on a structural function C(h) or γ(h), which describes the spatial correlation of data points (Figure 3).
Ordinary kriging (OK) is the most widely used kriging method and is based on local second-order stationarity (JOURNEL & HUIJBREGTS, 1978).Ordinary kriging is by excellence the recommended method for ore reserve estimation, because it was designed for block evaluation.Indeed, ordinary kriging was the first method that provided block computation.
OK estimate is based on weighted average of available data: where the weights {λ i , i=1,n'} are determined by the solution of a kriging system of equations.
In order to ensure unbiasedness of estimates, the expected value of the error which can be developed as: which is known as the non-bias condition (JOURNEL & HUIJBREGTS, 1978).
On the other hand, ordinary kriging makes estimation under the condition that the error variance is minimum.The error variance is: which can be written after ISAAKS & SRIVASTAVA (1989) as: Deriving each term on the right side of (9), we have: Minimization of the error variance (10) constrained to the unbiasedness condition (7) results in the following system of ordinary kriging equations: where µ is the Lagrange multiplier.
In the case of block kriging the right side: − , i.e. the average covariance between the i-th point and centers of sub-blocks.
Figure 4 illustrates the process of converting distances to covariances.Figure 4A shows how to compute the covariances between sample 1 and all others.Repeating this process for all samples we will take into account the mutual dependence among samples.This feature recognizes spatial clusters and assigns correct weights to the samples.In Figure 4B distances between samples and the center of the block are converted to covariances (point kriging).
The set of weights {λ i , i=1,n'} results from the solution of the ordinary kriging system (10).These weights applied to expression (6) give the estimated value.The error associated with the estimated value is given by: See that the kriging variance is homoscedastic, i.e. it is independent of the data used to obtain the estimator F*(x o ) (OLEA, 1991).Figure 5 illustrates why kriging variance should not be used as an approximation of uncertainty (ARMSTRONG, 1994).
As the data layouts are identical in both cases, the kriging variances are identical and as it happens, so are the kriged estimates (ARMSTRONG, 1994).

RADIAL BASIS FUNCTIONS
Radial basis functions are a generalization of the original multiquadric equations (HARDY, 1971).The basic hypothesis of the multiquadric analysis is that any smooth mathematical surface, and also any smooth arbitrary surface (mathematically undefined) may be approximated to any desired degree of exactness by the summation of a wide variety of regular, mathematically defined surfaces, particularly quadric forms (HARDY, 1977).The quadric forms are not only the simplest, but also the most efficient in converging on irregular surfaces (HARDY, 1977).Now, the quadric forms are in fact one of the available kernels of radial basis functions.It is important to mention that although radial basis functions were defined on global basis, local approximation can give reliable results.In this regard, accuracy of multiquadrics can be improved and computational effort can be dramatically reduced by transforming a large global problem into many small quase-local problems by domain decomposition (KANSA, 1990).
A radial basis function with centers in {x 1 ,x 2 ,. . .,x n' } as given in BEATSON & NEWSAM (1992) is represented as: where n' is the number of neighboring points and f is a symmetric radial function in R n .Indeed, this is the form originally proposed by HARDY (1971), but a more general form has been proposed by other authors (CARLSON & FOLEY 1992, GIROSI 1992, MYERS 1992), which includes a second summation for equation ( 12) as: φ where p j (x), j=0,..,k are the polynomial terms in the position coordinates of x.Solution for this system of equations requires the following set of constraints: Addition of polynomial terms does not mean a precision improvement, because it may degrade the accuracy of the interpolant on some data sets (CARLSON & FOLEY, 1992).Thus, some authors (MICCHELLI 1986, MADYCH 1992) have used the expression (12) with just a term of degree 0 added: where a o is also known as the constant term.Addition of the constant term requires the constant precision constraint (FOLEY, 1992): Addition of a constant term improves the precision of radial basis function, especially when n' (number of neighboring data-values) is small.Actually, a o represents the local mean of data values as implicit in ordinary kriging (WACKERNAGEL, 1995).Therefore, radial basis function when represented as ( 13) is very similar to ordinary kriging, as it will be shown hereafter.
There are several options for basis function ( ) φ , as described by KANSA (1990) and  ) 16) in ( 15) we have: which is known in geostatistics as the dual system (WACKERNAGEL, 1995).In fact the dual system has been introduced by HARDY (1977) to compare covariance and multiquadric methods.Now weights of radial basis functions in the dual form can be expressed as: (18) where µ is a constant introduced so that the equality becomes possible.
For the generalized multiquadric, when k is -1 we have the reciprocal multiquadric kernel and when k is equal to zero we obtain the original multiquadric kernel (HARDY, 1971).
Among available radial basis functions the multiquadric kernel is the most widely used, because it always guarantees a non-singular matrix of coefficients.Splines and cubic kernels call for coordinate transformation (usually normalization) because they rise rapidly with x and so they are not convenient for ore reserve estimation.In this regard, the performance of multiquadrics is sensitive to scaling, i.e. it is important to the condition number whether distances are expressed in centimeters or meters (KANSA, 1990).Now we can show that radial basis functions expressed as (13) are equivalent to the ordinary kriging.Expression (13) can also be written as: The coefficients {c 1 ,c 2 ,...,c n' } and the constant term a o result from the following system of linear equations: Transposing both sides of equation ( 19) it becomes: or ( ) Therefore, radial basis function estimator when represented as dual form ( 17) is similar to the ordinary kriging estimator (6), as well as their corresponding systems of linear equations ( 20) and ( 10).Thus, the introduction of a constant term to the radial basis function implies in the following condition: that is no other than the unbiasedness condition suggested by SIRAYANONE (1988).
Radial basis functions as well as the inverse distance weighting depend only on distance measurements and so they are easily generalized to higher dimensional spaces (BARNHILL, 1993).Therefore, these methods can be adapted directly to solve problems in ore reserve estimation, where the data-values are usually in the three-dimensional space.Ordinary kriging calls for a three-dimensional variogram.
See that the system of equations ( 20) corresponds to a point estimator rather than to an estimator of a spatial average such as the average grade of a block, but the extension to spatial averages is relatively easy in the kriging form of the estimator (MYERS, 1992).So, for block estimation using radial basis functions, replace the right side of the system of linear equations (20) by the average basis values, as follows: As a case study it was considered the exploration data from the Chapada Copper Deposit.This deposit is located at municipality of Mara Rosa, State of Goiás, Brazil, according to the location map of Figure 6.
The exploration data, considered in this paper, come from 141 drill holes (Figure 7) which totalizes 16315.80meters of continuous sampling.A total of 10504 samples analysed for copper were composited to a bench height of 10 meters resulting in 1596 composite samples and they will be the database for this study.Figure 8 shows copper grade distribution with its statistics.
So, given the boundaries of mineralization (Figure 7), the deposit is discretized into blocks of compatible dimensions with the exploration data, defined to be 100 x 100 x 10 m (width x length x height).The block model for the Chapada Copper Deposit is shown in Figure 9.
Geostatistical analysis carried out using copper exploration data resulted in the following covariogram model: Geometrical anisotropy (parallel to the coordinate axes) was recognized according to the anisotropy ratios of Table 2. Now we have the necessary elements to perform grade estimation in blocks or at points by ordinary kriging and make comparisons with radial basis functions.

Choosing a basis function
The first step for radial basis function interpolation is related to choosing an adequate kernel or basis function.This step corresponds to the variogram modeling for ordinary kriging estimation.Among several available options for basis functions we have to choose one, which could give unbiased results.Moreover, in geology we also have to consider the geological correlation between samples, i.e. there is a maximum distance within which samples present a correlation.Beyond this maximum distance samples are not correlated, when changes in composition, lithology and other physical and chemical properties usually occur and, therefore, cannot be predicted.Therefore, a maximum distance between samples and the point to be interpolated should be taken into account.For the case study, as we have a covariogram model, a maximum distance equal to 450 m (corresponding to the variogram range) will be used.
Our main concern when choosing a basis function is related to the metrics involved, because distances between samples are measured in meters (tens to hundreds of meters).So, some basis functions like cubic and splines are not convenient, because the rounding errors in the floating-point arithmetic can lead to innacurate results.
The cross validation method was adopted as a standard testing procedure to choose a basis function.This procedure consists of removing each  sample from the database in turn and estimating the value at that location using n' neighboring points.So, for each point we have the actual: F(x i ) and estimated: F*(x i ) values which can be compared, giving us an indication of how well the basis function fits into the data.Some statistics can be computed from the difference between true and estimated values such as: correlation coefficient and RMS error.RMS error (MADYSH, 1992) was computed using: which measures the dispersion of estimated values around true ones.Another way for evaluating basis function can be done by determining the condition number of a matrix of the system of equations ( 20).The condition number of matrix tells us how invertible the matrix is.If the condition number is infinite the matrix is singular, and ill-conditioned if it is too large (PRESS et al., 1989).A procedure available in PRESS et al. (1989) for singular value decomposition was used to compute the condition number of a matrix.
So, several simulations were carried out for the Chapada database in order to choose a multiquadric constant (c), which would give the minimum RMS error.In these simulations c values ranging from 1 to 10000 were used.Figure 10A shows RMS error as a function of the multiquadric constant c.In such figure we observe that the multiquadric constant should be chosen as small as possible.It is important to mention that the multiquadric constant only causes round off errors and it practically does not affect the condition number of the matrix as shown in Figure 10B, with magnitudes ranging from 10 5 to 10 6 .In Figure 10 both RMS error and condition number are represented as their averages coming from 1596 interpolated points.
According to obtained results from cross validation tests, the multiquadric constant must be as small as possible.Therefore, we have adopted a constant equal to 1, which best fits to the exploration data from Chapada.Actually, we will adopt hereafter in this paper the translated linear kernel:   ( ) . Now we are interested in comparative tests between different interpolation techniques: inverse distance weighting, ordinary kriging and radial basis function with the translated linear kernel.In case of the inverse distance weighting a power of distance equal to 2 was considered.Table 3 presents the statistics for cross validation tests.
Scattergrams for the cross validation results are presented in Figure 11.
The best result is achieved with ordinary kriging, followed by radial basis functions and finally by inverse distance weighting.When using radial basis functions the metrics involved is important and it will govern the accuracy of the estimates.So, in order to improve the estimates using radial basis functions some transformations on the translated   linear kernel were considered.These transforms are as follows: Square root: ( ) Transformed kernels are illustrated in Figure 12.Table 4 presents computed statistics on scattergrams as well as the average condition number for estimates using radial basis functions with transformed kernels.As we can see, transformed kernels present average condition number much lesser than the multiquadric kernel.It means that besides the inversion problem on matrix of the system of equations ( 20), the round off errors as seen in Figure 10A can be avoided.
Figure 13 presents scattergrams of cross validation tests using radial basis functions with transformed kernels.
Indeed transformed kernels improve estimates as shown in scattergrams of Figure 13, because the accuracy of the radial basis functions depends strongly on the metrics used.Moreover, transformed kernels effectively diminish the condition number of matrix (Table 4).So, as we can see it is nonsense to consider other kernels like cubic and splines, because they rise rapidly with x .Therefore, all the transformations were done over the translated linear kernel.See that the transformed kernels limit the influence of far samples as well as the variogram does.
Considering that ordinary kriging gives the best results, let us compare it with inverse distance weighting and radial basis functions with no transformation (Figure 14).
Figure 14 shows that both methods provide biased estimates in relation to the ordinary kriging.Therefore, they cannot be used for reliable ore reserve estimation.Now let us verify if the transformed kernels provide better results than kernels with no transformation.Figure 15 presents scattergrams of ordinary kriging estimates against radial basis functions using transformed kernels.
As we can see square root is not an effective transformation for the database of the study case, because of the metrics involved (hundreds of meters).On the other hand, natural and decimal logarithms are the most effective giving the best results.Both give almost unbiased and basically the same results, which proves that kernel transformation provides good estimates comparable to the ordinary kriging results.The dispersion around the bisector may be due to the differences between spherical covariance (kernel for ordinary kriging) and transformed kernel (natural or decimal logarithm).Moreover, there is also the correction for geometrical anisotropy done by ordinary kriging while radial basis functions do not.
Radial basis functions with appropriate kernels give reliable results being comparable to the ordinary kriging ones.Hence, radial basis function is another method available for ore reserve estimation.

CONCLUDING REMARKS
Radial basis function estimator with a constant precision term has a dual form absolutely equivalent to the ordinary kriging estimator.Hence, radial basis function arises as an alternative method to ordinary kriging, especially in cases where experimental variograms cannot be computed from data points.This method is much better than the inverse distance weighting, which has been used as an alternative to the ordinary kriging.Radial basis function is superior to the inverse distance weighting, because it can recognize clustered points giving adequate weights according to their spatial configuration, such as done by ordinary kriging.Accurate results from radial basis function estimates depend strongly on metrics involved and also on kernels.The linear kernel with a logarithmic transformation gives the best result for exploration data from the Chapada Copper Deposit, as confirmed by average condition number of matrices.

FIGURE 4 -FIGURE 5 -
FIGURE 4 -Converting distances to covariances: A) between sample 1 and all others and B) between samples and the block center.

FIGURE 6 -
FIGURE 6 -Location map of the Chapada Copper Deposit.

FIGURE 7 -
FIGURE 7 -Location map of drill holes from Chapada Copper Deposit and boundaries of mineralization.

FIGURE 9 -
FIGURE 9 -Three-dimensional model of blocks for the Chapada Copper Deposit and the block 3176.

FIGURE 10 -
FIGURE 10 -Influence of the multiquadric constant over estimates as measured by RMS error (A) and over condition number of matrices (B).

FIGURE 15 -
FIGURE 15 -Comparison of ordinary kriging estimates against radial basis functions using transformations by: square root (A), natural logarithm (B) and decimal logarithm (C).

TABLE 2 -
Anisotropy ratios for the copper semivariograms.

TABLE 3 -
Statistics of cross validation tests for the Chapada data base.

TABLE 4 -
Statistics of cross validation tests using radial basis functions with transformed kernels.