Geostatistical clustering as an aid for ore body domaining: Case study at the Rocklea Dome channel iron ore deposit,Western Australia

Abstract

An important step in mineral resource estimation process is the grouping of drill hole samples into domains that reflect zones of homogeneous properties for accurate grade estimation and practical exploitation purposes. In practice, this challenging task is performed through a subjective, time-consuming manual interpretation of the mineral deposit. Therefore, various interpretations are possible. The definition of domains can be viewed as a clustering problem consisting of grouping samples into clusters, herein called domains, so that samples belonging to the same cluster are more similar than those in different clusters. Several methods exist for this purpose; however, groups of samples created through traditional clustering tend to show poor spatial contiguity. Alternatively, spatially contiguous clusters can be obtained through geostatistical clustering where the spatial dependency between samples is considered. This paper is devoted to the application of geostatistical clustering to support domaining of an iron ore deposit located in Western Australia.

Keywords

Clustering domaining drill hole data geostatistics iron ore deposit mineral resource estimation

Introduction

A crucial step in mineral resource evaluation process is ore body domaining for robust grade estimation and practical exploitation purposes (Sousa 1989; Coombes 1997; Abzalov and Humphreys 2002; Emery and Ortiz 2004, 2005; Romary et al. 2012; Rossi and Deutsch 2014; Glacken and Trueman 2014; Duke and Hanna 2014; Abzalov 2016). Ore body domaining consists of separating an ore body into domains of similar characteristics so that blocks inside any of the zones are estimated from measured values of samples within the same zone. Ore bodies are rarely homogeneous because their internal structure commonly exhibits different degrees and types of zoning and layering (Abzalov 2016). Thus, the variables of interest (i.e. mineral grades) may have different spatial behaviour across the region under study in terms of mean, variance and spatial continuity. A domain in this context represents an area or volume within which the variables of interest can be assumed homogeneously distributed (Rossi and Deutsch 2014; Glacken and Trueman 2014; Duke and Hanna 2014). As pointed out by Rossi and Deutsch (2014), this does not mean that the variables of interest are constant within the domains; however, the variables of interest within a domain can be interpreted as a realisation of a stationary random function, simplifying the subsequent geostatistical modelling and estimation steps. In other words, a domain is equivalent to a geostatistical stationary zone.

Domaining is a first-order decision and will influence all subsequent steps in estimation. Thus, poor definition of domains can result in a deterioration of the accuracy of the resource estimates (Abzalov and Humphreys 2002; Emery and Ortiz 2004, 2005). Usually, domaining is performed through a subjective, time-consuming manual interpretation of the mineral deposit. Various interpretations are therefore possible. It is possible to generate domains using the geology logs to divide the region under study into geological units. Unfortunately, there are a number of problems with this approach (Haest et al. 2012b; Cáceres and Emery 2013; Adeli and Emery 2017). First, geological logging provides subjective interpretation of a given sample's composition (e.g. in terms of mineralogy, geochemistry), due to the visual nature of logging. It is commonly accepted that geological logging may not be consistent for a set of drill holes, particularly if they were logged by more than one person. Second, the geological units may not represent the best division of the space into geostatistical stationary zones. Third, this approach is qualitative and time-consuming. Instead of using geology logs, quantitative mineralogy generated from hyperspectral reflectance measurements can be used. It provides a more detailed and accurate characterisation of the ores (Haest et al. 2012b).

The definition of domains can be viewed as a clustering problem in which samples are grouped into clusters herein called domains, where samples belonging to the same cluster are more similar than those in different clusters according to a set of variables of interest. Various methods exist for this purpose; however, classical clustering (which does not include spatial information) tends to produce spatially scattered groups of samples, which is certainly undesirable from the standpoint of mining exploitation. Alternatively, spatially contiguous clusters can be obtained through geostatistical clustering where the spatial dependency between samples is accounted for. Most importantly, automation speeds up the domaining process, facilitating the use of large data sets. Automation also allows rapid updating of the model to include new data or test different model assumptions.

