Machine-Learning for Mineral Identification and Ore Estimation From Hyperspectral Imagery in Tin-Tungsten Deposits

This study aims to assess the feasibility of delineating and identifying mineral ores from hyperspectral images of tin-tungsten mine excavation faces using machine-learning classification. We compiled a set of hand samples of minerals of interest from a tin-tungsten mine and analyzed two types of hyperspectral images: 1) images acquired with a laboratory set-up under close-to-optimal conditions; and 2) scan of a simulated mine face using a field set-up, under conditions closer to those in the gallery. We have analyzed the following minerals: cassiterite (tin ore), wolframite (tungsten ore), chalcopyrite, malachite, muscovite, and quartz. Classification (Linear Discriminant Analysis, Singular Vector Machines and Random Forest) of laboratory spectra had a very high overall accuracy rate (98%), slightly lower if the 450 – 950 nm and 950 – 1780 nm ranges are considered independently, and much lower (74.5%) for simulated conventional RGB imagery. Classification accuracy for the simulation was lower than in the laboratory but still high (85%), likely a consequence of the lower spatial resolution. All three classification methods performed similarly in this case, with Random Forest producing results of slightly higher accuracy. The user’s accuracy for wolframite was 85%, but cassiterite was often confused with wolframite (user’s accuracy: 70%). A lumped ore category achieved 94.9% user’s accuracy. Our study confirms the suitability of hyperspectral imaging to record the spatial distribution of ore mineralization in progressing tungsten-tin mine faces.


Introduction
Building and updating 3D models to quantitatively characterize the distribution and heterogeneity of ore grade, vein thickness, orientation, and network geometry are critical for the mining industry as it provides a better understanding of the mineralization with important benefits for the management and optimization of the exploitation, thus reducing environmental impact. Typically, these 3D models are based on knowledge of the site's structural geology, with geophysical exploration contributing to provide a better understanding of the deposit. Image data acquired from the mine excavation itself could be used to confirm or update these models, monitoring ore grade and distribution. In this regard, we aim at integrating proximal ground-based conventional and hyperspectral imaging of the mine excavation face to provide the ore grade and 2D geometry of the mineralization in the section. Ideally, a 3D tomography could be derived from these 2D planes as the excavation front progresses.
Hyperspectral images produced by imaging spectrometers are 3D arrays in which each voxel holds a radiance spectrum that is processed to reflectance [1] [2]. Acquisition of single reflectance spectra with spectrometers in the visible, near-infrared, and shortwave infrared wavelength domains (400 -2600 nm) is a relatively simple and non-invasive technique that has been used in the laboratory and the field for decades ago [3] [4]. 1999 The particular optical and electronic properties of each material result, under images in the VISNIR (350 -1000 nm) range for the automated identification of metallic ore species, which was achieved with very high accuracy.
Deployment and operation of hyperspectral imaging systems in the field is always involved and doing so under gallery mine conditions is challenging. As a previous step to actually acquire, process and analyze hyperspectral images of a mine excavation face, we explore here the feasibility and interest of such an approach by using hand samples that feature the minerals of interest. For this purpose, we conducted the following studies: 1. Laboratory imaging spectrometry. We scanned the hand samples using hyperspectral cameras on a laboratory set-up to assess spectral separability and evaluate machine-learning classification methods under close-to-optimal conditions. 2.
Simulation of hyperspectral imaging of the mine face. We scanned the whole set of hand samples with the same field set-up, illumination, and distance to object that are expected to be used in the mine gallery, evaluating machine-learning methods to identify and map the distribution of materials in the resulting image.

