Automated classification of metamorphosed coal from geophysical log data using supervised machine learning techniques

International Journal of Coal Geology 214 (2019) 103284

Contents lists available at ScienceDirect

International Journal of Coal Geology journal homepage: www.elsevier.com/locate/coal

Automated classification of metamorphosed coal from geophysical log data using supervised machine learning techniques

T

⁎

Kane Maxwell , Mojtaba Rajabi, Joan Esterle School of Earth and Environmental Sciences, University of Queensland, QLD 4072, Australia

A R T I C LE I N FO

A B S T R A C T

Keywords: Altered coal Supervised machine learning Geophysical log data Random forest Coal resource

Identification of the extent of heat altered coals is important for coal mining and resource estimation because alteration directly affects key economic properties such as ash content, volatile matter, specific energy, sulphur, free swell index and beneficiation characteristics. The most reliable method to identify altered coal is through geochemical analysis of drill core samples; however this method is time consuming and costly. Cheaper alternative methods to identify altered coal is by macroscopic observation of non-core drill cuttings, but this is less reliable than core. Geophysical wireline log data is also used to identify altered coal, but with variable success. Over the past decades, Machine Learning methods have proven popular to automatically classify a variety of lithotypes from geophysical log data. However rarely has there been focus on predicting altered coal. In this paper we apply the most recent machine learning methods which include gradient boosted machines, random forests and artificial neural networks to automatically predict altered and non-altered lithotypes using geophysical log data. We use a massive data set comprising of 1230 samples from 263 boreholes from a highly intruded deposit in the Bowen Basin, Eastern Australia. To do this, we calculate each sample's distance from intrusion and predict their relative density from geophysical log inputs including gamma ray, caliper and compensated density. We then train our machine learning methods on 80% of the data to predict alteration class using the calculated distance to intrusion, predicted relative density, and geophysical log gamma ray, as primary inputs. We evaluated each of these machine learning methods on the remaining 20% of the data to determine the best performing model. Finally, using the best performing model we further split the altered and non-altered classification into lithotypes using a decision tree based on geological knowledge of the case study area. The results indicate that of the machine learning algorithms the random forest produced the best results with only 11 misclassifications across the entire data set of 1230 samples which represents < 1% error.

1. Introduction Magmatic intrusions, and associated contact metamorphism, in stratified coal basins can change the physical and chemical properties of coal seams at different degrees, depending on the morphology and type of intrusion (Cooper et al., 2007; Golab et al., 2007; Rimmer et al., 2009; Salmachi et al., 2016). In coal seams, contact metamorphism may directly affect their key economic properties such as ash content, volatile matter, specific energy, sulphur, free swell index and beneficiation characteristics (Presswood et al., 2016; Quaderer et al., 2016; Rimmer et al., 2009). Alteration may also influence the self-heating potential of coal (Shi et al., 2018), increase the chance of gas outburst (Jiang and Cheng, 2014; Yao and Liu, 2012) and affect coal and surrounding strata geotechnical characteristics (Lu et al., 2016; Wu et al., 2018). For these reasons, accurately defining the extent of altered coal

⁎

is important for safety, resource estimation and mine design. Laboratory analyses of coal samples from drill cores is considered a reliable method for the identification of fresh and altered coals. For example, several studies from coal deposits across the world have highlighted that altered coal has an increased density and ash, and decreased volatile matter and fixed carbon near intrusion (Finkelman et al., 1998; Fredericks et al., 1985; Rimmer et al., 2009). However, acquisition of core samples and laboratory analysis are not always possible as they are expensive and time consuming. Due to this, these analyses are often sparse making it difficult to identify the extent and degree of coal alteration. Where laboratory analysis is not available, fresh and altered coals are attempted to be differentiated by macroscopic observation of drill cuttings obtained from cheaper rotary noncore drilling. Attempts are then made to manually align these observations with geophysical log data which is commonly acquired. The

Corresponding author. E-mail addresses: [email protected] (K. Maxwell), [email protected] (M. Rajabi), [email protected] (J. Esterle).

https://doi.org/10.1016/j.coal.2019.103284 Received 27 June 2019; Received in revised form 10 September 2019; Accepted 12 September 2019 Available online 14 September 2019 0166-5162/ © 2019 Elsevier B.V. All rights reserved.

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Finally, we subdivide these altered and non-altered class predictions into lithotype groups using a decision tree. We specifically incorporate predicted relative density and distance to intrusion calculations in order to improve the accuracy of classification and to avoid reliance on the resistivity log.

