Unsupervised machine learning for data-driven classification of rock mass using drilling data
How can a data-driven system handle limitations in existing rock mass classification systems?

RMCSs have been essential for the rock engineering industry and have contributed to increased consistency in rock mass assessment and rock support worldwide for decades. Developed mainly in the 1970s, these systems predate modern data capture technologies such as comprehensive scan/image capturing of newly exposed rock surface profiles, Measure While Drilling (MWD) data, and geophysics from the excavation face. At that time, high-resolution datasets with extensive rock mass coverage and advanced statistical learning techniques were unavailable, and computational power was limited. Given the current advances in data availability and automation potential in Scandinavian tunnelling, these classification systems exhibit limitations that may result in suboptimal decisions. Modern data-driven approaches could address these limitations more effectively. The objective here is not to undermine the existing systems but to highlight the enhancements possible with today’s technology.

1. 1.

   Subjectivity in Assessment: These systems rely heavily on the subjective judgement of face engineers, leading to variations in the assessment and support decisions for identical rock conditions [4, 5, 6].

2. 2.

   Inconsistent Observations: Face engineers may focus inconsistently on specific features or areas, potentially overlooking variations in the rock mass, and might not be able to perceive the exposed rock in the last blasting round similarly [6, 4]

3. 3.

   Safety and Accessibility Limitations: High-risk conditions and physical barriers often prevent thorough inspection of exposed rock, particularly in large double-track tunnels. Several of these tunnels are inspected from the floor at a safe location, leading to a poor assessment. [7, 6]

4. 4.

   Constraints on Quantification: Quantifying representative values for the newly exposed rock mass, such as the Rock Quality Designation (RQD), within the limited time frame of a fast-paced tunnel cycle presents significant challenges. The RQD metric, which involves measuring core pieces longer than 10 cm for every meter of rock, is tough to assess accurately under these conditions. [8, 9].

5. 5.

   Conservative Over-Supporting: Typically, the poorest rock mass conditions dictate the support for the entire blasted area, which can lead to unnecessary over-support in better-quality sections [7, 6]

6. 6.

   Empirical Data Limitations: The original empirical data may not cover all geological conditions, construction geometries, or site-specific factors, resulting in an oversimplified approach to complex geologies [7, 6, 9]

7. 7.

   System Update Challenges: Updates to these systems are infrequent and labour-intensive, hindered by non-transparent processes and inherent biases. The process involves trial and error to adapt input configurations to experienced stability in existing sites [6, 7, 5].

8. 8.

   Complex Rules for Exceptions: The systems incorporate complicated rules and factors that adjust classifications for different scenarios, which can lead to errors if not applied correctly, E.g. forgetting to multiply the Jn value (number of joint sets) in a junction with three, you might end up with a rock support class that is too low [7, 6].

9. 9.

   Visual Assessment Limitations: Current systems focus on visually assessable rock mass, neglecting the stability of rock outside the immediate tunnel profile, which is crucial for overall stability [7, 6].

10. 10.

    Non-existing advance rock mass assessment: The existing system cannot effectively assess the rock mass quality in front of the excavation, making it less useful for decision support on advance support and excavation method. You might say you can classify a drilled rock core from the face in RMR or Q-class, but such a process severely impacts the efficient tunnel factory and is only a point value. A system which describes advance support classes (face bolts, stability grouting, spiling bolts, etc.) from data ahead of the tunnel remains elusive [6, 7].s

11. 11.

    Concervative Support: Rock support classes are defined by inspecting primarily stable conditions, majorly in civil infrastructure tunnels with a high safety factor and where the rock support is conservative, which may not reflect the actual stability needs [7, 6].

12. 12.

    Mismatch in support classes: There are criticisms regarding how the combination of defined classes and linked rock support description can address the right type of rock support in general. [10, 7, 11].

13. 13.

    Failure Modes: Existing systems may not adequately consider various failure modes in their classifications, impacting the accuracy of rock support assessments [7, 6]

Existing rock mass classification systems face several limitations: they are challenging to update, inherently conservative, subject to user bias, lack sufficient details, are unsuitable for forecasting, and do not assess rock mass where support is most needed. Additionally, the conservative nature of the industry often places higher trust in the subjective assessments of experienced human experts over automated systems [6, 12]. This resistance to automation may hinder updates that could incorporate more systematic data collection. Transitioning some or all rock mass assessments to a transparent, reproducible, bias-free, and easily updatable data-driven decision support system could address these issues effectively.