Although geostatistical clustering has been known for a very long time, very few references to this subject can be found in the literature. One way to perform clustering of multivariate spatially continuous indexed data consists of using existing classical clustering methods (e.g. hierarchical clustering, spectral clustering) by modifying the dissimilarity or similarity measure between sample locations to explicitly accounting for the spatial dependency (Olivier and Webster 1989; Bourgault et al. 1992; Fouedjio 2016a,b,c, 2017a,b). This type of approach relies on the joint spatial continuity structure of the variables of interest for providing spatially contiguous clusters. An alternative approach consists of using existing traditional clustering methods according to a set of spatial contiguity constraints (Romary et al. 2012, 2015). These constraints are given by a graph structuring sample locations in the study domain, such as the Delaunay triangulation. The dissimilarities between linked sample locations are computed from attribute values including geographical coordinates as in classical clustering methods.

In this paper, geostatistical clustering based on the approach developed by Fouedjio (2016a, 2017a) is applied to the Rocklea Dome channel iron ore deposit located in Western Australia. The aim is to explore the capability of a geostatistical clustering approach to assist ore body domaining. The available data consist of exploration drill hole samples which are clustered on the basis of geochemical variables (mineral grades) and mineralogical variables (hyperspectral analysis variables).

The remainder of the paper is structured as follows. The geostatistical clustering approach applied in this case study is outlined in the section ‘Geostatistical clustering’. Its application to the Rocklea Dome channel iron ore deposit is presented in the section ‘Application to Rocklea Dome channel iron ore deposit’. Comparisons are made with a traditional clustering technique. Concluding remarks are outlined in the section ‘Conclusion’. Throughout this paper, the terms ‘cluster’ and ‘domain’ will refer to the same thing.

Geostatistical clustering

Consider a set of p variables of interest defined on the region under study typically equal to 2 or 3. All these variables are assumed to be measured at a set of n distinct sample locations . In other words, data are available for each variable at all sampling locations (i.e. isotopic data). A variable of interest could be continuous (e.g. mineral grade) or binary (e.g. presence/absence of a factor). In the case of categorical variables, a categorical variable with c levels will be transformed into c−1 binary variables, each corresponding to one of the c levels of the categorical information. In the case of partial heterotopic data where some variables share only some sampling locations, practitioners must decide between excluding incomplete observations and imputing the missing values (Barnett and Deutsch 2015).

The aim of geostatistical clustering is to group samples into spatially contiguous clusters so that samples belonging to the same cluster have a certain degree of homogeneity while samples in different clusters are as different as possible according to the set of variables of interest. Before the clustering, a processing step is required. Continuous variables of interest need to be standardised. It may also be useful to transform continuous variables of interest to normal scores if their distributions are skewed. This preliminary processing makes the variables comparable and spatial pattern easier to recognise. Geographical coordinates need to be scaled to the same absolute range so they have equal impact on the clustering; especially in 3D where the range of the third coordinate is often relatively small compared to those of first and second coordinates.

Geostatistical (dis) similarity function

The key ingredient in clustering is the (dis) similarity measure which measures the degree of proximity between individual observations. In traditional clustering, the (dis) similarity measure is defined as a function of the Euclidean distance in the space of variables of interest. However, in the geostatistical setting, this type of (dis) similarity measure cannot reflect the spatial continuity structure of data, even if geographical coordinates are also considered as variables of interest (Romary et al. 2015; Fouedjio 2016c).

The first step consists of building a similarity measure that takes into account the spatial dependency between samples. The joint spatial continuity structure of data is commonly described using direct and cross variograms (Wackernagel 2003) defined at any pair of target locations . The joint spatial continuity structure can be captured using the following non-parametric kernel estimator (Fouedjio 2016a,b,c, 2017a,b): (1) where is a non-negative kernel function with constant bandwidth parameter ; takes the value 1 for and 0 for .

The non-parametric kernel estimator of the joint spatial continuity structure defined in Equation (1) is a weighted average of the products of the increments of two variables. It is defined at any pair of locations and including any pair of sample locations. The role of the kernel function in Equation (1) is to weight data locations according to target locations so that data locations close to target locations receive more weight than remote data locations (Figure 1). Regarding the choice of the kernel function , the form of the kernel function is less important than its bandwidth parameter. We use the Epanechnikov kernel shown in Figure 1(b) and whose support is compact and which has optimal properties in density estimation (Wand and Jones 1995). In addition, the use of a kernel function with compact support considerably reduces the computational burden, since the number of terms to calculate decreases in Equation (1). The bandwidth parameter λ is chosen so that the support of the kernel function centred at each sample location contains at least 35 sample locations (Fouedjio 2016b, 2017a). Thus, for each sample location its distance to the 35th closest neighbour is computed; then, the maximum resulting distance is taken as the value of the bandwidth parameter λ. Simulation studies performed by Fouedjio (2016b, 2017a) revealed that this choice is reasonable. The rationale behind this choice is to have a sufficient minimum number of neighbouring samples to estimate the spatial continuity structure reliably.