combination of geophysical logs which are best able to identify altered coal are density and resistivity (Johnson et al., 1963; Roslin and Esterle, 2016; Thomas, 2012) however, resistivity logs are uncommon in the open-cut coal mining industry and most boreholes do not have resistivity log information. In the past decades, different statistical and machine learning approaches have been used to automatedly predict lithotypes directly using a combination of geophysical well log data in petroleum wells (Al-Anazi and Gates, 2010; Corina and Hovda, 2018; El Sharawy and Nabawy, 2016; Kuroda et al., 2012). Generally, identification of coal intervals from other non-coal lithologies is not a difficult issue, and several studies have already revealed that different machine learning approaches can accurately differentiate coal from other non-coal lithologies due to its distinguishing low density and gamma ray ranges (Horrocks et al., 2015; Roslin and Esterle, 2016; Schmitt et al., 2012; Xie et al., 2018; Zhou and O'Brien, 2016). However, altered coals are more difficult to distinguish due to their higher density which makes them easily confused with other high ash coal or non-coal lithologies. For example, Roslin and Esterle (2016) using a clustering method showed that the density and gamma ray ranges of altered coals can be coincident with dull inertinite rich coal, and hence required resistivity log to correctly distinguish altered coal. Among different types of machine learning methods, Random Forests (RF), Gradient Boosted Machine (GBM) and Deep Neural Network (DNN) algorithms have recently proven the most accurate at classifying lithofacies from geophysical log data (Bhattacharya and Mishra, 2018; Gu et al., 2018; Horrocks et al., 2015; Xie et al., 2018). Therefore we specifically choose these most accurate methods for predicting altered and non-altered lithology groups. To do this, we use a massive dataset including geochemical analyses of 1230 coal samples and geophysical log data in 263 boreholes in the Leichhardt Seam of the Bowen Basin in Eastern Australia (Fig. 1), of which 80% of the database is used for model constructions and 20% is used for testing the modelling results. In order to predict altered coals from geophysical logs we first use geochemical data to label samples as altered or non-altered, we then predict relative density and calculate distance from intrusion for each sample. We then use these predictions and distance calculations along with other geophysical log data to predict altered and non-altered material using GBM, RF and DNN to find the best performing model.

2. Study area The study mine is located in the eastern part of the Bowen Basin within the Nebo Synclinorium, in Queensland (Fig. 1). The main mining targets in the basin are the Late Permian Moranbah and Rangal coal measures and their stratigraphic equivalents (Fig. 1). The Rangal Coal Measures (typically 150 m thick) contain low to moderate ash coal seams and are the focus of this study. Within the Rangal Coal Measures, the primary seams in the region are the Vermont, Leichhardt and Phillips Seams, which have variable character (Roslin and Esterle, 2015). The only economic seam at the study mine is the Leichhardt Seam which is mined and marketed as low volatile pulverized coal injection coal for use in steel making. A major geological feature at the study area is the presence of Cretaceous age lamprophyres and dolerites (Ritchie, 2010) which heavily intruded the eastern flank of the deposit causing alteration of coal. The primary morphology of these intrusions are stratified sills and subvertical in-coal climbing sills/dykes with typical thicknesses of 0.6 m (Fig. 2). This morphology is typical of that described in Chen et al. (2014), and Yao and Liu (2012). Altered coal in the study area is differentiated due to its much lower volatile matter content, slightly higher ash, and much higher density. The low volatile matter of the altered coal is seen as favourable as it can be utilised in coal blending to lower the average volatile matter content of other coals. However, marked increase in relative density decreases beneficiation properties resulting in low product yield. Altered and nonaltered coals are defined in resource estimates through either geological or geochemical characterisation. Where geochemical analysis is available, coal is currently classed as either altered or non-altered based on an 8% volatile cut off value according to mine site guidelines. Where geochemical analysis is not available, coal is categorised into altered and non-altered based on drill core or cuttings observation and subsequent adjustment to geophysical log signatures based on internal

Fig. 1. Approximate location and generalised Permian coal stratigraphy of the study area relative to broader Bowen Basin. 2

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 2. Photograph of high-wall at the study mine showing the four main morphology types of intrusion including (1) coal roof sill; (2) coal floor sill; (3) in coal sill; (4) climbing sill/dyke.

mine site guidelines. The unit thicknesses of these categorisations are used in subsequent 3D geological modelling process and ultimately determine the volume of altered coal and non-altered coal resource. Because drill-holes with geochemical analyses are sparse, a large majority of the resource volume is dependent on accurate geological categorisation of coal in non-core holes, aided by geophysical wireline logs.

trees are not too large, it is robust to correlated predictors, and is insensitive to outliers (Breiman, 2001). Disadvantages are that predictions in production environments can be slow and interpretability (as compared to single decision trees) is less intuitive.

3. Theory of machine learning methods

Gradient Boosted Machine (GBM), like the RF algorithm is an ensemble method comprised of many base (weak) learners. The method used in this paper was proposed by Friedman et al. (2001) and the specific algorithm follows the method as outlined in Friedman et al. (2001). Distinct to RF, the GBM builds many (usually) shallow decision trees sequentially, and each successive tree aims to correct errors in the previous ensembles of decision trees, which is often referred to as boosting (Friedman, 2002). The measure of error is determined by the operator and is termed a loss function. An example of a loss function for regression problems might be squared error, and for classification might be accuracy. The way in which successive ensembles of trees try to further minimise the chosen loss function is called gradient descent (Friedman, 2002). In this case, the goal of gradient decent is to optimise the parameters of each successive decision tree in such a way that the combined predictions of all decision trees result in the lowest error, and optimally has zero error. In order to control how much error is minimized after each step, a parameter called learning rate is introduced. Most important parameters to the GBM algorithm are the depth of the decisions trees, the number of decision trees to build, the learning rate the sampling technique, and the loss function. For example, when the depth of trees is very deep, and the number of decision trees is large the algorithm can be computationally intensive and can easily over-fit the data. Overfitting can somewhat be mitigated by small learning rates but can result in long training times, so appropriate balance between the number of trees and learning rate is required (Friedman, 2002). Sampling without replacement can reduce training times but can introduce randomisation into the gradient decent process and will not allow for perfect fit so also should be well considered. Finally, since the loss function can be dependent on the research area poor selection can result in sub-optimal predictions. GBM has recently shown superior performance to all other algorithms for lithofacies prediction (Society of Exploration Geophysicists, 2016; Xie et al., 2018; Zhang and Zhan, 2017). Key advantages of the algorithm are that it can be highly accurate and can generally outperform random forest if carefully tuned, it can take most any loss function and use any base learner algorithm making it highly flexible, pre-processing of features (such as normalisation and imputation) are usually not required. Primary disadvantages are that it can be prone to overfitting, there are many parameters which can affect the performance of the algorithm and choosing these can be difficult or computationally intensive.

3.2. Gradient boosted machine

Generally, the majority of machine learning approaches fall into two categories including unsupervised and supervised methods (Tan et al., 2018). When classifying lithotypes, unsupervised approaches are used when the reliability of lithology identification by geologists is poor or absent and geophysical well log data is abundant. Supervised approaches, which we used in this study, are used when the reliability of lithotype descriptions is high (for example from fully cored holes with geochemical analysis and lithotypes described by highly trained or skilled geologist) and there are enough examples of these data. As mentioned above, we used three supervised approaches which have recently proven to be the most accurate at classifying lithofacies from geophysical log data (Bhattacharya and Mishra, 2018; Gu et al., 2018; Horrocks et al., 2015; Xie et al., 2018). 3.1. Random forest Random Forest (RF) algorithm is based on a decision tree that is non-linear algorithm and can be used for both regression and classification (Breiman, 2001). A single decision tree can suffer from high error and large variance and is prone to overfitting, however the RF algorithm aims to mitigate these issues by building and combining the predictions of numerous decision trees from random samples of data to find the optimum prediction (Breiman, 2001, 2017). The aim of this random sampling is to reduce the variance of individual trees and to create trees that are uncorrelated which improves predictive performance (Breiman, 2001). The method of combining individual decision trees (termed base or ‘weak’ learners) into a single learner (termed a meta or ‘strong’ learner) is known as ensemble (Dietterich, 2000). The random sampling method is termed bootstrap aggregation (Breiman, 1996) which aims to create new slightly different data sets but with the same distribution as original. In addition to bootstrap aggregation, random features are also selected to train individual trees which distinguishes it from other tree-based methods (Breiman, 2001). The RF algorithm has recently been successful for lithofacies prediction from geophysical log data and has been shown to perform similar to Gradient Boosted Machine and better than other methods such as Support Vector Machines and Naïve Bayes algorithms (Bhattacharya and Mishra, 2018; Xie et al., 2018). The key advantages of the RF algorithm are that it can achieve good performance with minimal tuning, it is less prone to overfitting, it provides an estimate of error (from out of bag samples) without need for cross validation, it can be fast to train if number of 3

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

prediction of altered materials and interpretation of predicted lithology groups using a decision tree. In this section, we briefly describe the materials and explain each of the pre-modelling steps that we used to prepare data for modelling.

3.3. Deep neural networks Neural networks were motivated by the concept and structure of the biological brain however their architectures and processing are distinct (Svozil et al., 1997). They comprise of interconnect nodes (referred to as neurons) of input, hidden and output layers. The input layer read the data in, and nodes in the hidden layer updates weights and apply a function in attempt to make the output as close to the true input value as possible (Svozil et al., 1997). Deep Neural Networks (DNN) are distinguished from simpler networks by fact that they have many hidden layers (Schmidhuber, 2015). They may also be referred to as multi-layer perceptron or feedforward artificial neural networks (Schmidhuber, 2015; Svozil et al., 1997). The DNN algorithm as implement in this study is a multilayer feedforward network which uses stochastic gradient descent with back-propagation (Candel and LeDell, 2019). Gradient decent is used to optimise the weight values of hidden layers in a forward direction, while back-propagation allows us to update these weights in a backward direction. While neural networks have previously been successful in predicting lithofacies (Dubois et al., 2007; Horrocks et al., 2015; Xie et al., 2018) the primary disadvantage of this method is that they have many parameters which need to be specified by the operator. This makes it more difficult to find the optimum combination of parameters for the problem space and using similar values from previous studies may not be appropriate due to natural heterogeneity of geological environments and availability of geophysical log data.

4.1. Computer software and programming languages All data manipulation, exploratory data analysis, machine learning modelling and data visualisation was written in R programming language (R Core Team, 2017) using the graphical user interface RStudio (RStudio Team, 2016). All machine learning classification models were created using the H2O framework which is an “open source, inmemory, distributed, fast, and scalable machine learning and predictive analytics platform” (The H2O.ai team, 2019). This framework was accessed through the R package (‘h2o’) (The H2O.ai team, 2018) and was chosen as the primary platform for its speed (due to distributed processing), access to current ‘state-of-art’ algorithms, and ability to easily compare between machine learning algorithms. 4.2. Data sources Our original data set to conduct this study comprised of 3996 coal samples with geochemical analyses from 537 wells with different types of geophysical logs from Leichhardt Seam of the Bowen Basin. Following our data preparation step, which is described in the next section, we used different types of geophysical logs from 263 boreholes and 1230 coal samples to construct our machine learning models. The spatial location, geologist lithological descriptions, geochemical and petrophysical properties of samples were compiled from multiple data sources. Well collar locations were provided in the mine site survey datum and were captured by a qualified surveyor using a differential global positioning system with an accuracy of 0.2 m. Lithological descriptions of the study coals were captured by a qualified geologist using internal standard of the study mine, in conjunction with the

4. Materials, methods and pre-modelling steps As illustrated in Fig. 3, the full methodology for classifying altered and non-altered lithology groups comprised of multiple steps including pre-modelling and modelling steps. Pre-modelling steps included data preparation, calculation of distance from intrusions, and classification of samples. Modelling steps included prediction of relative density,

Fig. 3. The proposed multistep methodology for predicting alteration class and lithology groups in this study. The grey boxes indicate the pre-step modelling while the white boxes show the modelling attempts. 4

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

intercept data based on the geologist's log descriptions and consider these as reliable because intrusion is easily differentiated compared to coal and altered coals. For example, intrusions have distinctly highdensity signature compared to coals and can easily be identified and corrected accurately to geophysical logs. Due to this the identification of intrusion could be automated however for our study we simply rely on geologist log data. It must be noted that because the prediction of altered and non-altered coals is reliant on these calculations, the quality of the predictions will be affected by any inaccuracies in geologist logs. Further, an alternative method could be to calculate a samples distance from intrusion which has been geologically modelled in 3D however the accuracy of these distance calculations would be reliant on the robustness of the 3D model. Therefore, in our case we consider the calculation from borehole intercept data less onerous. In order to calculate the sample distance from intrusions (D), we used a nearest neighbour function based on the libnabo library, preferred for its computational efficiency (Elseberg et al., 2012), which was subsequently adjusted to account for sample thickness Eq. (2). Following the calculation of D, the average value of select geochemical and geophysical properties for each sample in the study mine were plotted relative to their calculated distance to intrusion (Fig. 4). This was done to gain an understanding of the scale and impact of intrusion on these properties to assist in determining important parameters for subsequent machine learning prediction.

CoalLog Geology and Geotechnical Training Manual (Green, 2015). Coal geochemical data comprised of proximate analyses (AS 1038 Part 3) and relative density (AS 1038.21.1.1-2002) conducted on borehole core by laboratories accredited by the National Association of Testing Authorities, Australia (NATA). For each well, geophysical log data was provided at 1 cm sample increments and was acquired primarily by a single vendor. Provided geophysical logs comprised of caliper, short spaced, long spaced compensated and VECTAR processed density, gamma ray, short and long spaced laterlog resistivity and porosity calculated from density. Depth matching of core sample intervals to geophysical log data was done manually by qualified geologists using internal standards of the study mine, which are based on guidelines in Firth and Elkington (1999). 4.3. Data preparation Sample data with geochemical analysis was merged with the 1 cm increment geophysical log data matched on depth. Where no geophysical data were available, these samples were removed from our database. In addition, geophysical log values within 10 cm of sample boundaries were removed to reduce or eliminate boundary effects due to limitations in log resolution Zhou and Esterle (2008). Geophysical log data is highly prone to errors and, hence, we defined theoretical rock and tool limits (Table 1) in order to remove erroneous data due to logging environments, incorrect tool calibration and data processing error. Compensated density (CODE) is one of the main input parameters for our modelling workflow, however, is not available in some of the study boreholes. Hence, we calculated this parameter using Eq. (1) (modified from Samworth, (1992)). In the next step, the 1 cm resolution log data were averaged across the sample interval (excluding null values). Finally, data without ash, relative density, volatile matter compensated density and gamma ray were removed. After the data preparation step, we have 1230 coal/sample data with geochemical analyses (i.e. proximate analyses and relative density) from 263 wells (Table 2). All coal samples had available caliper, short spaced, long spaced, compensated and gamma ray geophysical log data, and very limited VECTAR processed, resistivity and porosity logs. Therefore, we rely on these most common geophysical logs for most of our modelling.

CODE =

4 1 DENL − DENB 3 3

D = ⎛Dmid − ⎝ ⎜

TKig

⎞ + TK co 2 ⎠ 2 ⎟

(2)

Where D is the distance to intrusion, Dmid is the distance of the sample mid-point, TKig and TKsample are spatial mid-point coordinates and thickness for each intrusion interval and sample, respectively. From Fig. 4, it is observed that approaching intrusion, relative density, ash, total sulphur, moisture, geophysical log density and compression velocity show increasing trends, while volatile matter, fixed carbon, specific energy, geophysical log porosity, and geophysical log resistivity decrease. Phosphorous and geophysical log gamma ray do not appear to show significant increasing or decreasing trend. Transition from non-altered to fully altered occurs at ~1.5 m (around two times the thickness of intrusion) which is consistent with findings by Dow (1977). The increase in moisture is in line with Finkelman et al. (1998), and Goodarzi and Cameron, (1990) and decrease in resistivity and increase in geophysical log density is consistent with (Zhou and Esterle, 2008). In making these observations we assume that of all the geophysical log parameters, density, resistivity and porosity should contribute most to discriminating altered coal.

(1)

Where CODE is the compensated density, DENL is long spaced density and DENB is short spaced density. 4.4. Calculate sample distance to intrusion

4.5. Classification of each sample to altered or non-altered (Label altered coal using ash density relationship)

Calculating a samples distance to intrusion is considered important for two purposes, to incorporate as a parameter in modelling, and to determine the response of coal quality and geophysical log parameters as they approach intrusion. In our process, we use the intrusion

As discussed previously, the existing method for geochemically discriminating altered coal at the study mine is by using a raw volatile matter (ad) cut-off value of 8%. However, as seen in Fig. 5, this method may not always be accurate because high ash material (e.g. rock) has similar volatile matter to altered coal. Therefore, this method cannot be used across all ash ranges. Therefore, we prefer to use both ash and density to determine altered and non-altered material which can be used across all ash ranges. This can be achieved since unaltered material follows a linear ash relationship with the reciprocal of density and altered material does not. To achieve this differentiation, a linear separator (hyper-plane) is simply drawn through the minimum density ash relationship for normal, unaltered material (Fig. 6).

Table 1 Upper and lower limits for geophysical predictors based off theoretical limits of common strata tool detection and well diameter. Data outside of these limits were removed as erroneous. Tool code

Units

Description

Lower limit

Upper limit

CADE CODE DENB DENL DPRS FE1 FE2 GRDE HDEN

MM G/C3 G/C3 G/C3 PERC OHMM OHMM API G/C3

Caliper Density (Compensated) Density (Short Spaced) Density (Long Spaced) Density Porosity (Sandstone) Resistivity (Shallow) Resistivity (Deep) Gamma Ray (From Density) Density (Vector Processed)

30 1 1 1 0.01 0 0 0 1

110 5 5 5 200 50,000 50,000 500 5

5

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Table 2 Number of wells sample and log data available for analysis after data preparation (removal of outliers and aggregation). The acronyms are explained in Table 1. Number of data Wells

Samples

CADE

CODE

DENB

DENL

DPRS

FE2

FE1

GRDE

HDEN

263

1230

1230

1230

1230

1230

95

64

61

1230

568

5. Modelling results and discussion

to estimate other coal properties such as ash, yield and volatile matter (Zhou and Esterle, 2008) and can be used to help discriminate between altered and non-altered coals (Roslin and Esterle, 2015). A common method to establish the density of coal is by measuring its relative density according to AS1038.21.1.1–2002 (Standards Australia

5.1. Prediction of relative density from wireline log data Coal density is an important physical parameter which can be used

Fig. 4. Change in geochemical and geophysical properties approaching intrusion at study area. Change from completely fresh to completely altered at approximately 1.5 m (vertical dashed line) which is approximately 2 times the average thickness of dyke. 6

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 5. Plot showing separation of material into altered and non-altered class based on an 8% volatile matter cut-point. A linear separating plane cannot be fit through this data exactly and therefore cannot be used as a separator for all ash ranges.

Fig. 6. Separation of material into altered and non-altered class based on ash density relationship for 1230 coal samples in the study mine. 7

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 7. Relative density against each one of short spaced (DENB) long spaced (DENL) and compensated (CODE) geophysical density. Poor linear correlation is observed across all variables, therefore linear modelling of relative density using one of these parameters alone will produce poor predictions.

distribution which should produce an unbiased evaluation. After these data partitions were established, each algorithm was trained using caliper, compensated density and gamma ray log data as inputs and evaluated by calculating R2 values for 5-cross folds on the training data set and on the test set. The model with the highest R2 on the test set and train set was deemed to be the best model. The best model was then trained using the full (100%) data set and predicted relative density was added as a new variable to the data set. As shown in Fig. 9, the RF algorithm had the highest R2 for the training and test data, followed by the GBM algorithm. Fig. 10 shows an abstraction of the best-fit model based on RF algorithm for prediction of relative density which has 50 trees and maximum depth of 20. The LM performed poorly on both the training and testing data sets. It was found that for the RF algorithm, R2 was similar across all training folds (Appendix 1) indicating that the algorithm did not overfit to the training data. Reduced performance (smaller R2) of all algorithms on the test set may be due to relatively small sample size. Since RF and GBM's perform markedly better than LM we assume that these algorithms can exploit nonlinear relationships between relative density, caliper, compensated density and gamma ray log data and/or that they can better account for variations in caliper which affect compensated density.

International, 2002). Zhou and Esterle (2008) have previously shown that when good quality geophysical log data are available, relative density (according to AS1038.21.1.1–2002) can be accurately predicted using simple linear regression from different types of geophysical density logs such as short spaced (DENB), long spaced (DENL), compensated (CODE) and VECTAR processed (HDEN) data. Therefore, in order to see if there are significant correlations between relative density and geophysical density logs, as suggested by Zhou and Esterle (2008), we plotted CODE, DENB and DENL against relative density (Fig. 7). We do not use HDEN in our case due to its limited availability in the data set (Table 2). In contrast to Zhou and Esterle (2008), there are poor linear correlations between these parameters and relative density, hence simple linear modelling will produce poor predictions. We assume that these poor correlations are due a combination of factors including poor borehole conditions (irregularities in the borehole wall), differing borehole diameters, and resolution mismatch between the density tool and sampling intervals. In attempt to mitigate these factors and improve the prediction of relative density, we use compensated density (CODE), caliper and gamma ray log as input variables into a multi-linear model (LM). In addition, we also test two machine learning methods, including Random Forest (RF), Gradient Boosted Machine (GBM) to determine if these more complex non-linear methods improve over linear modelling. In our modelling procedure, the data set was randomly partitioned into an 80% training and a 20% test set and the distributions of the data were checked to ensure they have similar distribution (Fig. 8). This was done so that evaluation of models on the test set is from same/similar

5.2. Prediction of altered materials In this step, the aim was to predict altered and non-altered material by using a binary classification algorithm. Since altered coal cannot be

Fig. 8. Distribution of all data for each of the test (20%) and train (80%) data partitions plotted on the reciprocal of relative density against ash chart. These plots clearly show that there are a similar pattern/distribution between test and train dataset. 8

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 9. Regression results for the training and testing data set for the random forest (RF), gradient boosted machine (GBM) and multi-linear model (LM) algorithms using compensated density (CODE), caliper (CADE) and gamma (GRDE) as predictors. For the training data RF has the highest R2 (0.94) followed by GBM (0.87). LM is worst performing model with an R2 of 0.5. For the testing data RF and GBM have the highest R2 (0.76) followed by LM (0.66).

set of resistivity, these were included so that where available these could be used. Conversely, although porosity showed a clear relationship with distance to intrusion it was not used because it derived from the density log. To predict the alteration class, GBM, RF, and DNN classification algorithm were compared. Each algorithm has a unique

classified purely based off relative density, additional variables are needed for accurate classification. As such, in addition to the predicted relative density, the calculated distance to intrusion, gamma ray (GRDE), and short and long spaced laterlog resistivity data (FE1 and FE2) were included as input variables. Although there is a limited data

Fig. 10. Abstraction of the optimal random forest architecture for predicting relative density. The model uses caliper (CADE), Compensated density (CODE) gamma ray log (GRDE) as inputs. The number of trees is 50 and the number of variables randomly sampled is 1. 9

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

have both high precision and recall.

Table 3 Hyperparameter search ranges and results. Model

Hyper parameter

Search range

Optimal value

GBM

Maximum tree depth Number of trees Learning rate Activation function Hidden layers

5—20 ≤5000 ≤0.05 Rectifier, tanh, maxout (20,20), (50,50), (30,30,30), (25,25,25,25) 0.0000—1 × 10−4 0.0000—1 × 10−4 1–5

5 50 0.1 Maxout (20,20) 4.6 × 10−5 1 × 10−5 2

≤5000 ≤5000

20 50

DNN

RF

L1 regularisation L2 regularisation Number of variables randomly sampled Maximum tree depth Number of trees

Precision = Recall =

F1 = 2.

tp tp + fp

tp tp + fn

precision. recall precision + recall

(3) (4)

(5)

Where tp is true negative, and fp is false negative. In our case, a true positive is when altered class is predicted correctly whereas a false positive is when the class is predicted as altered but is non-altered. Similarly, a true negative is when unaltered class is predicted correctly, and a false negative is when a class is predicted as non-altered but is altered. Given these measures, the best model was chosen to be the one with lowest F1 across the training and test set. Finally, after the best model was determined it was retrained across the full data set and predictions were added to the data set. A plot of the locations of any misclassified results with respect to ash and density was also constructed to visualise areas of misclassification (Fig. 11). Random grid search for each of the algorithms resulted in 156 GBM models, 1000 DNN models and 137 RF models (Appendix 2). Using the best model for each algorithm we found that RF, as compared to DNN and GBM, had the highest precision, recall and F1 measure across both classes for each of the train and test set Table 5). We also found that these evaluation measures were similar for each of the cross-validation folds (Appendix 3). For the training data set, RF precision was 0.981 which matched that of the test data set. The RF was also able to perfectly classify non-altered material on the training data set and had high (0.97) precision on the test data set. In all, for RF there were only 11 misclassifications across both the training and testing data set of 1230 samples, representing error of < 1%. Addition to this, RF misclassifications were mostly very close to the class separating hyperplane except for few false negatives in the test set (Fig. 11). It should be noted that the best-fit RF model for prediction of altered materials has 50 trees, with maximum depth of 20, and the number of variables randomly sample is 2. The GBM was the next best classifier with precision of 0.91 for both the test and training data set for altered coal class. DNN was the worst performing algorithm with precision of 0.86 and 0.88 for the test and train data sets respective. Since the maximum number of models was limited to 1000 and the DNN reached this maximum, allowing the number of models to increase beyond 1000 may produce better result. Likewise, for the GBM algorithm, reducing the early stopping criteria a lower tolerance may have resulted in more modes being built resulting in better performance. However, given that RF produced the best results in the shortest time, and that this algorithm as the simplest and easiest to tune, attempting to further improve the DNN and GBM performance need not be considered.