An accurate and objective understanding of rock mass stability is crucial to optimise rock support, blasting design, and excavation methods in tunnelling and mining. Rock materials are inherently heterogeneous, and the quality of the rock mass can vary significantly over short distances. However, historically, we have not been able to comprehensively describe this complexity to facilitate optimised decision-making in the fast-paced tunnel cycle. Consequently, rock mass quality is grouped into practical and meaningful target classes.

Existing classification systems are impractical in environments that are inaccessible to humans. These include hazardous areas of weak rock in current tunnelling projects, production mining environments operating close to a safety factor 1.0, and extraterrestrial locations planned for future human bases. In these scenarios, automated systems are necessary for thorough rock mass assessments.

A finely tuned rock mass classification on a sufficiently small scale is essential for optimising decisions. Current broad-spectrum decisions are predominantly conservative, leading to the excessive use of steel and concrete for rock support. Moreover, there is a significant gap in our ability to assess rock mass quality ahead of the excavation, which is vital for planning advance support and selecting the appropriate excavation method. Addressing these challenges and the limitations outlined in Section 1.1 requires the development of new, more adaptive systems that can operate autonomously and provide accurate assessments in real time to enhance safety and efficiency in challenging and dynamic environments.

To date, no purely data-driven rock mass classification system for tunnelling and underground mining encompasses the following properties: (a) the capability to classify rock mass stability of the visually exposed rock mass and outside the tunnel profile, including the area ahead of the tunnel face; (b) independence from existing classification systems; (c) the use of comprehensive data, patterns and decisions from larger rock volumes than individual drillholes, such as entire blasting rounds or a 1 m slice of the tunnel; (d) practical classification of rock mass into stability classes that aid decision-making, beyond merely correlating mechanical properties like UCS, E-modulus, and single features.

The necessity for criterion (a) stems from assessing the rock mass’s stability surrounding and ahead of the tunnel face. Criterion (b) avoids inheriting limitations from previous systems. Criterion (c) acknowledges that the rock mass quality can vary significantly over short distances. It is the combined signature of the rock mass for a bigger volume that should be used to make practical decisions regarding rock support and excavation methods. Classifying from single drillholes might introduce so much noise and variations when making relevant decisions, based on the forecasting information from all the drillholes, that an extra interpretation step is needed to make a decision such as "should I install spiling bolts or not". Criterion (d) highlights that decisions on rock support and excavation methods cannot rely solely on single mechanical properties or features.

To our knowledge, the datasets derived from face seismics and Measure While Drilling (MWD) in hard rock tunnelling are the only ones providing the necessary spatial detail to effectively act as signatures of the rock mass beyond the tunnel profile. Several studies have linked high-resolution datasets to rock mass quality without explicitly aiming to develop a new data-driven classification system. We have categorised these studies into two groups based on their use of these datasets for forecasting purposes. Using geophysical data, Dickmann et al. [13] and Dickmann and Hecht-Méndez [14] have automated the characterisation of rock mass into stability zones ahead of the tunnel and linked rock support to ground treatment using tunnel seismic prediction (TSP). Sapronova et al. [15] employed unsupervised clustering to group similar data into clusters representing different geological conditions, subsequently labelled through supervised learning. These studies show promising results in characterising the rock mass ahead of the tunnel face. However, the development of a complete rock mass classification system is pending, and accuracy remains to be enhanced. Moreover, the use of geophysics at the excavation face involves human intervention. It could significantly affect the drill and blast cycle more than automatically collected MWD data.

Sapronova et al. [16] utilises correlation values between MWD features from single drillholes as feature vectors to predict Q-system classes, which contradicts points b and c regarding independence from existing systems and analysing larger rock volumes. Similarly, Hansen et al. [17] employs statistically derived MWD values from all drillholes in a blasting round to predict Q-classes, yet still breaches point b by relying on existing classification systems. Fernández et al. [18] applies machine learning to single-hole MWD data to detect discontinuities using a calculated discontinuity index, thus contravening point c about using summarised data decisions from larger rock volumes. van Eldert et al. [19] predicts Q-values and rock support from single-hole MWD data, violating points b and c, although the study’s calculation and visualisation of a fracture index from single holes, segmented manually into different rock mass quality zones, moves towards a purely data-driven system. However, segmentation must be automated to reduce human bias, and explicit support classes must be linked to fracture indices. He et al. [20] also uses single-hole MWD data to predict UCS and friction angles with machine learning, infringing upon points c and d, which call for broader data integration and practical decision-making utility. Lastly, Galende-Hernández et al. [21] clusters MWD data and links these to RMR values using expert-based fuzzy rules, breaching point b.