Given the set of estimated direct and cross variograms defined in Equation (1), one can derive a similarity measure that takes into account the spatial dependency between samples. The similarity measure at two sample locations and is given as follows (Fouedjio 2016b, 2017a): (2) a normalising factor. The term represents the normalised dissimilarity at sample locations and . The similarity measure has the following properties: (i) , (ii) , (iii) and (iv) .

Geostatistical spectral clustering

Given the geostatistical similarity measure specified in Equation (2), a similarity matrix at all pairs of sample locations can be computed and then used as an input to a classical clustering algorithm like spectral clustering. This latter belongs to a class of partitional (non-hierarchical) clustering algorithms that relies on the eigendecomposition (spectral decomposition) of similarity matrices to partition samples into clusters (von Luxburg 2007). Advantages of using spectral clustering include its flexibility in terms of incorporating diverse types of similarity measures, the superiority of its clustering solution compared to traditional clustering methods such as k-means and its well established and appealing theoretical properties.

The spectral clustering algorithm starts by computing the Laplacian matrix from the similarity matrix at all pairs of sample locations. For a given number of clusters k, the spectral clustering algorithm finds the top k eigenvectors of the Laplacian matrix. This operation consists in performing the spectral decomposition (or eigendecomposition) of the Laplacian matrix, hence the term ‘spectral’. The spectral decomposition is the factorisation of a matrix into a canonical form, whereby the matrix is depicted in terms of its eigenvalues and eigenvectors. The top k eigenvectors of the Laplacian matrix define a k-dimensional projection of the data. Then, a standard clustering algorithm such as k-means can be applied in that projected space to derive the final clusters of samples. Specifically, the geostatistical spectral clustering performs the following steps (Fouedjio 2016b, 2017a):

Form the similarity matrix at all pairs of sample locations using Equation (2);

Compute the Laplacian matrix , where is a diagonal matrix;

Compute the k largest eigenvalues of the matrix and form the matrix whose columns are the associated k first eigenvectors of ;

Normalise the rows of the matrix to norm 1;

Cluster the rows of the matrix with the classical k-means clustering algorithm into clusters ;

Assign sample location to the same cluster the tth row of the matrix has been assigned.

The optimal number of clusters is usually determined by evaluating clustering with different numbers of clusters and selecting the one that provided the best performance. Various internal cluster validity indexes can be found in the literature (Hui and Zhongmou 2013). Fouedjio (2016c, 2017a) suggested the use of Caliński–Harabasz index (Caliński and Harabasz 1974) for the geostatistical spectral clustering. This index relies on the between-cluster variation and within-cluster variation. Ideally, we would like to find a solution where between-cluster variation is high and within-cluster variation is low. Thus, the optimal number of clusters is the one that maximises this index defined as follows: (3) where is the overall between-cluster variance and is the overall within-cluster variance; is the vector corresponding to the tth row of the matrix whose columns are the corresponding top k eigenvectors; is the average of points in cluster ; and is the overall average; is the number of points in cluster .

The determination of the optimal number of clusters does not require to repeat all steps of the geostatistical spectral clustering for each value of the predefined number of clusters k. The only step which is carried out for each value of k is the k-means clustering. After creating clusters, it is important to know the contribution that each variable made to the resultant clusters. This information can greatly improve the interpretation of the clusters. It is possible to calculate the importance of a variable in the definition of the clusters using the ratio of the between-cluster variance to the within-cluster variance of this variable. The importance of a variable is then given by: (4) where is the mean of the variable in the cluster ; is the mean of the variable over all samples; is the measurement of the variable at sample location . The more a variable is important, the more its between-cluster variance is high and its within-cluster variance is low. Thus, variables with larger ratios have greater importance.

Application to Rocklea Dome channel iron ore deposit

The geostatistical clustering approach described in the section ‘Geostatistical clustering’ is applied to the Rocklea Dome channel iron ore deposit located in Western Australia. Comparisons are made with traditional clustering which serves as baseline.

Deposit description