set of hyperparameters and appropriate selection can result in increased classification performance (Horrocks et al., 2015; Xie et al., 2018). These hyperparameters may be selected based on a prior knowledge (from previous similar studies), by trial and error or by grid search which is the process of building models from several combinations of hyperparameters in a certain search range. A full grid search is one which builds modes over all combination of a given range of hyperparameters and can be computationally intensive and time consuming. Random grid search is the process of randomly selecting hyperparameter combinations from a very broad search space and has shown to be an effective way to reduce computation time, and can find combinations which yields results similar to or better than those from a restricted full grid search (Bergstra and Bengio, 2012). As such, we used random grid search on the previous partition of 80% of the data and used a broad hyperparameter search space (Table 3) to create numerous models for each algorithm. For each algorithm, the hyperparameter search was stopped when either there were 1000 individual models built, when run time reached 3600 s, or when five consecutive models had error within 0.1% of each other. For each algorithm, five-fold cross validation was conducted to evaluate performance on the training data set and to determine if algorithm is likely to under or overfit. This was done by comparing the evaluation metrics (Eqs. (3)–(5)) for each validation fold. Algorithms were considered stable (not overfitting or underfitting) when these evaluation metrics were similar. Each model was also evaluated by making predictions on the training (20% partition) in order to determine how the algorithm might perform on completely unseen data sets. For each algorithm and for each of the training and testing data sets, a confusion matrix was constructed and evaluation metrics: precision, recall and F1 were calculated in order to evaluate the performance per class. A confusion matrix allows us to easily visualise the predicted class against the actual class. An example of how a confusion matrix should be interpreted in the context of the study is provided in Table 4. Precision is defined as the number of true positives divided by the number of true positives plus the number of false positives (Eq. (3)). Recall is defined as the number of true positives divided by the number of true positives plus the number of false negatives (Eq. (4)). The F1 measure reflects both precision and recall (Eq. (5)) and a good classifier should