We can think of two natural strategies for optimising rock mass classification systems in a data-driven way: (a) employing modern data collection and learning techniques to improve the existing systems in various ways (e.g. improve assessment of input variables or extending the systems with new features), or (b) developing entirely new, purely data-driven systems. This study adopts the latter approach, hypothesising that rock mass can be intricately grouped, clustered, and classified using spatially extensive, high-resolution Measure While Drilling (MWD) data, which serves as a signature of the rock mass. Our objective is to address several limitations in existing systems in this context. MWD data is a cost-effective and easily retrievable data source in global tunnelling and mining operations [22], with the added advantage of not impacting the tunnel cycle. We analysed the natural clustering of MWD data, organised as tabular data samples for every meter of tunnel excavation across 15 hard rock tunnels, totalling 23,000 meters and involving approximately 500,000 blasting drillholes of infrastructure tunnel data:

• •

  Information extraction involved calculating six statistical features from about 5000 values for each of eight MWD parameters, yielding 48 values in total.

• •

  To enhance clustering, we reduced the dimensionality of the 48-value vectors using nonlinear manifold learning techniques such as UMAP and linear Principal Component Analysis (PCA).

• •

  We employed unsupervised machine learning techniques (HDBSCAN, Agglomerative Clustering, K-means) to identify natural groupings and structures within the data, creating clusters. We further optimised and explored the hyperparameters in these algorithms using multi-objective optimisation to ensure effective and meaningful clustering.

• •

  We mapped and structured the clusters according to the physical properties of the rock mass, such as rock type and quality, and the distributions of key MWD parameters essential for rock mass assessment and clustering to determine if the clustering meaningfully differentiated the rock mass.

The subsequent sections will outline the dataset, methods used, results, and their analysis. The methodology section covers dimension reduction, clustering, decision metrics, organising experiments, hyperparameter optimisation, and linking clusters to physical properties of the rock mass. The discussion section explores the implications of these results. Finally, the conclusion and outlook sections summarise our findings and suggest directions for future research.

Dimension reduction prior to clustering improves computational efficiency and clarity by reducing noise and redundancy, thereby enabling clearer and more significant groupings in the reduced feature space [25]. In addressing the complexities of rock mass categorisation, we have applied dimension reduction techniques. Specifically, we utilised UMAP (Uniform Manifold Approximation and Projection) [26] for its non-linear capabilities and PCA (Principal Component Analysis) [27] for linear reduction. These methods aim to preserve and highlight local and global data structures, facilitating pattern recognition and grouping by clustering algorithms. Common in fields like bioinformatics, this sequential approach of reduction followed by clustering accommodates heterogeneous data [28]. Aware of potential distortions in cluster relevance due to manifold learning [29], we also perform clustering directly on the MWD-features, and we ensure a detailed inspection of their distributions in all experiments (not the dimension-reduced components) and alignment with existing dataset labels. Additionally we employed dimension reduction with UMAP to effectively visualise the high dimensional data with clusters in 2D and 3D plots.

UMAP is particularly useful for clustering, as it can capture the non-linear relationships in the data, which can be important for rock mass classification. UMAP (Uniform Manifold Approximation and Projection) is a dimension reduction technique that can be used in exploratory data analysis and visualisation, as well as in machine learning. It is a non-linear technique that preserves the global and local structure of the data, making it particularly useful for clustering. UMAP works by constructing a high-dimensional graph representation of the data and then optimising a low-dimensional version of this graph to be as structurally similar as possible. This results in a low-dimensional data representation that maintains much of the original structure. Conversely, PCA is selected for its efficiency and simplicity in linear dimension reduction. It works by identifying the directions (principal components) that maximise the variance in the data, which, for a dataset with a vector length of 50, is an effective strategy for reducing dimensionality while retaining as much information as possible. PCA’s strength lies in its ability to provide a clear overview of the data’s linear structure, making it a suitable linear counterpart to UMAP in our study. Using both UMAP and PCA, we can explore the data from different perspectives and potentially improve the performance of the clustering algorithms.