The Rocklea Dome channel iron deposit is located in the Hamersley Province of Western Australia (Figure 2(a)). The Tertiary channel iron was deposited atop Archaean to Proterozoic metasedimentary and metavolcanic rocks. Along with the iron oxy-hydroxide minerals, major mineral groups include kaolinite, smectite, carbonate and quartz. Haest et al. (2012b) provide a detailed overview of the geology of the Rocklea Dome and the formation of the channel iron deposit, which are briefly summarised here (Figure 2(b)). The basement geology of the Rocklea Dome area consists of an Archaean monzogranite pluton and cross-cutting mafic and ultramafic intrusives that form part of the Pilbara Craton. The pluton is overlain by Archaean to Proterozoic metasedimentary and volcanic rocks of the Hamersley Province, which envelope the central dome of monzogranite. A meandering tertiary palaeochannel overlies the Archaean and Proterozoic basement and locally contains channel iron ore deposits. Channel iron ore was drilled along 8 km strike length of a palaeochannel on the eastern side of the Rocklea Dome, which was referred to by Haest et al. (2012,b) as the Rocklea Dome channel iron deposit. The deposit is approximately 6400 m long and 1000 m wide.

Data set description

The data available for this study includes exploration drill hole samples with measurements of geochemical variables (mineral grades) and hyperspectral variables (mineral abundance and composition). Holes have been drilled a depth of up to 60 m, with core samples collected at 1.0 m intervals. The geochemical variables of interest comprise IRON (Fe), ALUMINA (Al2O3), SILICA (SiO2), POTASSIUM (K2O), CALCIUM (CaO), MAGNESIUM (MgO), TITANIUM (TiO2), PHOSPHORUS (P), SULFUR (S), MANGANESE (Mn) and LOSS ON IGNITION (LOI). The hyperspectral variables of interest comprise ferric oxide abundance (Fe_ox_ab), hematite–goethite ratio (Hem_Goe) and kaolinite abundance (Kao_ab). The geochemical data were collected by Kalassay Ltd, using a Bruker Pioneer X-ray fluorescence instrument, equipped with an end window 4 kW rhodium X-Ray tube. Hyperspectral analysis was conducted using CSIRO's HyChips system, which collects sample spectra from a mm area. For a detailed description of the geochemical and hyperspectral analyses, the reader is referred to Haest et al. (2012). Mineral abundance (i.e. ferric oxide abundance, kaolinite abundance) and compositional (i.e. hematite–goethite ratio) variables were computed in The Spectral Geologist software (TSG ) using a multiple feature extraction method (Laukamp et al. 2010).

After the preprocessing step, 2637 out of 7520 core samples have been removed. They correspond to samples having missing entries for all geochemical measurements. These core samples are related to the upper, poorly mineralised, ∼10 m of drill holes where geochemical measurements have not been taken. The data used for the domaining correspond to a set of 4883 samples from 156 drill holes. The map of drill hole sample locations of interest is shown in Figure 3. Figure 4 presents the spatial plot of main geochemical variables (Fe, Al2O3, SiO2) and mineralogical variables (Fe_ox_ab, Hem_Goe, Kao_ab). The univariate distributions of main geochemical variables and mineralogical variables are given in Figure 5. One can observe that the histogram of the main economic grade component (Fe) has a bimodal shape clearly indicating for the presence of the different groups of data. Table 1 reports some descriptive statistics of geochemical and mineralogical variables. Most of the variables have skewed distributions. Thus, before the clustering, all variables will be transformed to normal scores. This transformation also allows the variables to be directly compared and clarifies spatial patterns in the data. The relative order is maintained so high normal scores correspond to high raw values and vice versa.

Table 1

General statistics of geochemical and mineralogical variables.

Variable	Unit	Mean	Stdv	Min	Max	0.25Q	Median	0.75Q	Skewness
Fe	%	31.770	17.764	0.000	62.000	14.000	35.000	48.500	−0.222
Al2O3	%	10.890	7.376	0.320	39.430	5.350	8.820	15.130	0.933
SiO2	%	28.140	19.738	2.380	92.810	11.140	22.460	42.220	0.805
K2O	%	0.180	0.492	0.000	9.504	0.022	0.048	0.136	7.691
CaO	%	1.498	4.192	0.010	40.820	0.100	0.190	0.500	4.448
MgO	%	1.246	2.244	0.000	20.010	0.280	0.470	1.020	3.886
TiO2	%	0.627	0.517	0.003	3.742	0.227	0.482	0.890	1.321
P	%	0.025	0.038	0.000	0.845	0.006	0.019	0.035	10.070
S	%	0.051	0.362	0.000	14.145	0.005	0.011	0.023	22.46
Mn	%	0.082	0.151	0.000	6.160	0.000	0.040	0.120	17.409
LOI	%	11.280	4.624	0.000	45.280	9.690	11.200	11.900	2.513
Fe_ox_ab		0.208	0.070	0.000	0.355	0.163	0.225	0.259	−0.676
Hem_Goe		911.300	7.179	871.900	974.100	908.800	911.300	914.000	0.647
Kao_ab		0.129	0.075	0.017	0.412	0.068	0.112	0.174	0.840