5.3. Assign lithology groups The predicted class from the RF algorithm (determined as the best model in 5.2) is divided into interpreted lithology groups including altered coal, fresh coal, carbonaceous rock, rock and altered coal-rock mix based on the predicted density and gamma ray log (GRDE). This

Table 4 Example of a confusion matrix in the context of the study. For each class (altered or non-altered) the number of predictions are tallied and compared against the actual class. Actual class

Predicted Class

Altered Non-Altered

Altered

Non-Altered

Number of altered class predicted correctly (true positive) Number of results predicted as non-altered but were altered (false negative)

Number of results predicted as altered but were non-altered (false positive) Number results which were precited as non-altered and were non-altered (true negative)

10

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 11. Location of mis-classified material for the RF algorithm with respect to ash and relative density for each of the train and test set. For the training data set, the algorithm predicted four results as altered, but these were non-altered (false positives). For the test data set the algorithm predicted six points as non-altered but these were actually altered (false negatives), and one result as altered which was actually nonaltered (false positive).

matter and geophysical log resistivity markedly decrease. These changes materially impact the beneficiation characteristics and product specification of the coal in the Leichhardt Seam in the studied mine. Predicting the impacted coal from lab derived relative density or wireline log density alone is difficult because the density ranges of altered coal are like high ash coal or carbonaceous material. In addition, the existing techniques to classify altered material based on volatile matter alone are not robust because they do not account for all ash ranges. As such we proposed an alternative method to classify altered and non-altered material based on ash density relationship which is more robust than existing volatile matter cut off method. We also showed that the altered and non-altered material can be predicted with a high degree of accuracy using machine learning algorithms from wireline log data in combination with calculated sample distance from intrusion and predicted relative density. We found that of machine learning algorithms random forests algorithm produced the best result with only 11 misclassifications across the entire data set of 1230 samples, representing a very low (< 1%) error. We showed that after making these predictions, we could accurately separate samples into broad lithology groups which very well match expected ash density ranges representative of true lithology groups. These high accuracy classifications will enable more robust resource definition of the case site by differentiating altered coal material from fresh coal material. The methodology also provides a semi-quantitative approach to classification of altered coal and provides an alternative to the existing manual qualitative approach currently undertake at the case site.