While other dimension reduction techniques exist, both linear—such as Linear Discriminant Analysis (LDA) and Singular Value Decomposition (SVD)—and non-linear—such as t-SNE and autoencoders [27, 30, 31], we opted for UMAP and PCA due to their effectiveness, simplicity, and efficiency in handling the complexities of rock mass data. t-SNE, although powerful for manifold learning, is computationally intensive and challenging to tune and interpret. LDA depends on predefined class labels, which we aim to avoid as initial inputs. SVD is primarily utilised for matrix factorisation in contexts like natural language processing and autoencoders, despite their capability for non-linear dimension reduction, demanding intricate tuning, and being complex to implement. Therefore, UMAP and PCA were selected for their balanced performance, ease of use, and interpretability, which align well with the objectives of our study.

Following dimension reduction employing PCA, UMAP, and direct analysis without dimension reduction, we conducted clustering to understand the structural variations within the MWD signature of the rock mass. This approach is essential for identifying inherent groupings within the MWD data. We selected three distinct clustering algorithms based on their methodological diversity and applicability to our data type:

• •

  K-means Clustering: This algorithm was chosen for its efficiency in handling big data sets and simplicity, making it highly interpretable. K-means partitions the dataset into K distinct, non-overlapping clusters by minimising the variance within each cluster [32]. This method is particularly effective for identifying spherical clusters in feature space.

• •

  Agglomerative Clustering: As a hierarchical clustering technique, this algorithm was employed to provide insights into the possible hierarchical structure of the dataset [33]. It progressively merges pairs of clusters that minimally increase a given linkage distance. Agglomerative clustering is useful for our study as it allows the examination of cluster structures at different scales.

• •

  HDBSCAN: Selected as an advanced density-based algorithm, HDBSCAN extends DBSCAN by converting it into a hierarchical clustering algorithm and introducing a stability-based cluster selection technique. This choice is justified by its ability to handle variable cluster densities, which is crucial for datasets with complex spatial relationships like those found in geotechnical data [34]. HDBSCAN determines the core clusters based on the areas of high density and treats points in low-density regions as outliers. Unlike K-means, which forces every point into a cluster even if it does not logically belong to any, HDBSCAN allows points to remain unclassified if they do not fit well with any cluster (given by the cluster value -1). The ability of HDBSCAN to identify and separate outliers is beneficial in geotechnical data, where anomalous readings can indicate critical phenomena like unstable rock masses or unusual material properties. By effectively isolating these outliers, HDBSCAN can help focus attention on potential areas of concern that might require further investigation or monitoring, enhancing the safety and reliability of geotechnical assessments

The clustering approach does not rely on a predefined number of clusters; instead, it employs the Silhouette, Davies-Bouldin, and Calinski-Harabasz scores to optimise the cluster count. This method ensures that the optimal number of clusters is determined based on the data characteristics, avoiding bias from preset values. The implementation varies across algorithms; for instance, HDBSCAN does not require specifying a cluster number, whereas K-means and agglomerative clustering initially require a set number of clusters. However, as detailed later, the number of clusters is subsequently optimised, mitigating any potential issues from initial settings.

Four distinct feature sets were utilised to analyse natural clustering tendencies in MWD data. Each set is identified by its name at the beginning of the sets described below.

• •

  All. All 48 MWD-parameters, plus the two geometric parameters overburden thickness and tunnel width (50 in total), to investigate in what way extra feature information, in addition to the MWD rock mass signature, impacts the clustering result.

• •

  MWD. All 48 MWD parameters are used to inspect clustering results when only the MWD signature is used.

• •

  MWD_rock. Using domain and data knowledge to reduce the feature set to 30 by removing the 12 MWD-parameters dealing with WaterFlow, which might not impact the rock mass stability, and the correlated standard deviation parameter (correlated to variance).

• •

  MWD_median. The median of the 48 MWD parameters, giving eight features, is analysed to determine whether it adequately represents the stability of the rock mass.

The evaluation of clustering algorithms is a multifaceted process, presenting unique challenges compared to supervised learning. Unlike supervised learning, where precision and recall can directly measure performance, clustering evaluation must focus on the data’s intrinsic structure without relying on predefined labels. This involves assessing whether the algorithm effectively identifies meaningful groups within the data that reflect some form of ground truth or underlying assumptions about data similarity [24]. The evaluation incorporates established clustering scores, qualitative assessments, and visual inspections, demonstrating the thoroughness and complexity of the process.