The Spearman correlation coefficients between the variables are shown in Figure 6. Note that the Spearman correlation coefficient, which measures the correlation of ranks, is invariant under a monotonic transformation of the data such as normal scores transformation. Many of the variables are weakly correlated, but there are some strong correlations. Pairs of variables with strong positive correlations are: IRON (Fe)–ferric oxide abundance (Fe_ox_ab), ALUMINA (Al2O3)–TITANIUM (TiO2) and CALCIUM (CaO)–MAGNESIUM (MgO). The variable pair: IRON (Fe)–SILICA (SiO2) shows a strong negative correlation. The strong positive correlation between the pXRF-derived IRON (Fe) content and the hyperspectrally derived ferric oxide abundance (Fe_ox_ab) confirms the observations by Haest et al. (2012,b). The strong positive correlation between the geochemically derived ALUMINA (Al2O3) and TITANIUM (TiO2) contents may be related to heavy minerals associated with detrital material in clay-rich layers of the channel. Haest et al. (2012) described extensive calcrete in the northern part of the Rocklea Dome channel iron deposit, which may explain the high correlation between CALCIUM (CaO) and MAGNESIUM (MgO). The strong negative correlation between IRON (Fe) and SILICA (SiO2) is likely attributable to the intercalation of Al-clay-rich sedimentary layers (high in silica, low in iron) with iron-oxide-rich horizons which are relatively low in silica. This is supported by the positive correlation between SILICA (SiO2), ALUMINA (Al2O3) and POTASSIUM (P).

Results and discussions

Figures 7, 8 and 9 show respectively a three-dimensional view, a plan view and a vertical view of domains produced by classical k-means clustering and geostatistical spectral clustering for a different number of domains ( ). These plots illustrate that domains created through classical clustering turn out to show poor spatial contiguity compared to ones obtained by geostatistical clustering. Moreover, as the number of domains increases, the spatial contiguity decreases. Geostatistical clustering provides domains with much clearer spatial contiguity. It produces homogeneous regions which are spatially contiguous both along the length of the drill holes (depth direction) and between drill holes (Northing and Easting directions).

Figure 10(a) plots the number of domains versus the Caliński–Harabasz index defined in Equation (3). In addition to the Caliński–Harabasz index, two other well-known internal cluster validity indexes are plotted (Figure 10(b, c)), namely, silhouette index (Rousseeuw 1987; Kaufman and Rousseeuw 1990) and Davies–Bouldin index (Davies and Bouldin 1979). The silhouette index relies on the pairwise difference of between-cluster distances and within-cluster distances. In turn, the Davies–Bouldin index is based on a ratio of within-cluster and between-cluster distances. The maximum in the plot of the Caliński–Harabasz index (or the silhouette index) versus the number of domains is taken to indicate the underlying number of domains. The minimum in the plot of the Davies–Bouldin index versus the number of domains is an indication of the relevant number of domains. It appears that the relevant number of domains according to the Caliński–Harabasz index as well as the two other indexes is two. It is important to keep in mind that during the data processing step, a part of the set of core sample data have been removed. This part is related to the upper, poorly mineralised, ∼10 m of drill holes where geochemical measurements have not been taken.

