1. Introduction
Mineral deposits such as iron ore, bauxite, and phosphate are characterized by containing, in addition to the main elements (Fe_{2}O_{3}, Al_{2}O_{3}, P_{2}O_{5}, etc.), other elements or chemical species with effects on economic viability, industrial processes, or mine planning. It is common to estimate multiple elements, possibly correlated and sometimes with a combination of contents that must sum to a particular figure (e.g., 100%). According to ^{Aitchison (1981)}, data that add up to a constant are termed compositional data (CODA), and they carry information that is relative and not absolute. This condition of summation to a constant implies that the estimates should also sum to a constant.
When working with multi-element deposits and, furthermore, having to deal with CODA, which are not necessarily physically correlated (spurious correlation; ^{Pearson, 1897}), it is not possible to use traditional methods to achieve closure of the mass balance of the multiple chemical species or physical variables. Thus, to overcome this inconsistency, it is typically necessary to perform post-processing, such as proportional distribution of the error of closure between the different granulometric fractions for each of the estimated elements.
Classical geostatistical methods such as ordinary kriging (^{Matheron, 1963}) may be appropriate for the best local estimate of a single variable, ignoring its spatial interdependence with other correlated attributes. Each variable is estimated separately (in the case of ordinary kriging) with its specific parameters of spatial continuity, which leads to different weights being obtained for each attribute and a failure to obtain estimates that satisfy the constant-sum constraint. In the case of ordinary cokriging (^{Marechal, 1970}), which takes into consideration the correlation between multiple variables, closure can only be ensured when working with an intrinsic coregionalization model (ICM). However, the model used as reference for the direct and cross variograms is rarely adjusted adequately for all variables. When a linear coregionalization model (LCM) is used, the complexity of modeling the variogram increases with the number of variables and also fails to ensure the closure balance.
^{Aitchison (1986)} developed two transformations to deal with CODA, ensuring that any operation applied to the transformed data sums to a constant after these data are back-transformed to the original space. These transformations are known as additive logratio transformations (alr) and centered logratio transformations (clr). ^{Egozcue et al. (2003)} defined new transformations, called isometric logratio transformations (ilr), which are used in this article. A fundamental feature of the methods mentioned is the use of nonlinear transformations (logarithms).
Pawlowsky et al. (1995) and ^{Odeh et al. (2003)} applied alr cokriging and univariate ordinary kriging, respectively, to predict composition at unsampled locations. They used the inverse transformations (agl) to back-transform the elements of the estimated variates. However, this back-transform is biased (the average data transformed by a nonlinear function cannot be back-transformed by a linear function (OK) without generating a bias) and a solution for an unbiased back-transform is unknown (^{Pawlowsky-Glahn and Olea, 2004}).
^{Bragulat et al. (2002)}, ^{Bragulat and Sala (2003)}, and ^{Boezio et al. (2012)} used kriging and cokriging of logratio transformations applied to mineral deposits. They also used the inverse transformations in the estimated variables, but a problem appears when this type of transformation is used in the estimation process, since the average kriged block values cannot be back-transformed without biasing the estimated grades. To solve this problem, ^{Pawlowsky-Glahn and Olea (2004)} suggested a numerical approximation to generate unbiased estimates in the inverse transformations of CODA. This approximation is obtained through the use of the Gauss-Hermite procedure (^{Lark and Bishop, 2007}; ^{Ward and Muller, 2012}; Delgado et al., 2012).
This paper presents an alternative way to deal with transformed data (ilr) and avoid bias in the back-transformation, i.e., geostatistical simulations. The main idea is to back-transform simulated points on a closely spaced grid (instead of estimated blocks) to the original space (at point support), postponing the averaging into larger volumes (blocks) to the end of the process. Furthermore, it is proposed that closure of the sum on chemical and granulometric variables be ensured at each simulated block, thereby reproducing the correlations between them. A complete case study of a major bauxite deposit illustrates the methodology.
2. Methodology
2.1 Compositional data analysis (CODA)
A composition of D parts is a vector x=[x_{1},x_{2},...,x_{D}] all of whose components are strictly positive numbers and carry only relative information. This information is conditioned to sum of a constant and represents parts of a whole, for example, unit (1), percent (100%), or parts per million (ppm). ^{Pawlowsky-Glahn and Buccianti (2011)} define the sample space containing the compositional data as the D-simplex.
where the components of each vector in S^{D} are called the parts of the composition. The operation that defines the closure of a composition in a constant k is given by
where
C(Z) is the closure operation;
k is the closure constant (generally 100%);
Z_{i} is the value of the ith sample.
2.2 Isometric logratio transformation (ilr)
Before defining the transformation (ilr), it is necessary to understand the concept of an orthonormal basis. As in any Euclidean space, there are an infinite number of orthonormal bases in S^{D} that can be obtained by various methods, for example, the Gram-Schmidt procedure mentioned by ^{Egozcue et al. (2003)} or the singular value decomposition (SVD) procedure described by ^{Pawlowsky-Glahn et al. (2010)}. Pawlowsky et al. (2005) proposed a new method for obtaining an orthonormal basis, known as sequential binary partition (SBP).
The SBP is defined by ^{Egozcue et al. (2005)} as a hierarchy of parts of a composition for obtaining particular orthonormal coordinates. In the first order of the hierarchy, all parts are divided into two binary groups (+1 and -1). In the following steps, each group is divided into two new groups, and the process continues until all groups have a single part. The number of binary partitions at the end of the process is D - 1 (where D is the number of dimensions, corresponding to the number of variables per fraction). Table 1 shows an example of the SBP applied to a composition of five parts.
Order | P_{1} | P_{2} | P_{3} | P_{4} | P_{5} | r(+) | s(−) |
---|---|---|---|---|---|---|---|
1 | +1 | +1 | +1 | −1 | −1 | 3 | 2 |
2 | +1 | +1 | −1 | 0 | 0 | 2 | 1 |
3 | +1 | −1 | 0 | 0 | 0 | 1 | 1 |
4 | 0 | 0 | 0 | +1 | −1 | 1 | 1 |
As proposed by ^{Egozcue et al. (2003)}, the isometric logratio transformation of the ith composition is defined by
where
ilr_{i} = ilr transformation for ith composition;
r = sum of the positive 1's (+1) in the SBP;
s = sum of the negative 1's (-1) in the SBP;
(x_{i1} x_{i2}...x_{ir})^{1⁄r} = geometric mean of the variables that were selected with (+1) in the SBP;
(x_{j}_{1} x_{j}_{2}...x_{j}_{s})^{1⁄s} = geometric mean of the variables that were selected with (-1) in the SBP.
This new transformation will have D - 1 dimensions for each composition analyzed, depending on the number of original variables.
The next step consists in using geostatistical simulation methods for each transformation (ilr). In this specific study, the turning bands algorithm (^{Matheron, 1973}) was used to run simulations in multi-Gaussian space. Various alternative simulation methods are available in literature (^{Deutsch and Journel, 1998}), but the one chosen here proved to be efficient for the purpose of the study. At the end of the process, each simulation is back-transformed to the space of the original data by an inverse isometric logratio transformation given by
where
ilr^{-1} = back-transformation;
x = simulated value for the transformation (ilr);
ψ = matrix constructed from the SBP;
C = closure operation (equation (2)).
The construction of the matrix ψ is based on the SBP that was initially defined. Each partition will have its own matrix depending on the number of variables. This new matrix is calculated as follows:
where ψ_{i+} and ψ_{i-} represent the values of the matrix ψ defined as +1 and -1 in the SBP, and r_{i} and s_{i} represent the sums of +1 and -1 obtained in the same partition. For the example presented in Table 1 corresponding to an SBP of five parts, the matrix ψ is defined in Table 2.
2.3 Case study
The case study corresponds to a data set from a bauxite deposit located in the Brazilian Amazon (Figure 1).
The variables correspond to three granulometric fractions (percentages of the total mass retained at given sieves during screening tests). These variables are defined as recoveries at the following fractions: +14# (REC14), +400# (REC400), and -400# (REC-400). Each variable is defined as the percentage of the mass retained on each sieve, and the sum of the variables for each analyzed sample should be 100%. However, there are some errors associated with sampling (^{Abzalov, 2011}) that prevent that the sum of the variables analyzed from closing to a constant. In this particular case, these errors were not greater than ±3%. Therefore, to start the analysis of the CODA, the closure operation given in equation (2) was applied.
The isometric logratio transformation was subsequently applied for each of the compositions of the three analyzed variables. This transformation led to a two-dimensional sample space, in which the variables were called ilr_{1} and ilr_{2}. Each variable was independently simulated, considering their spatial continuity models and search parameters. The total number of simulations was 30 for each variable (number of realizations sufficient to map uncertainty due to the standardization of the variance of the means). The final estimate was taken to be the E-type (average) of these 30 simulations. Figure 2 shows a suitable procedure for working with CODA without generating bias by back-transforming block estimations. Note that the average of the simulated blocks (50x50x0.5)m is obtained after the punctual simulations (10x10x0.5)m are back-transformed by the inverse function ilr^{-1} (step 7).
3. Results
Each simulation generated was validated by the reproduction of the basic statistics of the original data, the variogram model, and the correlations among the variables. Table 3 shows a statistical summary of the original data and the upscaled results for two realizations selected randomly (Nos. 4 and 24). Note that these simulations, like the others, satisfactorily reproduce the general characteristics of the analyzed variables, not exceeding the minimum and maximum values of the original data and with a relative error in the average not exceeding 5%.
Variable | Original data (%) | Upscaled realization No. 4 (%) |
Upscaled realization No. 24 (%) |
||||||
---|---|---|---|---|---|---|---|---|---|
Min. | Max. | Mean | Min. | Max. | Mean | Min. | Max. | Mean | |
REC14 | 3.21 | 97.52 | 67.63 | 4.06 | 96.94 | 67.55 | 4.23 | 97.23 | 66.97 |
REC400 | 0.55 | 44.86 | 9.36 | 1.24 | 50.34 | 9.57 | 0.79 | 62.14 | 9.67 |
REC-400 | 1.04 | 94.14 | 23 | 1.6 | 85.08 | 22.87 | 1.05 | 92.82 | 23.36 |
Figure 3 shows the non-ergodic correlogram model of the original data (red) and those for the 30 realizations (green) corresponding to the variables Rec14, Rec400, and Rec-400. Note that for all variables, the model satisfactorily reproduced the ergodic fluctuations. For modeling spatial continuity, a non-ergodic correlogram was used (^{Srivastava, 1987}). Table 4 shows correlation matrices between the original variables (a) and the E-type simulations (b). Note that the correlation of the E-type simulations between the variables Rec14 and Rec400 showed a small increase of 0.1. This small increase is a characteristic resulting from the smoothing effect generated by the E-type model (as in kriging).
(a) Original data | |||
Variable | Rec14 | Rec400 | Rec-400 |
Rec14 | 1 | −0.67 | −0.95 |
Rec400 | −0.67 | 1 | 0.42 |
(b) E-type simulations | |||
Variable | Rec14 | Rec400 | Rec-400 |
Rec14 | 1 | −0.77 | −0.95 |
Rec400 | −0.77 | 1 | 0.54 |
Rec-400 | −0.95 | 0.54 | 1 |
A final validation was carried out through checking the local average reproduction (swath plot) by comparing the block grade means versus the declustered data means for each variable respectively. The plots check the E-type model derived from 30 simulations. Figure 4 shows the local averages of the variable Rec14 along the East-West, North-South, and vertical (Z) directions. Note that the model and data mean show good adherence along all directions.
Finally, the closure of the estimated masses retained on each sieve was analyzed. This closure is given by the sum of the percentages retained at the three granulometric fractions at each simulated point. Figure 5 shows the histograms for the closure values at each block for three randomly selected simulations (Nos. 8, 15, and 21). Note that in each case, the closure was at 100%; that is, the sum of the percentages of the total mass at each simulated node was guaranteed to be constant.
4. Conclusion
Simulations of the transformation (ilr) have shown it to be an alternative tool for dealing with compositional data on multi-element mineral deposits for several reasons. It avoids bias in direct kriging of blocks of nonlinearly transformed data by retaining the E-type of multiple simulations. The simulations are satisfactorily validated by the data statistics: they satisfactorily reproduce the basic statistics of the original data, the model of spatial continuity, and the correlations between variables, and they exhibit good adherence between E-type block values and the local data average. All the simulations ensured the granulometric closure of the masses retained on each sieve (100%) at each grid node or block (after upscaling).