The experimentation process in this study is structured as a pipeline comprising scaling, dimension reduction, clustering, and evaluation, implemented in Python and hosted on GitHub. Apart from UMAP [26], all components, including scalers, dimension reduction algorithms, clustering algorithms, and metrics, are sourced from Scikit-learn [24].

The study conducted over a thousand experiments, primarily focusing on hyperparameter optimisation, with approximately 20% undergoing detailed scrutiny. For efficient management and inspection of these experiments, we utilised mlflow [41] and hydra [42]. Mlflow was used to store and visualise the results of the experiments, along with the most important parameters and artefact files of features, true labels, and cluster values. Hydra was used to organise and store the detailed configurations of the experiments and to make it easy to run the experiments with different configurations. All configuration values were quality-checked using Pydantic [43] in each experiment. The experiments were run in a makefile, which is a simple way to run the experiments in an efficient and reproducible way [44]. Result plots and compilations were generated reproducibly by retrieving data from mlflow using unique experiment IDs. For further insights into scientific ML-based reproducible and transparent experimentation, see [45].

Optimising hyperparameters is essential for effective clustering, given the sensitivity of results to these settings in dimension reduction and clustering algorithms. Proper hyperparameter configuration significantly impacts clustering outcomes. As outlined in Section 3.3, evaluating multiple metrics is crucial in unsupervised clustering. In this study, the Silhouette score, Davies Bouldin index, and Calinski Harabasz index were selected as objective scores for optimisation. We employed Bayesian multi-objective optimisation, using the Optuna package [46], to adjust the hyperparameters in several stages.

For hyperparameter sampling, we adopted the multi-objective version of the Tree-structured Parzen Estimator (MOTPE) [47]. MOTPE effectively approximates a Pareto front and is more efficient than evolutionary samplers like NSGA-II, which directly optimise a Pareto front. A Pareto front represents the set of solutions considered non-dominated, meaning no solution can improve one objective without worsening at least one other objective [48]. This concept is fundamental in multi-objective optimisation as it visually and numerically illustrates the trade-offs between competing objectives, helping decision-makers choose the most suitable solutions according to their preferences. We provide details on the optimisation process due to its significance.

Approximating the Pareto Front: MOTPE extends the traditional TPE for multiobjective optimisation by defining vector dominance and using a nondomination ranking system. Observations are split into two groups via density functions based on their dominance relations to a reference set $Y^{*}$, ensuring $p(y\succ Y^{*}\cup y\parallel Y^{*})=\gamma$. Additionally, the Hypervolume Subset Selection Problem (HSSP) is utilised to select observations that maximise the hypervolume indicator, facilitating an effective and diverse approximation of the Pareto front.

Selection of Pareto Optimized Experiments: In MOTPE, the selection of experiments is guided by the Expected Hypervolume Improvement (EHVI) criterion, which prioritises candidate solutions based on their potential to extend the current Pareto front. During each optimisation iteration, candidates are sampled from a distribution modelled by $l(x)$, and the one with the highest EHVI is selected for evaluation. This greedy approach ensures continuous refinement of the Pareto front by balancing exploration and exploitation throughout optimisation. MOTPE builds models to estimate the probability distribution of objective values given the hyperparameters and then samples new hyperparameters from areas expected to yield better objective values. The four most important parameters were defined for each algorithm (HDBSCAN, K-means, UMAP, etc.), and ranges of values were set up for the process of choosing new parameters.

1. 1.

   Sampling. The MOTP sampler samples a new set of parameters from the hyperparameter space. This sampling is influenced by the past performance of parameter sets, aiming to explore regions with potentially better trade-offs between objectives.

2. 2.

   Evaluation. Using a parametric pipeline including scaling, dimension reduction and clustering, a new set of hyperparameters was tested in each iteration, and the three objective metrics were reported from each run. According to each metric’ optimised goal, the target metrics were set to be (ref. order above) maximised, minimised and maximised.

3. 3.

   Updating. Based on the outcomes of the evaluations, Optuna updates its understanding of the hyperparameter space. For TPE, it involves updating the internal probabilistic models.

4. 4.

   Selection: In the next iteration, Optuna uses the updated information to sample new parameters again, effectively iterating towards better solutions over time. In multi-objective optimisation, the ’better’ solution must consider the trade-offs between competing objectives.