As highlighted by Guibal (2014), the statistics and variography should back up the definition of domains. Table 2 reports some descriptive statistics of geochemical and mineralogical variables corresponding to the two optimal domains. It is clear that the statistical contrast between the two domains is substantial. Domain 1 (sky blue points in Figure 7(b)) is characterised by high IRON (Fe), PHOSPHORUS (P) and MANGANESE (Mn) grades and by low ALUMINA (Al2O3), SILICA (SiO2), POTASSIUM (K2O), CALCIUM (CaO), MAGNESIUM (MgO), TITANIUM (TiO2) and SULFUR (S) grades. By contrast, Domain 2 (red points in Figure 7(b)) is characterised by low IRON (Fe), PHOSPHORUS (P) and MANGANESE (Mn) grades and by high ALUMINA (Al2O3), SILICA (SiO2), POTASSIUM (K2O), CALCIUM (CaO), MAGNESIUM (MgO), TITANIUM (TiO2) and SULFUR (S) grades. For some geochemical components (Al2O3, SiO2, K2O, CaO, MgO, TiO2 and S), the standard deviation of the samples within domain is likely proportional to their mean grade; while for other mineral grade elements (Fe, P and Mn), the standard deviation of the samples within domain is likely inversely proportional to their mean grade. Domain 1 contains 2305 drill hole samples while the domain 2 contains 2578 drill hole samples. Domain 1 is located mainly at the middle of the deposit, i.e. the central zone, containing most of the main channel iron ore resource. Domain 2 is located at the extremities of the deposit (zones north-west and south).

Table 2

Geostatistical spectral clustering: descriptive statistics of the variables of interest corresponding to the two optimal domains.

		Domain 1 ( )				Domain 2 ( )
Variable	Unit	Mean	Stdv	CV	Median	Mean	Stdv	CV	Median
Fe	%	44.617	10.913	4.089	48.000	20.278	14.545	1.394	16.000
Al2O3	%	7.322	5.015	1.460	6.230	14.080	7.682	1.833	12.990
SiO2	%	15.690	11.151	1.407	12.110	39.275	19.085	2.058	38.950
K2O	%	0.042	0.149	0.282	0.023	0.303	0.638	0.475	0.113
CaO	%	0.420	1.298	0.323	0.150	2.463	5.461	0.451	0.240
MgO	%	0.536	0.720	0.745	0.330	1.881	2.868	0.656	0.725
TiO2	%	0.399	0.353	1.130	0.294	0.832	0.553	1.503	0.719
P	%	0.032	0.020	1.642	0.031	0.019	0.048	0.399	0.008
S	%	0.029	0.240	0.122	0.012	0.070	0.443	0.158	0.010
Mn	%	0.112	0.097	1.159	0.100	0.055	0.182	0.304	0.010
LOI	%	11.265	2.071	5.439	11.460	11.289	6.056	1.864	10.590
Fe_ox_ab		0.254	0.043	5.945	0.256	0.167	0.065	2.578	0.174
Hem_Goe		910.669	3.709	245.537	910.660	911.817	9.205	99.060	912.450
Kao_ab		0.092	0.054	1.716	0.078	0.161	0.077	2.110	0.150

Figure 11 shows the importance of each variable in the formation of two optimal domains. It appears that the six most important variables are IRON (Fe), ferric oxide abundance (Fe_ox_ab), SILICA (SiO2), kaolonite abundance (Kao_ab), ALUMINA (Al2O3) and TITANIUM (TiO2). This result stresses the importance of the mineralogical information for domaining. Figure 12 presents histograms of the six most important variables. One can see that distributions of these variables differ markedly between the two domains. Moreover, the histogram of the main economic grade component (Fe) is now unimodal in each domain. Domain variography of the six most important variables is given in Figures 13 (vertical variograms) and 14 (horizontal variograms). One can observe that the spatial continuity structure of these variables differs from one domain to another. In both domains, these variables show a high continuity along the length of drill hole (vertically) while exhibiting less continuity between drill holes (horizontally).

The geology of the Rocklea Dome channel iron deposit can be separated into four geological units (Haest et al. 2012b): (1) the bedrock, represented by basalts of the Archaean Fortescue Group, (2) the iron ore, represented by the Tertiary channel iron deposit (CID) and overlying reworked iron-oxide-rich material (vitreous and ochreous goethite), (3) clay-rich (mainly kaolin) layers surrounding and intercalated within the CID and (4) Quaternary detritals, calcrete and silcrete. Only the first three of these geological units were included in the modelling presented in our paper. In contrast to the remaining three geological units, the geostatistical spectral clustering method suggests that two domains are sufficient to provide a spatially contiguous ore grade model, with the six most important variables being the geochemically derived ALUMINIUM, IRON, SILICA and TITANIUM contents and the hyperspectrally derived ferric oxide and kaolinite abundance. Even when considering the full set of evaluated variables (Figure 6), the bipartite nature of the local geology becomes evident. The iron ore is characterised by a correlation of geochemically derived IRON, PHOSPHORUS and MANGANESE as well as high ferric oxide abundance. Compared to the iron ore, the clay-rich interlayers, as well as the bedrock, are high in ALUMINIUM, SILICA, POTASSIUM, TITANIUM and kaolinite abundance. The hyperspectrally derived parameter of the hematite–goethite ratio, which was used by Haest et al. (2012b) for iron-oxide speciation, is only applicable in samples that have a high ferric oxide abundance. In summary, for the purpose of ore grade modelling, two domains are sufficient. For geometallurgical modelling, the iron ore domain could be divided into more subdomains, such as separating the vitreous goethite-rich layers from the ochreous goethite-rich ones (Haest et al. 2012b).

