Identification of intrusive lithologies in volcanic terrains in British Columbia by machine learning using random forests: The value of using a soft classifier

Northern British Columbia, Canada, predominantly consists of a series of intraoceanic island arc terrains that were accreted onto ancestral North America in the Mesozoic (e.g., Bond and Kominz, 1984; Gabrielse et al., 1991; Monger and Price, 2002; Nelson and Colpron, 2007; Johnston, 2008). The Kliyul Cu-Au prospect is located in the northern end of one such terrain, known as Quesnellia (Figure ¹). Quesnellia consists of high-potassium calc-alkaline to shoshonitic submarine volcanics and sediments known as the Takla Group (Lord, 1948; Monger, 1977). The Hazelton Group andesitic volcanics, which are abundant at the Kemess porphyry Cu-Au deposit north of Kliyul (Rebagliati et al., 1995; Gordee et al., 2004) unconformably overlay the Takla group. The Hogem Batholith is located to the south of Kliyul, and it is a large, fertile, Late Triassic to Early Cretaceous igneous body that hosts other porphyry deposits such as the Lorraine (Woodsworth, 1976; Nelson and Bellefontaine, 1996).

A new geologic map of the Kliyul property (Figure ²) was created following a 350 m gridded rock sample campaign of the entire property, combined with targeted 1:10,000 scale mapping in August 2017 by First Quantum Minerals and AuRico Metals personnel. This map, showing features as described below, will constitute the reference geology map with which Random Forests (RF) outputs will be compared. The Kliyul property consists of volcanosedimentary strata, a series of discrete to large intrusive bodies, and several prominent faults (Table 1; Figure ²). Quaternary glacial-fluvial cover fills the low-lying valleys. Volcanosedimentary strata at Kliyul consist of two units of the Takla Group volcanics, Goldway Peak (GP) and Kliyul Creek (KC) Groups. These are equivalent to the same named units in the British Columbia Geological Survey government map of the area (Schiarizza, 2004a, 2004b; Schiarizza and Tan, 2005). The oldest unit is KC, which consists of four subunits: intermediate to mafic volcanics and volcanoclastics (KCv), carbonate-rich sediments and volcanics (KCc), felsic volcanics and volcanoclastics interlayered with carbonate-rich sediments (KCfc), and siliceous sediments (KCs). Strata have dominantly low-angle northeast dips. Overlying KC and all of its subunits are the GP volcanics. GP consists of more mafic to intermediate augite-phyric volcanic flows and breccias. A series of calc-alkaline, fine- to coarse-grained intrusions with variable pre-, syn- to postmineralization and alteration relationships cross-cut the KC and GP volcanosedimentary strata (Figure ²). Age dating is poorly constrained, but historic and recent efforts estimate a Late Triassic to Early Jurassic timing (~217 Ma). Three apparently long-lived structures influence the map patterns at Kliyul: (1) the north–northwest-striking KC fault with an apparent dextral sense of shear with sinistral and dextral fault splays; (2) the east-northeast-striking valley fault with poorly constrained, but apparent normal, north-block-down kinematics; and (3) the north-striking dextral Dortatelle fault that had a strong asymmetric influence of rock types to the west versus east of the fault, with a mylonitic-like shearing of the rocks in the southwest.

We use RF ( Breiman, 2001), a supervised machine-learning algorithm, to classify each pixel in the study area according to the lithology. A number of algorithms have been tested and applied to geologic mapping enterprises. Examples include variants of artificial neural networks (e.g., Barnett and Williams, 2009; Cracknell et al. 2014; Grunsky and Kjarsgaard, 2016), support vector machines (Yu et al., 2012; Cracknell et al., 2014), and RF (Cracknell and Reading, 2013, 2014; Cracknell et al., 2014; Rodriguez-Galiano et al., 2014; Harris and Grunsky, 2015; Kuhn et al., 2016, 2018). Whereas all of these choices have been successfully deployed in various scenarios, RF has been consistently identified as a good choice of classifier on grounds of performance, ease of use (in terms of few and simple tuning parameters and insensitivity to variable scaling), and ability to handle data with high intraclass variability and high interclass similarities (Breiman, 2001; Cracknell and Reading, 2013) by multiple studies. These include using RF alone (e.g., Cracknell et al., 2014; Harris and Grunsky, 2015; Kuhn et al., 2016, 2018) or in comparison with other algorithms (Waske et al, 2009; Cracknell and Reading, 2014). The accessibility of RF makes the methods described in this study widely applicable in the geoscientific disciplines, minimizing the barrier to entry imposed by requirements for specialized computing skills and equipment.