In choosing the final best set of hyperparameters for each pipeline this study focused on the Silhouette score and the Calinski Harabasz index, when several Pareto optimal solutions were found.

This study primarily investigates the natural clustering of MWD-data to establish a foundation for a data-driven rock mass classification system. The second step involves examining the physical properties of the rock mass within each cluster. Although this research initiates the process and demonstrates promising results and guidelines, further investigation into the properties of each cluster and refinement of cluster compositions are necessary to enhance the system’s industrial applicability. That is a task for future studies. The properties of each cluster have been analysed in multiple ways across three high-performing experiments:

• •

  The cumulative distribution function (CDF) plots for three important MWD-parameters—Normalised penetration, feeder pressure, and rotation pressure (torque)—are presented for each cluster, as identified in studies by [49, 50]. These parameters, which intuitively reflect the physical properties of the rock mass (e.g., higher penetration suggests softer rock, and higher feeder pressure indicates better rock quality), are analysed to determine if clustering effectively discriminates between rock masses with distinct physical characteristics. A good spread of distinct distribution suggests effective separation, whereas overlapping distributions imply less distinction.

• •

  A table is provided that lists the median values and sample counts of the three MWD-parameters across clusters, sorted by increasing normalised penetration.

• •

  To align with existing classification systems, rock types and Q-class labels, determined by majority vote, are assigned to each cluster. This alignment, while not a primary objective, aids in integrating the new system with established classifications and enhancing its acceptance.

• •

  In one experiment, clusters are sequenced based on the typical compressive strength of the assigned rock types (rock types are well aligned with the clusters in this experiment).

Approximately 50 experiments were required for three cluster metrics to converge, enabling the identification of optimal hyperparameters for each pipeline. Optimisation proved essential for generating informative clusters as default algorithm settings were inadequate. The goal was to maximise the Silhouette and Calinski-Harabasz scores while minimising the Davis-Bouldain index. The best hyperparameters identified in each run were logged in mlflow and applied in the final clustering experiments. The results of this hyperparameter optimisation are displayed in an interactive Pareto plot using a Plotly function [40], implemented in the Optuna package [46]. The plot in Fig. 3 shows the trials for an ’MWD’ featureset optimisation involving UMAP and Agglomerative Clustering. Increasing trial numbers are indicated by progressively darker shades of blue, and improved Pareto optimal solutions are marked in darker reds. This presentation method effectively compares experiment outcomes against all objectives and illustrates the typical scatter development from the optimisation process. The MOTPE-sampler methodically explores parameter configurations, gradually converging on a Pareto front of optimal solutions.

The optimisation process, along with reporting and plotting, was executed for each distinct combination of feature set, dimension reduction technique, and clustering algorithm, leading to a comprehensive database of experiments tracked in mlflow. The four most effective experiments for a combination, as determined by the Pareto optimisation process, had their hyperparameters and unique IDs recorded. Each experiment was then thoroughly examined for its clustering results, and at least one experiment per configuration is detailed in two consolidated tables, Table 2 and Table 3 , in Section 4.2. The optimised hyperparameters for three of these experiments are listed in Table 7 in Appendix A. Any parameters not specified revert to the algorithm’s default settings.

Optimisation of hyperparameters significantly enhanced the clustering results, as evidenced by improved scores in the clustering metrics for experiments 0, 3, and 7, detailed in Table 2. The SC, DBI, and CHI metrics showed marked improvements, particularly in pipelines using the HDBSCAN algorithm. Through optimisation, HDBSCAN pipelines transitioned from 25% unclustered samples and over 800 clusters to effective clustering. In contrast, pipelines with the agglomerative clustering algorithm exhibited scores for default parameters that were closer to those of optimised parameters. Additionally, the occasional higher CHI scores for default configurations underscore the necessity of optimising multiple metrics and conducting a comprehensive evaluation of various approaches.

In Table 2 and 3, the results of the best-performing experiments, trained with optimised parameters, are presented. The results are grouped by featureset and then sorted after the Silhouette score.

In Table 4, 5, and 6, we present detailed results from experiments 0, 3, and 7. The clusters are arranged by increasing values of the Penetration MWD feature. It is important not to focus on the absolute values of the MWD features as they are scaled, but rather on their relative sizes.