Conclusion

The application of geostatistical clustering for assisting domaining of the Rocklea Dome Channel Iron Ore Deposit revealed that it is distributed in two distinct domains. These domains are mainly characterised by four geochemical variables (IRON, ALUMINA, SILICA and TITANIUM) and two mineralogical variables (ferric oxide abundance and kaolinite abundance). These variables have distinctive statistics and spatial structures in the two domains. From this case study, the advantages of using geostatistical clustering in support of ore body domaining can be acknowledged. Ore body domaining can be achieved in a less subjective, simple and effective way using geostatistical clustering. By taking into account the spatial continuity structure of the data, geostatistical clustering can produce homogeneous groups of samples (domains) which are spatially contiguous along and across the drill holes. The optimal number of domains can be determined automatically through some cluster validity indexes. This reduces the risk of creating more domains than required (over-domaining) which can compromise the integrity of the grade distribution. The importance of variables to the domaining can also be provided. This information can greatly improve the characterisation of the domains and can also be used as a guide by the mining company to optimal data collection strategies. Geostatistical clustering is easily reproducible, making the incorporation of new data relatively straightforward. Although the iron ore deposit presented in this paper is not complex, geostatistical clustering as a data-driven approach can effectively help to domain complex mineral deposits if there are enough data to support this complexity.

Footnotes

Disclosure statement

No potential conflict of interest was reported by the authors.

References

Abzalov

2016. Applied mining geology. Vol. 12. Cham: Springer

Sage UK: London, England

Abzalov

, Humphreys

2002. Geostatistically assisted domaining of structurally complex mineralisation: method and case studies. The AusIMM 2002 Conference, Auckland, New Zealand. Vol. 150. p. 345–350

Adeli

, Emery

2017. A geostatistical approach to measure the consistency between geological logs and quantitative covariates. Ore Geol Rev. 82C:160–169. doi: 10.1016/j.oregeorev.2016.11.028

Barnett

, Deutsch

CV.

2015. Multivariate imputation of unequally sampled geological variables. Math Geosci. 47(7):791–817. doi: 10.1007/s11004-014-9580-8

Bourgault

, Marcotte

, Legendre

1992. The multivariate (co)variogram as a spatial weighting function in classification methods. Math Geol. 24(5):463–478. doi: 10.1007/BF00890530

Cáceres

, Emery

2013. Geostatistical validation of geological logging. In: Ambrus J, et al., editors. Proceedings of the Third International Seminar on Geology for the Mining Industry. Santiago: Gecamin Ltd. p. 73–80.

Caliński

, Harabasz

1974. A dendrite method for cluster analysis. Commun Statist. 3(1):1–27.

Coombes

1997. Handy hints for variography. Towards 2000 – National Conference on Ironmaking Resources and Reserves Estimation, Perth, WA. p. 127–130

Davies

, Bouldin

DW.

1979. A cluster separation measure. IEEE Trans Pattern Anal Mach Intell. PAMI-1(2):224–227. doi: 10.1109/TPAMI.1979.4766909

10.

Duke

, Hanna

2014. Geological interpretation for resource modelling and estimation. In: Edwards AC, editor. Mineral resource and ore reserve estimation – the AusIMM guide to good practice. Melbourne: The Australasian Institute of Mining and Metallurgy. p. 145–153.

11.

Emery

, Ortiz

2004. Defining geological units by grade domaining. Tech. rep., Universidad de Chile.

12.

Emery

, Ortiz

2005. Estimation of mineral resources using grade domains: critical analysis and a suggested methodology. J South Afr Inst Min Metall. 105:247–256.

13.

Fouedjio

2016a. A clustering approach for discovering intrinsic clusters in multivariate geostatistical data. In: Perner P, editor. MLDM 2016, LNCS. Vol. 9729. Cham: Springer. p. 491–500