RF extends upon the classification and regression tree methodology (Breiman et al., 1984), addressing the high variance (Friedman, 1997) that can be associated with single decision trees (Murthy, 1998) by constructing an ensemble or forest of pseudounique decision trees (Figure ³). The performance of an RF is controlled by the strength of individual trees in the forest and the correlation between trees, with a forest of strong decorrelated trees being the ideal case (Breiman, 2001; Hastie et al., 2009). To minimize correlation between trees, randomness is introduced into the selection of subsets shown to each tree and in the selection of variables used to split decision nodes. Bagging (Breiman, 1996) is used to produce a random subset of input training data, of equal size to the training set provided to each classification tree. This has been shown, on average, to include approximately 63.7% of unique instances from the training set, with the remaining, “out of bag” instances, held aside for testing (Breiman, 2001). During learning, a randomly selected pool of variables of predetermined size is provided to each decision node. From this pool, the variable that produces the maximal improvement in homogeneity of the child nodes relative to the parent node, as measured by Gini impurity (Breiman et al., 1984 and described by Kuhn et al., 2018), is used to split that node (Figure ³). This split is performed at all nodes in each classification tree and, subsequently, all trees in the RF. Final classifications are assigned as the modal decision of the forest (or mean in the case of regression).

In this study, we explore a pragmatic approach for a relatively data-limited situation. We train an RF classifier using soil sample points defined by geochemical and geophysical data. This RF is subsequently used to classify all of the remaining samples in the Kliyul project area. The soft classification metric, CMP, and the more typical classification results are used to predict lithology and identify intrusions in the project area. Information entropy is used to define, objectively, uncertainty and complexity associated with classification. We use the FQM geologic map, which benefits from rock and soil sampling, in addition to detailed location and transect mapping, as a benchmark with which to compare classification results.

This study incorporates all available geochemical and geophysical datasets that encompass the project in addition to shuttle radar topography mission (SRTM) elevation data. Geochemical data include a suite of 49 measured elements (ICP21 + MS61 + XRF5000) for rock and soil samples collected at a 350 × 350 m spacing over the project. These have been preprocessed such that all data exhibiting a measurement below the detection limit were assigned a value of half of the detection limit for that element. Data sets comprising an excessive number of samples below the detection limit (in excess of 25% or where spatial distribution resulted in large areas with insufficient real data) were omitted. Geophysical data (Supplement 1; the full supplement can be accessed through the following link: S1.pdf) comprised reduced-to-pole (RTP) total magnetic intensity data and derived data sets (vertical and horizontal derivatives, tilt derivative, horizontal gradient magnitude, analytic signal, and analytic signal of the vertical integral), and airborne radiometric data (potassium, thorium, uranium, and total count), flown at a 100 m line spacing. The magnetic and radiometric data were gridded at a 100 m cell size using a bidirectional spline algorithm. Airborne electromagnetic data flown at 100 m line spacing were gridded at 100 m. IP/resistivity data have been acquired in this area; however, their limited extent did not provide adequate coverage of the entirety of the study area, and thus they were not included in this exercise. All data sets were resampled to a 50 × 50 m grid, a compromise between the higher resolution afforded by the magnetic and radiometric data and the lower resolution geochemical and electromagnetic data and compiled as a matrix in the form of easting, northing, V₁, V₂,…, V_n (where V is a given variable/data set). To preserve conditional independence, where pairs of variables correlated in excess of 0.8 (Pearson’s correlation coefficient), one of the pairs was removed. Where a variable was a data set exhibiting any indication of poor quality, for example, excessive missing data or readings below detection limit, that data set was omitted. Where multiple correlations were present, data sets were removed in such a way as to maximally reduce the dimensionality of the data set. If neither of these criteria were encountered, but only where a correlation in excess of 0.8 was shown, data sets were removed based on subjective geologic utility. The data sets analyzed as important for best possible results (as described in the following section) are shown in Table 2.