We observe distinct patterns for experiment 0 using the ’All’ feature set. The number of samples in each cluster is evenly spread, correlating with a low Gini score. Each cluster’s median feature values are distinct, indicating differing physical properties. Had the median values been similar across multiple clusters, it would suggest poor clustering relative to rock mass properties. The overburden feature generally decreases as penetration increases, possibly due to reduced rock stress and more dayrock in areas with lower overburden and higher penetration. Cluster 0, with only 22 samples, appears to be an outlier. The rock quality label from the Q-system is assigned to each cluster by majority vote. This is an example of the opposite of the good discrimination seen in the MWD-features. The Q-classes do not follow the clusters, as seen from the B and C-class samples spread out in all clusters. We can visually inspect the same bad alignment in Fig. 4 and the low adjusted rand score for experiment 0 in Table 3. This indicates that the natural clustering of the rock mass based on MWD features does not align with the Q-classes. We have not assigned the rock type label to these clusters since the overburden thickness and tunnel width feature do not have any intuitive causal relationship with rock type, thus using our domain knowledge. The cumulative distribution plots in Fig. 7 illustrate a broad distribution for Feeder pressure. At the same time, penetration and overburden tend to form three groups, confirming the pattern observed in the median values analysis.

For experiment 3, cluster 3 appears to be an outlier with only 23 samples. The median values are also here clearly separable for each cluster (seen from Table 5 and CDF’s in Fig. 8), with a tendency of three groups for penetration and a well spread of values for feeder pressure and rotation pressure (torque). Rock types have been assigned to each cluster by majority vote. An adjusted rand score of 0.3 indicates a reasonable alignment with the rock type labels. The lowest penetration value for the strong Hornfels rock type and the highest penetration for the weaker Rhomb porfyric rock aligns well with the physical world.

Twenty-nine samples are not classified in experiment seven (see Table 6). These are visible as purple-coloured points in the scatter plot of clusters in Fig. 5 and highlighted by the blue-coloured outlier distribution in for the CDF plots of Feederpressure, Penetration and Rotationpressure in Fig. 9. What is apparent in experiment 7, compared to the other experiments, is the increased separation of cluster properties seen by all the Median values for MWD-features in Table 6 and CDFs in Fig. 9. This indicates a clear separation of rock mass properties in the ten defined clusters. We have assigned rock type labels to the clusters by majority vote. In Fig. 5b, we have coloured the scatter points with rock type labels, clearly illustrating the alignment of rock type to the clusters. This alignment is confirmed by the highest adjusted Rand score of the experiments of 0.48 (a score of 1.0 would be a perfect match between clusters and rock type). Such an alignment is also confirmed by the study of Hansen et al. [45] which forecasted rock type from MWD-data with high predictive accuracy (above 96%).

As for experiment 3, the strong Hornfels are correctly assigned to a cluster with the lowest penetration and the next-to-highest torque. Regarding penetration, we have an overall pattern for big clusters going from Hornfels to Granite, Gneiss and Shale, which seems like an intuitive order. Rhombphorphyry is a lava rock with distinct properties for different lava flows. The samples of Rhombporfhyry are from several lava flows. The fact that Rhombphorphyry has both low and high penetration can be explained by that fact. Augen Gneiss is a special kind of Gneiss with specific properties, with its eye-shaped mega crystals of feldspar in a matrix of quartz and other minerals, quite different from the other classicly layered Gneisses in the other samples. Evaluation only based on penetration is too simple for this strong rock type. Inspecting Feederpressure (note the pronounced CDF-distribution in Fig. 9) and Rotation pressure we see high values, reflecting the strongness of this rock type. We, therefore, can explain its label on cluster number two.

MWD-data, represented as extracted vectors of statistical metrics from thousands of values, can be clustered to well defined clusters. For developing a rock mass classification system, it is essential to thoroughly examine and understand the characteristics of each cluster in relation to the specific problem. The arrangement of classes may vary from weak to strong, or potentially in other dimensions depending on the problem requirements. To ensure the clusters are relevant, one must tune the clusters to your problem to exhibit desired properties. The resulting clusters will vary depending on the selected features. To obtain meaningful clusters, you can:

Utilise a core set of MWD features that serve as a signature for the rock mass and naturally form effective clusters. Repeat the following steps until the clusters align with your problem in terms of cluster properties, cluster count, and cluster definition (well-defined and separable). Subsequently, label your clusters with descriptive names that reflect the categories in your new classification system.