was done based on a rule-based system (Fig. 12) according to a prior knowledge of common lithology groups in the deposit area. In order to visually validate the results, these predicted lithology groups were plotted on the reciprocal density vs. ash % chart (Fig. 13). As a rule of thumb, fresh coal, carbonaceous material and rock should sit on the non-altered side of the ash/density hyper-plan and should have an ash range of < 25%, 25–50% and > 50%, respectively. Similarly, altered coal and altered coal rock-mix should sit on the altered side of the ash/ density hyper-plan and should have an ash range of < 50%. We found that predicted fresh coal, carbonaceous and rock groups broadly sit in their expected density and ash ranges (Fig. 13). Pure coal is very well predicted with only one sample plotting in the altered coal range and only minor confusion between coal and carbonaceous material occurs close to the 25% ash range. Rock samples are also well predicted with only two samples falling outside the expected the nonrock range (> 50%). Altered coal and altered coal/rock material are very well predicted with all sitting on the correct side of the ash/density hyperplane, and all but one sample in the expected (< 50%) ash range. While the classification of these lithology groups is based on a simple decision tree using a prior knowledge of the deposit, an alternative automated methodology such as clustering could be used to separate groups. 6. Conclusions In this study, we made the data driven observations that distance from intrusion has marked impact on coal quality and petrophysical qualities. In particular, relative density and geophysical log density markedly increase near intrusion, ash slightly increases, and volatile