Pixels comprising samples for which a lithology had been assigned (Table 1) were taken as training data. The remaining data were held for classification by the trained RF (a flowchart is given in Supplement 2; the full supplement can be accessed through the following link: S2.pdf). RF is prone to bias in favor of a numerically dominant class (Hastie et al, 2009). To mitigate this tendency, classes were balanced to 50 samples per class, through either bootstrapping (where fewer than 50 samples were available) or random decimation (where more than 50 samples were available). The composition of intrusions is highly variable and comprises felsic, intermediate, and mafic units. However, due to the impractically small number of samples for each, this class was combined. We confirm that the training data adequately sample the intrusion class and that training data samples and intrusion types are colocated. Data partitioning considerations in the spatial context, in a general sense applicable to machine learning, are outlined by Vucetic et al. (1999). The main objective of this study is to determine if RF can identify intrusions from data typical of early-stage exploration, despite the heterogeneous class definition (and also noting that we aim to identify an indicator lithology at this early stage, not the ore deposit itself). It is understood that better, and worse, results could be obtained through other means of retaining or adding synthetic samples to achieve a balanced class size. For the purpose of the present study, this process was restricted to the balancing described earlier, similar to the method used by Kuhn et al. (2019) using only observed data values as opposed to imputing from an estimation of a data set’s probability density function or other common methods. In this study, at an early stage in the exploration cycle, the notable limiting factor is the resolution of the input data; however, it would be expected that the results would improve if additional training data were available (Cracknell and Reading, 2013). Samples comprising the training set are shown in Figure ⁴.

Following the balanced sampling process, RF is used to rank the importance of variables according to the mechanism described by Breiman (2001). Under this strategy, each variable is permuted and shown to the RF (Demsar et al., 2013). The variable that produces the largest difference in classification accuracy is most important, and the variable that, when permuted, produces the smallest difference is least so. In this study, we seek to produce a classifier that benefits from the additional dimensionality permitted by a machine learner as compared to more conventional methods (e.g., successive comparison of scatterplots or a manually weighted geographic information system [GIS] analysis) while producing a result that maximizes interpretability by the end user. As such, we seek to reduce the number of input variables to the minimum required to produce the best result. This is possible because RF tends to produce peak accuracy with the inclusion of a given number of variables, after which the results remain stable or may even marginally deteriorate; additional variables at this point are redundant. Prior studies deploying RF in a similar fashion for lithology classification (e.g., Cracknell et al., 2014; Kuhn et al., 2016, 2018, 2019) have shown a tendency for this to occur at between 8 and 15 variables, although this may change as additional studies are performed.

For this study, we used an RF comprising 500 classification trees with no pruning or growth restrictions (see also Supplement 2, the full supplement can be accessed through the following link: s2.xlsx for other implementation parameters). RF is not prone to overfitting with additional trees, instead reaching a stable error minima (Breiman, 2001). This number of trees is well above what is likely required and represents a safe choice for geoscientists looking to replicate the methodology of this case study without any risk of using insufficient trees. These parameters were duplicated during ranking. All variables were ranked via the process described previously. These variables were then successively used to build a classifier that was in turn assessed for accuracy. We used tenfold cross-validation (James et al., 2013) in conjunction with a backward recursive process to determine the cross-validation accuracy of the RF, drop off the lowest-ranked variable, and rerank the variables. This procedure was repeated for all variables. Best results of 83.5% cross-validation accuracy were achieved with the inclusion of the top-ranked 15 variables (Table 2). Objectively determined variables required for best results comprised geophysical and geochemical data types. Radiometrically derived K and the RTP total magnetic intensity were both ranked highly, whereas common rock-forming elements Fe and Ca were the highest-ranked geochemical variables. Spectral remote sensing data did not produce an effect on the results nor did any magnetic derivatives (see the “Discussion” section). Figure ⁵ shows examples of top-ranked geophysical and geochemical data sets.

Using a confusion matrix (Table 3: a comparison of actual and predicted class labels for each training data point, in this case made during tenfold cross-validation), we can see that several classes were consistently classified correctly. The KCv class performed the most poorly, being commonly misclassified as intrusive or GPa. Intrusions were most often misclassified as KCv. The RF used to produce this result was selected as the final classifier for use on all remaining data. Although the relationships shown in this confusion matrix (based on cross-validation results) do not necessarily translate to classification on blind data, this does give some indication of the potential strengths and weaknesses of the selected classifier.

The results of the RF classification show 73% overall consistency (Figure ⁶) with the FQM geologic map (Figure ²). Results were more consistent with the FQM geologic map where a lithology was present in larger domains, whereas results appear less robust for narrower zones of a given class, such as where lithologies appear to wrap around the topography due to low-angle bedding dips. This is in part a function of the resolution of the geochemical input data and is discussed in the next section. Information entropy (H) identifies many of the class boundaries in the project area, whereas normalized information entropy (H_norm) shows a high per-pixel uncertainty throughout the project area (Figure ⁶ and ⁶). Both variants highlight some areas for which RF is classified with low uncertainty, including, for example, the large intrusion in the north of the project. CMPs provide a more detailed description of the information summarized in the calculation of H and H_norm. It can be seen (Figure ⁷) that there are clear domains classified as intrusive in the north, GPa in the east, SIVC in the west, and KCv through the center of the project. It can be seen, however, that as demonstrated by H_norm (Figure ⁶), multiple classes were similarly probable particularly in the center of the project area. It is worth noting that many discrete zones of elevated probability of intrusive class, the target of economic significance, are identified. In many of such cases, these areas were not classified in the final RF lithology map as intrusive. The value of using the soft classifier, elevated CMP, is clear from this analysis.