14.

Fouedjio

2016b. Discovering spatially contiguous clusters in multivariate geostatistical data through spectral clustering. In: Li J, et al., editors. ADMA 2016, LNCS. Vol. 10086. Cham: Springer. p. 547–557

15.

Fouedjio

2016c. A hierarchical clustering method for multivariate geostatistical data. Spat Stat. 18:333–351. doi: 10.1016/j.spasta.2016.07.003

16.

Fouedjio

2017a. A spectral clustering approach for multivariate geostatistical data. Int J Data Sci Anal. 4(4):301–312. doi: 10.1007/s41060-017-0069-7

17.

Fouedjio

2017b. A spectral clustering method for large-scale geostatistical datasets. In: Perner P, editor. MLDM 2017, LNCS. Vol. 10358. Cham: Springer. p. 248–261

18.

Glacken

, Trueman

2014. Overview – mineral resource estimation. In: Mineral resource and ore reserve estimation - the AusIMM guide to good practice. Edwards AC, editor. Melbourne: The Australasian Institute of Mining and Metallurgy. p. 263–276.

19.

Guibal

2014. Variography, a tool for the resource geologist. In: Mineral resource and ore reserve estimation – the AusIMM guide to good practice. Edwards AC, editor. Melbourne: The Australasian Institute of Mining and Metallurgy. p. 85–90.

20.

Haest

, Cudahy

, Laukamp

, Gregory

2012a. Quantitative mineralogy from infrared spectroscopic data – II. Three-dimensional mineralogical characterization of the Rocklea channel iron deposit, Western Australia. Econ Geol. 107(2):229–249. doi: 10.2113/econgeo.107.2.229

21.

Haest

, Cudahy

, Laukamp

, Gregory

2012b. Quantitative mineralogy from visible to shortwave infrared spectroscopic data – I. validation of mineral abundance and composition products of the Rocklea Dome channel iron deposit in Western Australia. Econ Geol. 107(2):209–228. doi: 10.2113/econgeo.107.2.209

22.

Hui

, Zhongmou

2013. Clustering validation measures. Chapman and Hall/CRC data mining and knowledge discovery series. p. 571–605.

23.

Kaufman

, Rousseeuw

1990. Finding groups in data: an introduction to cluster analysis. New York: Wiley.

24.

Laukamp

, Caccetta

, Chia

, Cudahy

, Gessner

, Haest

, Liu

, Ong

, Rodger

2010. The uses, abuses and opportunities for hyperspectral technologies and derived geoscience information. AIG Bull. 51:73–76.

25.

Olivier

, Webster

1989. A geostatistical basis for spatial weighting in multivariate classification. Math Geol. 21:15–35. doi: 10.1007/BF00897238

26.

Romary

, Ors

, Rivoirard

, Deraisme

2015. Unsupervised classification of multivariate geostatistical data: two algorithms. Comput Geosci. 85, Part B:96–103. doi: 10.1016/j.cageo.2015.05.019

27.

Romary

, Rivoirard

, Deraisme

, Quinones

, Freulon

2012. Domaining by clustering multivariate geostatistical data. In Abrahamsen P, et al., editors. Geostatistics Oslo 2012, Quantitative geology and geostatistics. Vol. 17. Dordrecht: Springer. p. 455–466

28.

Rossi

, Deutsch

CV.

2014. Mineral resource estimation. Dordrecht: Springer

Sage UK: London, England

29.

Rousseeuw

1987. Silhouettes: a graphical aid to the interpretation and validation of cluster analysis. J Comput Appl Math. 20:53–65. doi: 10.1016/0377-0427(87)90125-7

30.

Sousa

AJ.

1989. Geostatistical data analysis – an application to ore deposit. In: Geostatistics, Vol. 2, Proceedings of the Third International Geostatistics Congress, September 5–9, 1988. Avignon: Kluwer Academic Publishers. p. 851–860.

31.

von Luxburg

2007. A tutorial on spectral clustering. Stat Comput. 17(4):395–416. doi: 10.1007/s11222-007-9033-z

32.

Wackernagel

2003. Multivariate geostatistics: an introduction with applications. Berlin, Heidelberg: Springer.

33.

Wand

, Jones

1995. Kernel smoothing. Monographs on statistics and applied probability. London: Chapman and Hall

Sage UK: London, England