Table 5 Confusion matrix results, precision, recall and F1 results for each one of Random Forest (RF), Gradient Boosted Machine (GBM) and Deep Neural Networks (DNN) for the training and test datasets. Note that the results are ordered by the best performing. Actual class Training data

Predicted Class

Testing data

Altered

Non-altered

Precision

Recall

F1

Altered

Non-altered

Precision

Recall

F1

RF

Altered Non-altered

207 0

4 760

0.98 1.00

1.00 0.99

0.99 1.00

Altered Non-altered

41 6

1 211

0.98 0.97

0.87 1.00

0.92 0.98

GBM

Altered Non-altered

205 2

20 744

0.91 1.00

0.99 0.97

0.95 0.99

Altered Non-altered

40 7

4 208

0.91 0.97

0.85 0.98

0.88 0.97

DNN

Altered Non-altered

174 33

28 736

0.86 0.96

0.84 0.96

0.85 0.96

Altered Non-altered

36 11

5 207

0.88 0.95

0.77 0.98

0.82 0.96

11

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Fig. 12. Decision tree used to split predicted altered and non-altered classes into lithology groups.

Fig. 13. Predicted lithology groups and their position on the reciprocal density ash plot.

Declaration of Competing Interest

study. Appendix A. Supplementary data

None.