Many of the broader domains in the RF lithology map resulting from the RF classification are consistent with the FQM geologic map. The most notable shortcomings in the RF lithology map are due to a lack of resolution in regions where the mapped lithotypes have more detail than the geochemical sample spacing that is reflected in the pixel size used in this study. Large intrusions were identified in final classification; therefore, many of the subtler intrusive features were either not identified or mapped without adequate resolution to be of substantial use in targeting or follow-up work. There are several likely causes for this, foremost of which is the resolution of input data, a result of the sampling interval of each data set, and the effects of interpolation of those data. The limitations of gridding geochemical data are well understood as a potential source of error in the result. If attempting to predict the nature of a region not represented by a sample, a value must be imputed or interpolated. Interpolation of geochemical data is a technique used pervasively throughout the mineral exploration industry. As such, this study provides a faithful indication of performance on the data likely to be encountered in a typical industry case example. More conventional means of interpreting magnetic data (e.g., Salem et al., 2008; Isles and Rankin, 2013) could potentially be used to improve further RF interpretation. Whereas in this case all magnetic derivatives were omitted due to poor ranking on the ability to define lithology, these could be used to better define structure and boundaries throughout the project area.

Additionally, it is entirely plausible that the use of more sophisticated imputation strategies, in place of simple bootstrapping, could improve training class definition and thus the results. Particularly where available real samples facilitate the accurate prediction of a given class’ true distribution. This principle extends to several aspects of the study methodology. This study was intentionally designed to demonstrate a rapid, effective approach easily adopted by the wider geoscientific community, and it produced good results. It is possible that further tuning of RF parameters, or more advanced feature ranking and selection (e.g., Lundberg and Lee, 2017) might be useful for studies in a more data-rich context.

The use of an ensemble of unique classifiers (classification trees) is the key feature of RF. Not only does this have the important benefits of improved accuracy and the ability to respond well in dealing with high intraclass variability and interclass similarity, but this also facilitates a detailed analysis of classification results. Plotting individual CMPs indicates where each class was most likely to occur (Figure ⁷), and in this study is an approach that leads to valuable insight for the exploration program going forward. Information entropy responds to uncertainty and complexity: the distribution of probabilities and the number of classes possible at a given sample. In this case, observing H (Figure ⁶) highlights boundaries between classes and other regions of more complex geology. This coincides with similar observations of the behavior of H (e.g., Kuhn et al., 2016, 2019). The term H_norm (Figure ⁶), a measure of how closely each instance approaches its own maximal uncertainty, is high throughout much of the study area. This suggests that there is a high likelihood of the classification being incorrect (e.g., Kuhn et al., 2016, 2018).

Given the stated goal of identifying and defining intrusions, potentially of economic significance to the project, we focused on the CMP of this unit (Figure ⁷). Even where this is not the majority class and therefore not the result of final RF classification, an elevated CMP can be seen (Figure ⁸). Contouring intrusion CMPs at multiple levels from above the background level (CMP of 0.125, a random response) to several times the background level facilitates a trade-off between correctness and completeness. Many areas exhibit a sharp transition between background levels and an intrusion CMP of three times higher than background level (Figure ⁸). Furthermore, the location, form, and trend of these zones of elevated CMP better define the intrusive lithologies than the final RF classification. When compared to the FQM geology map, we can see that most intrusions have been captured (Figures ⁸ and ⁹), whereas the false-positive rate is low. This result was achieved in areas where no sample of the intrusive class was used in training data, indicating that this method may be viable when using sparse and/or incomplete training data. In this case, intrusions were predicted where no detailed mapping has taken place and their presence in the broader volcanosedimentary packages is likely. This includes the prediction of discrete but substantial intrusive bodies under regions mapped as glacial till at the surface (Figures ² and ⁸).