Geological settings
The studied samples were collected in the San Finx tin-tungsten mine, which is located in Lousame, A Coruña, Galicia, NW Spain. San Finx is a typical case of Late Palaeozoic, granite-related hydrothermal deposit, associated to the European metallogenic belt. Tin-tungsten ore is directly related to granites formed from Late Devonian to Permian during the Variscan or the so-called Hercynian orogeny, resulting from the continental collision of Laurussia and Gondwana [36]. The San Finx mineralization consists of a subvertical set of quartz lodes prevalent in the N 50 o E direction, which at its eastern end shows pegmatitic features (e.g., the occurrence of K-feldspar). This lode field reaches more than 1 km in length, and the lode thickness varies between 0.5 and 1.5 m [37]. Host rocks consist of schists, migmatites and granites, while the metalogenesis consists of Sn-W-Ta-Nb-Mo-Cu-As-Au-Ag-Bi. The exploited ore minerals are cassiterite and wolframite. Other minerals occurring in the quartz lodes are arsenopyrite, pyrite, scheelite, chalcopyrite, molybdenite columbite-tantalite, muscovite and tourmaline, among others [38,39]. The San Finx deposit is declared as a Site of Geological Interest (http://info.igme.es/ielig/LIG-Info.aspx?codigo=GM044).

Laboratory imaging spectrometry
In order to achieve close-to-optimal conditions, we scanned samples from the mine with Specim FX10 and FX17 cameras (Table 1) consecutively mounted on the same Specim LabScanner Setup 40 x 20 and with care not to move the samples between each scan. We used Specim's Lumo software for scanning, which renders hyperspectral images of radiance for the whole scanned area, and two images of, respectively, the dark (internal shutter) and white references, which we used to calculate images of spectral reflectance for each camera.  [48]. The co-registration processing proceeded in two steps. First, we run a firstorder polynomial warp with homologous points automatically extracted by the SIFT method (HomologousPointExtraction in OTB; [49]). Second, the co-registration was refined by calculating subtle local X and Y shifts that optimized local correlation between the corresponding bands of both cameras (FineRegistration in OTB). Geometric transformations were calculated using corresponding bands of the same wavelength in both cameras and then applied to all bands of FX17 (gridBasedImageResampling in OTB).
Hyperspectral images were displayed and sets of pixel spectra for each mineral were interactively extracted with the EnMap Box plugin [50] of QGIS [51], which required combining photo-interpretation and direct visual inspection of the samples. Finally, we imported the spectra into R, where we made comparative graphics along with reference spectra extracted from spectral libraries by USGS [52] and JPL [53] [54] and run statistical analysis. We calculated matrices of dissimilarities using the Jeffries-Matutsita metric [55] to measure the spectral resolving power of each camera (FX10, FX17, and combined FX10_FX17). We also simulated the RGB values of a conventional camera (Canon 60D) from the FX10 values using camera spectral sensitivity curves [56] and calculated the RGB dissimilarity matrix with these values as well to evaluate the power of conventional photos to identify the minerals of interest. Finally, we run three different classification methods: Linear Discriminant Analysis (LDA), Singular Vector Machines (SVM), and Random Forest (RF) for each of the four datasets (FX10, FX17, combined FX10_FX17, and Canon 60D), which basics are briefly introduced in section 3.3. To this end, we used R packages MASS [57], e1071 [58], RandomForest [59], and caret [60].

Ground-based panoramic hyperspectral imaging of simulated mine face
In order to prepare for the deployment of the hyperspectral equipment and analyze hyperspectral imagery acquired under conditions closer to those in the mine, we assembled a panel of 52 cm x 57 cm with the samples, along with a white reference (Sphere Optics SG 3141 95%), a Multi-Component Wavelength Standard reference (Labsphere WCS-MC-020) and a color chart ( Figure 2). The image was acquired in a dark room by a Specim FX17 camera mounted on a Specim Rotary scanner on top of a tripod at a distance of 2.5 m from the lens to the panel. This distance was selected as a realistic choice in the gallery. Illumination was provided by a Fresnel FilmGear 650w 3200k studio light with a diffuser behind and above the camera. All systems were powered by batteries as they would be under gallery mine conditions. Image resolution was 2.6 mm/pixel and covered an area of 166.4 cm (height) x 106.1 cm (width), from which the subscene corresponding to the samples in the panel was extracted. Despite the diffuser, some illumination unevenness could be appreciated, which we solved by applying a Rolling Ball Background Subtraction [61] in Fiji [62] to the first Principal Component and then applying the inverse transformation. As in the case of laboratory images, the resulting hyperspectral image was displayed in QGIS, and sets of pixel spectra for each mineral and other targets (Table S1) were extracted by interactively digitizing training and validation polygons, a task that required direct visual inspection of the samples as well. Finally, we input the file of polygons along with the hyperspectral image to the classifiers (LDA, SVM, and RF), and their respective accuracy assessments. Each classifier produced a predicted map of the panel materials according to the pre-defined set of categories (Table S1). We used tools in EnMap Box [30], which are based on scikit-learning, for SVM and RF classification, while R packages MASS [57] and raster [41] were used for LDA, as in the case of the laboratory spectra. As combining results of different methods produces more robust results, we also combined the classifications by selecting, for each pixel, the class that had been predicted by most methods. For those pixels in which each method was predicting a different class, we selected the class predicted by the classifier with the highest accuracy.

Linear Discriminant Analysis (LDA)
Linear Discriminant Analysis, like Principal Component Analysis (PCA), is often used as a dimensionality reduction technique before machine learning applications. LDA seeks to project the dataset in a space of fewer dimensions with maximum separability among classes, by maximizing the relationship between the within-class variance and the among classes variance. Unlike PCA, which is an unsupervised algorithm (it does not require training polygons for each class as it maximizes total variance for each principal component), LDA is supervised: it transforms the original space into components that linearly maximize the separation among classes, based on estimates of the mean and covariance matrix of each class, which are calculated on a training set of labelled observations [55] [63]. Classification is performed by applying a Bayes rule assuming multi-variate Gaussian distributions of common variance, which simplifies to the nearest centroid rule modulo the prior class probabilities in the transformed space, as the covariance matrix becomes the identity matrix. For the present work, we performed LDA using R with package MASS [57] for both the laboratory spectra and the hyperspectral image of the panel.

Support Vector Machines (SVM)
As comprehensively explained in [63], SVM are a generalization of the Maximal Margin Classifier (MMC). Given a p-dimensional space of descriptors and a simple binary (2 classes) problem with n training observations (the method can be generalized to >2 classes), the MMC seeks an n-1 hyperplane that linearly separates the two classes with the maximum margin (that is: with the maximum distance from the nearest point of each class to the hyperplane). Maximizing the margin increases the chances of having a correct classification of the rest of the observations (those not included in the training set). On some occasions, the distribution of the observations is such that the correct separation of the 2 classes according to MMC implies that the maximal margin hyperplane still lies very close to some observations and the distance between margins is very narrow. In these cases, a more tolerant approach defining a wider margin even at the expense of some errors in the classification of the training set would result in fewer errors with newer observations. This method is named Support Vectors Classification (SVC) by [63] and linearSVC in the scikitlearn library [64]. The degree of tolerance is based on the relative distances from observations to their class margin: observations at the correct side of their margin or on the margin itself are at a distance 0; observations beyond the margin of their respective class but within the limit of the hyperplane are at distances 0< e <1, while those at the wrong side of the hyperplane are at distances >1. The user sets the total tolerance (named "regularization" in the SVM jargon), that is, the allowed sum of all distances from observations to margins. For C=0 there is no tolerance, so hyperplane and margins will be set according to MMC. At increasing values of C, the margins can be progressively widened, with progressively more observations being left at the wrong side of the margin or even at the wrong side of the hyperplane, and actually, the hyperplane itself changes as different observations are included. Parameter C thus controls the number of observations that are actually involved in computing the hyperplane and its margins. It is important to note that unlike other classification approaches such as LDA, SVM classification depends only on those observations lying on the margins or beyond them, and these observations are known as support vectors (as each observation is a vector of p descriptors). An optimal value for C is calculated by k-fold cross-validation in the training set for a predefined range of C values.
There are cases in which no hyperplane can separate the classes, no matter the value of C. In these cases, rather than fitting non-linear functions, the SVM approach is to increase the dimensionality by adding new dimensions that are non-linear transforms of the p original ones. As the parameters defining the hyperplane in SVC (or linear SVC) are found based on the inner product between all pairs of training observations, SVM classification (named SVC in scikit-learn) uses a kernel-based approach for this purpose. Popular kernels in SVM classification are the polynomial and, in particular, the radial basis function (RBF) because of its flexibility. For RBF, we have, for any two training observations x,y (vectors in space of p descriptors): where γ is a positive constant that can be thought of as the inverse of the variance: for low values of γ, the class will be very wide, and for small values, the class will be narrow.
As for the case of C, γ can be set by cross-validation within a grid-search strategy as is the case of EnMap Box. R package caret, instead, uses an analytical formula to get reasonable estimates of γ and fix it to that value [60].

Random Forest
Random forest classification [65] is a development of classification trees. Classification trees are produced by binary and hierarchically splitting the training set by thresholding the descriptors [63]. At each level, all descriptors are sequentially tried and the best threshold for the best predictor is selected to optimize a class purity metric such as the Gini index. Once the tree is grown in the top-down direction, it is then simplified ("pruned" in the classification trees jargon) in the bottom-up direction according to an error minimization rule. One problem with classification trees is that while the fit to the training set can be very good, the prediction of new cases is very dependent upon the specific training set that has been used.
Results are greatly improved by Random Forests, which introduce two modifications: a set of n trees is constructed (with no pruning) instead of just one, and only a subset of m descriptors is used. The training set is resampled with replacement in n bootstrapped sets and a classification tree is produced for each set, but only a random selection of m descriptors is considered at each level. Each tree is then applied to the test data, which yields as many predictions as bootstrapped training sets, and the final result is the set of most commonly selected ("voted") class. A value of m = sqrt(p) is a common recommendation, and the number of trees (n) should be large enough to ensure that all descriptors have been considered.

Validation.
Validation was conducted through the analysis of confusion matrices, which were built by cross-validation: Leave-one-out in the case of LDA, and 10-fold cross validation in the case of SVM and RF. We calculated standard metrics from the confusion matrices using R package caret: overall accuracy (OA), user's accuracy (UA), producer's accuracy (PA), and F1 (which is the harmonic mean of user's and producer's accuracy [66][67][68]. Specifically, considering pixels in the digitized polygons,

Laboratory imaging spectrometry
Reflectance spectra of single-mineral polygons of cassiterite, malachite and muscovite measured in this work are similar to those found in the reference spectral libraries (Figure 3). Those of quartz targets are highly variable (as in the case of spectral libraries) because these spectra are very dependent upon the presence of other minerals or fluid inclusions as well as the effect on the optical properties of quartz of impurities and crystalline disorder. In the case of chalcopyrite, a set of spectra were much brighter than the Spectral separability (Jeffries-Matutsita index) among combined spectra of both cameras (450 -1650 nm) is very high (>1.99) for all mineral pairs. In accordance, LDA classification had very high accuracy (98%, Table 2). Spectral differences are well mapped in the LDA space ( Figure. 4), in which spectral samples are well ordinated in consistency with their respective mineral identity. Note that while spectra of wolframite and cassiterite overlap in the plane of the first two components, they are clearly separated by the 4 th component. Only 2 spectral samples had conflicting labelling between the "de visu" inspection and the LDA classification.  Figure 3. Ordination of laboratory spectra on planes defined by LD components 1 to 4. Colors correspond to classified spectra. Letters indicate the actual mineral of those spectra for which the classification was wrong.
Classification accuracy was very high with data from the hyperspectral cameras, slightly lower when spectra of each camera were taken separately (Table 3): 95% for the FX10 (450 -950 nm) and 94% for the FX17 (950-1650 nm). Instead, classification accuracy was much lower (0.745) if only simulated RGB values were used. Confusion among minerals (Table 3) slightly increases for chalcopyrite, cassiterite and malachite if only the FX10 spectra are considered, and for cassiterite and muscovite if only the FX17 spectra are considered. Confusion significantly increases if only RGB values are considered, in particular for cassiterite, which is always confused with wolframite, but also for chalcopyrite and wolframite. Interestingly, more modern classification methods such as SVM and RF achieved similar or even slightly lower performance than LDA (Table 4).

Ground-based panoramic hyperspectral imaging of simulated mine face
Each classification method applied to the FX17 hyperspectral image of the panel produces a map of the distribution of surface materials in the panel. Classification accuracy is lower than in the case of the laboratory spectra, but still high (Table 5)  The confusion matrices (Table 6 and S2) reveal that while a high accuracy was achieved for wolframite, cassiterite is too often confused with the former (in the case of the RF classification, a 37.5% of the all actual Cassiterite pixels are labelled as wolframite), but more rarely confused with any other material (10.6%). In other words, while the identification of wolframite was reliable, pixels identified as cassiterite can actually be either cassiterite or wolframite. As cassiterite and wolframite are the main ores in this mine, an operationally interesting product can result by lumping together both minerals in an "ore" category that reaches 94.8% of user's accuracy and 93.9% of producer's accuracy. Figure 4 displays the combined result: all three classification methods agree on 53.8% of the study area, while at least 2 methods agree in 87.5% of the study area. Within the training/validation polygons, the results of the three classifiers agree in 90% of the area and, in that case, the agreed class is always correct.

Discussion and Conclusions
Three factors can be called to explain the decrease of accuracy between results obtained in the laboratory and in the simulation of gallery conditions: narrower spectral range (450 -1780 nm in the laboratory while 950 -1780 nm in the gallery simulation), poorer illumination (dedicated illumination system in the laboratory stand vs. standard photography lighting system in the gallery simulation), and coarser pixel size because increased camera -target distance (0.224 mm vs. 2.6 mm). The narrower spectral range is a consequence of the logistic requirement of using one single hyperspectral imaging system in the gallery. While scanning with a compact system covering the 450 nm -1780 nm range (or more) would be better, our laboratory results indicate that, in this study, this has not been a substantial factor. Illumination could be more challenging in the actual gallery than in the simulated conditions, but our experience in this study indicates that no parts of the image were critically under-illuminated and that we were successful at correcting moderate illumination unevenness. Coarser pixels, instead, appear to be more responsible for the observed decreased accuracy. Large and uniform targets such as those of the standard reference have a crown of incorrect labels that put in evidence the effect of mixed pixels. While not so conspicuously, this effect is certainly present in other targets and must be more pervasive in those targets of smaller size, such as the case of cassiterite. Fragmentation probably makes the average size of ore minerals in the hand samples smaller than in the actual mine walls, and thus the effect of pixel coarseness at producing spurious mixtures might have been over-emphasized in this study. Notwithstanding, the fact that the minimum size of ore patches being worth extracting can be as small as 1 cm 2 , indicates the interest of using a hyperspectral system with higher spatial resolution for forthcoming studies in the actual mine gallery.
Spatial resolution can also be improved by co-registration of hyperspectral images to conventional digital photography ("RGB") and subsequent image fusion. Co-registered RGB images of high resolution are useful to add textural information from the high-resolution RGB image to the classification processing, which usually involves a previous segmentation step. Actually, a crude co-registration using 2D coordinates only was used in this study as an aid for photo-interpretation, but the correct co-registration requires 3D coordinates and thus generating a digital surface model [69], an involved task that will be worth addressing in the forthcoming study of the actual mine face.
Our results demonstrate the feasibility and interest of mapping materials of gallery mine faces in Sn-W deposits by analyzing hyperspectral images in the 950 nm -1780 nm range through machine learning methods. A sequence of these maps at given time intervals as the excavation progress would improve orebody assessment and document the structure of the deposit as a tomography, opening the way to detailed studies of the spatial distribution of ore mineralization and the evaluation of geologic models of the deposit. Our study under simulated gallery conditions encourages forthcoming studies based on deploying ground-based hyperspectral systems in the actual mine gallery, which still will have to face challenges in illumination and surface modelling.  Acknowledgments: We thank Valoriza Minería, in particular Ivan Losada, for providing the mine 13 hand samples, as well as Marc Campeny-Crego and Fernando Tornos for their useful comments. 14 Hyperspectral imaging was conducted by the Laboratory of Hyperspectral Imaging of GeoSciences 15 Barcelona, CSIC.

Supplementary Materials
16

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the 17 design of the study; in the collection, analyses, or interpretation of data; in the writing of the manu-18 script, or in the decision to publish the results. 19