Supplementary data to this article can be found online at https:// doi.org/10.1016/j.coal.2019.103284.

Acknowledgements The authors gratefully thank Dr. Karacan, Editor-in-Chief, and two anonymous reviewers of this article for their time and constructive comments. We also would like to thank Peabody for providing data for this study and for granting permission for these data to be used in this

References Al-Anazi, A., Gates, I.D., 2010. A support vector machine algorithm to classify lithofacies

12

International Journal of Coal Geology 214 (2019) 103284

K. Maxwell, et al.

Illinois Basin. Int. J. Coal Geol. 165, 243–256. Quaderer, A., Mastalerz, M., Schimmelmann, A., Drobniak, A., Bish, D.L., Wintsch, R.P., 2016. Dike-induced thermal alteration of the Springfield Coal Member (Pennsylvanian) and adjacent clastic rocks, Illinois Basin, USA. Int. J. Coal Geol. 166, 108–117. R Core Team, 2017. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria. Rimmer, S.M., Yoksoulian, L.E., Hower, J.C., 2009. Anatomy of an intruded coal, I: Effect of contact metamorphism on whole-coal geochemistry, Springfield (No. 5) (Pennsylvanian) coal, Illinois Basin. Int. J. Coal Geol. 79, 74–82. Ritchie, C., 2010. Lamprophyric Intrusions in the Rangal Coal Measures, Bowen Basin Classification, Geochemistry and Tectonic Significance. School of Earth Science/ University of Queensland. Roslin, A., Esterle, J.S., 2015. Electrofacies analysis using high-resolution wireline geophysical data as a proxy for inertinite-rich coal distribution in Late Permian Coal Seams, Bowen Basin. Int. J. Coal Geol. 152, 10–18. Roslin, A., Esterle, J.S., 2016. Electrofacies analysis for coal lithotype profiling based on high-resolution wireline log data. Comput. Geosci. 91, 1–10. Salmachi, A., Rajabi, M., Reynolds, P., Yarmohammadtooski, Z., Wainman, C., 2016. The effect of magmatic intrusions on coalbed methane reservoir characteristics: A case study from the Hoskissons coalbed, Gunnedah Basin, Australia. Int. J. Coal Geol. 165, 278–289. Samworth, J., 1992. The dual-spaced density log-characteristics, calibration, and compensation. Log Anal. 33. Schmidhuber, J., 2015. Deep learning in neural networks: An overview. Neural Netw. 61, 85–117. Schmitt, P., Veronez, M.R., Tognoli, F.M.W., Todt, V., Lopes, R.C., Silva, C.A.U., 2012. Electrofacies modelling and lithological classification of coals and mud-bearing finegrained siliciclastic rocks based on neural networks. Earth Science Research 2. Shi, Q., Qin, B., Bi, Q., Qu, B., 2018. An experimental study on the effect of igneous intrusions on chemical structure and combustion characteristics of coal in Daxing Mine, China. Fuel 226, 307–315. Society of Exploration Geophysicists, 2016. 2016-ml-contest. Standards Australia International, 2002. AS 1038.21.1.1–2002, Coal and Coke – Analysis and Testing Part 21.1.1: Higher Rank Coal and Coke – Relative Density – Analysis Sample/Density Bottle Method. Standards Australia International. Svozil, D., Kvasnicka, V., Pospichal, J.í., 1997. Introduction to multi-layer feed-forward neural networks. Chemom. Intell. Lab. Syst. 39, 43–62. Tan, P.-N., Steinbach, M., Karpatne, A., Kumar, V., 2018. Introduction to Data Mining, 2nd edition. Pearson, pp. 864. Team, RStudio, 2016. RStudio: Integrated Development Environment for R. RStudio, Inc, Boston, MA. The H2O.ai team, 2018. h2o: R Interface for H2O. The H2O.ai team, 2019. Welcome to H2O 3. Thomas, L., 2012. Coal Geology. Wiley-Blackwell, pp. 454. Wu, Q., Wu, Q., Xue, Y., Kong, P., Gong, B., 2018. Analysis of overlying strata movement and disaster-causing effects of coal mining face under the action of hard thick magmatic rock. Processes 6, 150. Xie, Y., Zhu, C., Zhou, W., Li, Z., Liu, X., Tu, M., 2018. Evaluation of machine learning methods for formation lithology identification: A comparison of tuning processes and model performances. J. Pet. Sci. Eng. 160, 182–193. Yao, Y., Liu, D., 2012. Effects of igneous intrusions on coal petrology, pore-fracture and coalbed methane characteristics in Hongyang, Handan and Huaibei coalfields, North China. Int. J. Coal Geol. 96-97, 72–81. Zhang, L., Zhan, C., 2017. Machine learning in rock facies classification: an application of XGBoost. In: International Geophysical Conference, Qingdao, China, 17-20 April 2017. Society of Exploration Geophysicists and Chinese Petroleum Society, pp. 1371–1374. Zhou, B., Esterle, J., 2008. Toward improved coal density estimation from geophysical logs. Explor. Geophys. 39. Zhou, B., O'Brien, G., 2016. Improving coal quality estimation through multiple geophysical log analysis. Int. J. Coal Geol. 167, 75–92.