Neither mapped nor predicted intrusions (Figure ⁹) show a consistent response to single observables, geochemical or geophysical (Figure ⁹ and ⁹), although a small number of units, SIVC, for example, can be separated on the basis of magnetic character, whereas a small number of other units are discernible in radiometric data. Our results indicate that use of the RTP was more effective than derived products such as the first or second vertical derivative, analytic signal, tilt derivative, or horizontal gradient magnitude. This is consistent with previous results (e.g., Kuhn et al., 2018) and can be explained by the fact that the method used in this study relies on a “per pixel” classification. Magnetic derivatives, ubiquitous tools for mapping geology, through the determination of texture, geometry, and edge definition, are effectively high-pass filters. As such, pixel values from the derived products are less diagnostic of lithotype than those of the RTP. Some intrusions observed in this data set can be observed in the magnetic data, and more still could be identified through more advanced analysis. Many intrusions, however, lack a magnetic response or do not show a consistent magnetic response, likely due to signal interference with surrounding country rock. Conversely, some lithologies lacking a diagnostic magnetic response do exhibit a characteristic response in K (radiometric) or Ca data sets (Figures ⁴ and ⁹).

Electromagnetic data were omitted during variable ranking due to an inability to classify lithology and as such were not used in the classification. Subsequently, they show limited correspondence to the classification results. Induced polarization, although not included due to the limited spatial extent, where present, does not correspond well with the classification results. This is not surprising given that these electrical/electromagnetic methods in this glacial till-covered environment are responding to many factors (e.g., variable weathering, ice, water content, and high-relief topography). This highlights the benefit of a machine-learning/data fusion approach, making use of many types of data (Figure ⁹). The situation in which elevated CMPs do not “win” the final classification is partly due to the overlapping response across class boundaries in the training data and extensive adjacent zones. Further contributing factors are the overly smooth response of interpolated data, the sampling interval of geochemical data relative to the width of intrusions, and the inherent resolution of the methods themselves, both as a function of the physical and chemical expression of the rocks and the element distribution in eroded and often transported soils. Taking note of CMPs is important because it demonstrates that, although the expression of subtle intrusions may not result in a majority decision in the RF classification, they can be detected as a subtler expression in the CMP, that is, soft classification metrics produced during the RF classification. As such, we assert the importance of using CMPs in any similar scenario in which the goal is the identification of an important rock unit (or any key indicative feature) as opposed to the overall accuracy of an RF lithology map. These products can be used in conjunction with or to further guide the deployment of more expensive and labor-intensive mapping and geophysical and geochemical acquisition campaigns.

Locating intrusions is a key component of exploration for porphyry-style Cu ± Au ± Mo deposits. In porphyry hosting terrains, such as British Columbia, Canada, where this study is located, lithologies are difficult to discriminate due to similarities in the numerous generations of volcanic, volcanoclastic, and intrusive units present, often with similar provenance.

In this study, we used systematic rock and soil sampling points as training data points. Additional detailed geologic mapping produced a new FQM geology map, against which the results of our RF classification could be compared. The similarity between many of the lithological classes present resulted in competing class probabilities and an erratic classification. Intrusions, particularly those of a more discrete nature, were inconsistently predicted by the aggregate classification. This was due to their limited extent, relative to the resolution of the underlying data samples, causing their expression to be overwhelmed by the spatially larger classes in which they are emplaced.

Closer examination of CMPs, rather than the aggregate classification, indicates that there were many locations where the probability of an intrusion being present was significantly elevated above background. Indeed, a large proportion of mapped intrusions were captured by areas of elevated intrusion CMP. Additionally, areas were identified showing an elevated probability of membership to the intrusion class that was not mapped in the FQM geology map.

For the task of identification and location of intrusions, intrusion CMP, that is, not the overall RF lithology map, is the more useful product. Use of this soft classifier has the potential to yield valuable insight, especially at the early stages of exploration. We anticipate that this understanding will find use as a rapid near-real-time tool by exploration teams wishing to predict intrusion locations in order to target and prioritize field activities. More generally, we encourage the use of RF CMPs to gain insight into the occurrence of the classes of greater significance in any data-driven research challenge.

We would like to thank First Quantum Minerals Ltd. for permission to access data. We thank Chris Wijns for support and discussion regarding data and results. Stephen Kuhn is supported by a Tasmanian Graduate Research Scholarship (TGRS) from the University of Tasmania. This research was conducted as part of the ARC Industrial Transformation Research Hub for Transforming the Mining Value Chain (project number IH130200004) at the Centre for Excellence in Ore Deposits and Earth Science, University of Tasmania. The views expressed herein are those of the authors and are not necessarily those of the Australian Research Council. Preprocessing, interpolation, and plotting were performed using Geosoft Oasis Montaj and ESRI ArcGIS.

Data associated with this research are confidential and cannot be released.