1. 1.

   Incorporate features known to be important for the problem, ensuring they are easily measurable at the excavation site. Potential features include rock cover, tunnel width, measured water inflow in drillholes, distance to a parallel tunnel, soil cover, distance to a fault, point load value, and distance to a lake. By adding the easily collectable features of rock cover and tunnel width to the core MWD clusters in this study, we significantly enhanced the clustering in both visual and metric terms.

   Further refinement may require additional features to achieve a sufficiently detailed cluster system, possibly to accommodate different failure mechanisms. The principal rock stress components are likely essential, given their significant impact on expected failure types. Unlike the time-consuming experiments needed to measure rock mechanical properties (E-modulus, UCS, Poisson’s ratio), which also face scaling issues, approximating rock stress might be feasible and sufficiently accurate for considering rock failure principles. The addition of overburden thickness, a critical variable in estimating vertical stress, exemplifies such a tuning step (see Section 4.2).

   When clustering your data, you are not constrained by the requirements of supervised learning, where the model must generalise based on features present in both training and prediction phases. This allows the use of any features in your dataset to form clusters, even those not available during real-world prediction. For example, post-blasting metrics like overbreak can be included in clustering but are unsuitable for supervised learning.

2. 2.

   Perform the clustering using the algoritmic pipeline outlined in this study.

3. 3.

   Investigate the properties of the clusters by analysing the distributions of cluster features. Begin by examining the distribution of already logged values, such as joint sets or point load tests, if they are not included in the cluster feature vector. Subsequently, gather additional data from the sample sites in the dataset.

You then need to map actions to the clusters, such as associating support classes with RMR and the Q-system. One method involves mapping empirical data of rock support from stable sites to the clustering value of that site (e.g. complementary to support classes for stability classes A, B, C…in the Q-system). Alternatively, you could incorporate descriptions of installed support from stable sites into the cluster feature. This approach would allow clusters to encompass information about specific rock masses and the corresponding stable rock support. When a new rock mass is encountered during excavation, the appropriate rock support can be directly inferred from the cluster. Other methods may also be viable, but exploring these is beyond the scope of this study.

To implement the clustering-based system, you need to build a supervised learning model for classifying a new rock mass encountered during excavation to the correct cluster label. Sapronova et al. [15] outlines a relevant approach. Typically, you would use pre-blasting data available for clustering. If the classification process identifies the rock mass as an outlier not covered by the system, the sample should be closely inspected. If no anomalies are found, it indicates the presence of a rock mass not included in the system, suggesting the need for system updates through reclustering.

Using these principles, you can develop a data-driven system tailored to your specific problem. This study focuses on rock mass stability, primarily concerning permanent or advance rock support. Alternatively, another system could concentrate on blasting design, grouting effort (grouting volume, pumping time), or water leakage, selecting features and designing clusters accordingly. Following the successful alignment of rock type to clusters, as detailed in experiment 7 in Section 4.3, we describe a concrete approach below to initiate the design of an intuitive classification system.

The clear alignment between rock types and clusters, suggests the potential for developing an intuitive rock mass classification system. Ground control engineers intuitively understand and accept the differences in required rock support between rock types, such as blocky versus tightly jointed granite, or a brittle, porous Rhombphorphyry versus a massive one. This understanding can extend to expectations of water problems or drill bit wear. By labelling the clusters (i.e. the rock mass classes) with descriptive names like "Massive Granite", "Blocky Granite", and "Tightly Jointed Granite", it becomes easier to determine the appropriate rock support for each class. Although this approach may result in a large number of classes, such descriptive labels are likely more intuitive and meaningful for engineers than existing alpha numeric systems. An engineer knows the implications of a blocky granite in terms of rock support, water leakage, etc., but may not fully understand what a Q-class C rock mass entails.

In this study we have clustered MWD-data, framed as vectors of extracted statistical metrics, collected from thousands of values in one meter sections of a full face blasting round. Conceptually the approach might be used on all kinds of MWD-data, from drilling boltholes, exploratory holes, grouting, and blasting holes, provided there is a sufficient quantity of MWD-values to extract the statistical vector that characterises the rock mass. For instance, the vector could be computed for $1\times 1\times 1$ meter cubes of data, $4\times 4\times 4$ meter cubes, a 1-meter long split profile section, or 1-meter sections in single holes. Determining the limits of this application is a subject for future research.