and model permeability in heterogeneous reservoirs. Eng. Geol. 114, 267–277. Bergstra, J., Bengio, Y., 2012. Random search for hyper-parameter optimization. J. Mach. Learn. Res. 13, 281–305. Bhattacharya, S., Mishra, S., 2018. Applications of machine learning for facies and fracture prediction using Bayesian Network Theory and Random Forest: Case studies from the Appalachian basin, USA. J. Pet. Sci. Eng. 170, 1005–1017. Breiman, L., 1996. Bagging predictors. Mach. Learn. 24, 123–140. Breiman, L., 2001. Random forests. Mach. Learn. 45, 5–32. Breiman, L., 2017. Classification and Regression Trees. Routledge. Candel, A., LeDell, E., 2019. Deep Learning with H2O. H2O.ai. Chen, J., Liu, G., Li, H., Wu, B., 2014. Mineralogical and geochemical responses of coal to igneous intrusion in the Pansan Coal Mine of the Huainan coalfield, Anhui, China. Int. J. Coal Geol. 124, 11–35. Cooper, J.R., Crelling, J.C., Rimmer, S.M., Whittington, A.G., 2007. Coal metamorphism by igneous intrusion in the Raton Basin, CO and NM: Implications for generation of volatiles. Int. J. Coal Geol. 71, 15–27. Corina, A.N., Hovda, S., 2018. Automatic lithology prediction from well logging using kernel density estimation. J. Pet. Sci. Eng. 170, 664–674. Dietterich, T.G., 2000. Ensemble methods in machine learning. In: International Workshop on Multiple Classifier Systems. Springer, pp. 1–15. Dubois, M.K., Bohling, G.C., Chakrabarti, S., 2007. Comparison of four approaches to a rock facies classification problem. Comput. Geosci. 33, 599–617. Elseberg, J., Magnenat, S., Siegwart, R., Nüchter, A., 2012. Comparison of nearestneighbor-search strategies and implementations for efficient shape registration. J. Software Eng. Robotics 3, 2–12. El Sharawy, M.S., Nabawy, B.S., 2016. Determination of electrofacies using wireline logs based on multivariate statistical analysis for the Kareem Formation, Gulf of Suez, Egypt. Environ. Earth Sci. 75. Finkelman, R.B., Bostick, N.H., Dulong, F.T., Senftle, F.E., Thorpe, A.N., 1998. Influence of an igneous intrusion on the inorganic geochemistry of a bituminous coal from Pitkin County, Colorado. Int. J. Coal Geol. 36, 223–241. Firth, D., Elkington, P., 1999. Log Analysis for Mining Applications. Reeves Wireline Services. Fredericks, P.M., Warbrooke, P., Wilson, M.A., 1985. A study of the effect of igneous intrusions on the structure of an Australian high voltage bituminous coal. Org. Geochem. 8, 329–340. Friedman, J.H., 2002. Stochastic gradient boosting. Comput. Stat. Data Anal. 38, 367–378. Friedman, J., Hastie, T., Tibshirani, R., 2001. The Elements of Statistical Learning. Springer Series in Statistics New York. Goodarzi, F., Cameron, A.R., 1990. Organic petrology and elemental distribution in thermally altered coals from Telkwa, British Columbia. Energy Sources 12, 315–343. Golab, A.N., Hutton, A.C., French, D., 2007. Petrography, carbonate mineralogy and geochemistry of thermally altered coal in Permian coal measures, Hunter Valley, Australia. Int. J. Coal Geol. 70, 150–165. Green, D.R., 2015. CoalLog Geology and Geotechnical Training Manual. ACARP. Gu, Y., Bao, Z., Rui, Z., 2018. Complex lithofacies identification using improved probabilistic neural networks. Petrophysics 59, 245–267. Horrocks, T., Holden, E.-J., Wedge, D., 2015. Evaluation of automated lithology classification architectures using highly-sampled wireline logs for coal exploration. Comput. Geosci. 83, 209–218. Jiang, J., Cheng, Y., 2014. Effects of igneous intrusion on microporosity and gas adsorption capacity of coals in the Haizi Mine, China. Sci. World J. 2014. Johnson, V., Gray, R., Schapiro, N., 1963. Effect of igneous intrusives on the chemical, physical, and optical properties of Somerset coal. Am Chem Soc Div Fuel Chem 7, 110–124. Kuroda, M.C., Vidal, A.C., Leite, E.P., Drummond, R.D., 2012. Electrofacies characterization using self-organizing maps. Revista Brasileira de Geofísica 30. Lu, C.-P., Liu, Y., Wang, H.-Y., Liu, P.-F., 2016. Microseismic signals of double-layer hard and thick igneous strata separation and fracturing. Int. J. Coal Geol. 160-161, 28–41. Presswood, S.M., Rimmer, S.M., Anderson, K.B., Filiberto, J., 2016. Geochemical and petrographic alteration of rapidly heated coals from the Herrin (No. 6) Coal Seam,

13