Simultaneous Prediction and Optimisation of Rock Fragmentation and Ground Vibration Using an ANN-RF Ensemble in Open-Pit Blasting

1. Introduction

The prevalent method employed in mining, quarrying, and civil engineering projects to break down rocks is blasting. Only around 20-30% of the energy produced is utilised for the desired rock fragmentation, while the remaining energy is wasted on ground vibration, airblast, noise, flyrock, and backbreak. [1,2,3,4,5,6]. The severity of these undesirable effects results in human discomfort, endangers human safety, and destroys nearby structures and equipment within the blast zone. [7,8,9]. Achieving an appropriate rock size is advantageous for reducing costs associated with loading, hauling, crushing, and grinding operations [10,11,12,13,14]. Consequently, accurate forecasting of rock fragmentation and ground vibration plays a vital role in the mining industry. This enables the optimisation of the particle size distribution of the blasted rock and minimises the environmental impact caused by ground vibration.

Several methods are used for measuring fragmentation including sieve analysis, shovel loading rate, oversize boulder count and image analysis. Strain waves from a blast propagating as elastic waves at ground level are known as ground vibrations. Ground vibration is evaluated based on its frequency, acceleration, displacement, and peak particle velocity (PPV) [15].

Blast impacts are typically influenced by two factors: controllable and uncontrollable parameters. Controllable parameters can be adjusted, such as blast design and explosive parameters. On the other hand, uncontrollable parameters cannot be modified, such as the geological and geotechnical properties of the rock [16,17,18]. Researchers have developed various empirical methods to estimate rock fragmentation and ground vibration [19,20,21,22,23,24]. These models, however, face limitations due to their inability to incorporate the effects of additional parameters comprehensively. Typically, empirical methods rely on a limited set of input variables and produce a single output, resulting in reduced accuracy. This approach makes it challenging to simultaneously consider all the relevant parameters, especially given the complex and non-linear interrelationships among them.

As a result, various artificial intelligence methodologies have been adopted in numerous engineering studies, including blasting, to address the shortcomings associated with traditional empirical approaches [25,26,27,28,29,30,31,32,33,34]. Zhou et al.[35] estimated the particle size distribution resulting from blasting operations by analysing 88 blasting datasets from two quarries located in Iran. They utilised several computational techniques, namely an adaptive neuro-fuzzy inference system (ANFIS) optimised with the firefly algorithm (FFA), genetic algorithm (GA), support vector regression (SVR), and artificial neural network (ANN). Among these models, the ANFIS-GA demonstrated superior predictive capabilities, achieving the highest coefficient of determination ($R^2$) value of 0.989 and the lowest root mean square error (RMSE) value of 0.974. Furthermore, Fang et al.[36] proposed an approach for predicting rock fragmentation outcomes by combining the boosted generalised additive model (BGAM) with the firefly algorithm (FFA). This integrated FFA-BGAM framework was benchmarked against other soft computing approaches, including FFA-ANN, FFA-ANFIS, support vector machine (SVM), Gaussian process regression (GPR), and k-nearest neighbors (k-NN). The comparative analysis indicated that the FFA-BGAM method provided superior performance in forecasting rock fragmentation compared to the alternative models.

Huang et al. [37] predicted rock fragmentation using auto-tuning model, called cat swarm optimisation (CSO) and particle swarm optimisation (PSO) algorithm. The input parameters considered for the study were spacing, burden, stemming, maximum charge per delay, rock mass rating, and specific charge. The CSO algorithm outperformed the PSO algorithm in estimating rock fragmentation with a RMSE of 0.847. Sensitivity analysis from the same study revealed that stemming had the most influence on fragmentation.

Amiri et al. [38] introduced an integrated approach combining artificial neural network (ANN) and k-nearest neighbours (k-NN) algorithms for the simultaneous prediction of blast-induced ground vibrations and air overpressure. They evaluated the performance of this combined ANN-kNN framework against an individual ANN model and two conventional empirical formulas established by the United States Bureau of Mines (USBM). Their findings revealed that the proposed ANN-kNN hybrid model delivered enhanced predictive accuracy for both ground vibration and airblast, outperforming the other tested approaches. Zhou et al. [39] employed Bayesian network (BN) and random forest (RF) methodologies to estimate ground vibrations induced by blasting, utilising data from 102 blasting operations. To enhance predictive accuracy by reducing the dataset dimensionality, a feature selection (FS) method was implemented, initially comprising eleven parameters. Subsequently, the FS approach identified five critical parameters, power factor, hole depth, maximum charge per delay, stemming, and distance from the blast face, as essential for accurate ground vibration predictions. Comparative analysis indicated that the RF model exhibited greater precision in forecasting ground vibrations than the BN model.

Yang et al. [40] employed an adaptive neuro-fuzzy inference system (ANFIS) optimised by genetic algorithm (GA) and particle swarm optimisation (PSO) for predicting ground vibration. Their findings indicated that the ANFIS-GA variant achieved superior predictive performance, resulting in a significant 61% reduction in root mean square error (RMSE) and a 10% improvement in the coefficient of determination ($R^2$), compared to the conventional ANFIS approach.

Chandrahas et al.[41] simultaneously predicted rock fragmentation and blast-induced ground vibrations utilising a dataset comprising 152 blast instances from Opencast mine I located at Ramagundam III Area of Singareni Collieries Company Limited, Telangana, India. They implemented several machine learning algorithms, including extreme gradient boosting (XGBoost), k-nearest neighbors (k-NN), and random forest (RF). Their study indicated that XGBoost outperformed the other models, demonstrating optimal values of mean absolute percentage error (MAPE) (22.5), RMSE (3.873), and $R^2$ (0.9125) for fragmentation prediction, as well as MAPE (18.4), RMSE (2.8890), and $R^2$ (0.9125) for ground vibration estimation. The inputs they analysed included spacing to burden ratio, maximum charge per delay, total explosive quantity, firing pattern, and joint angle degree. To the best of our knowledge, Chandrahas et al. [41] is the only researcher who has made a simultaneous prediction of rock fragmentation and ground vibration. This work is closely compared with his, as it also considers the simultaneous prediction of blast-induced rock fragmentation and ground vibration. Therefore the main contributions of this paper are the following:

• The research gathered a dataset comprising 120 blasting events, sourced from the blast records of the Jwaneng Diamond Mine in Botswana.
• This work takes into account ten input variables, which include parameters related to blast design, explosive properties, and rock mass characteristics.
• Three data-driven models are proposed for the simultaneous prediction of blast-induced rock fragmentation and ground vibration. These models are random forest (RF), ANN, and the ensemble ANN-RF.
• We optimise the architecture of the best-performing machine learning model using the Monte Carlo method. And create a solution surface from the optimised machine-learning model.
• Optimisation of input parameters to maximise fragmentation and minimise ground vibration is conducted using the gradient descent method from the created solution surface. The solution surface can be used for prediction, optimisation and finding the inverse solution by setting the desired value of fragmentation and ground vibration and searching in the solution space for the values of the corresponding input parameters.
• Feature importance analysis is performed using the optimised random forest model and the results are confirmed from the created solution space.

2. Materials and Methods

This section explores the materials and methodologies employed in this study, focusing on the datasets collected from the Jwaneng Diamond Mine, the machine learning techniques applied, the feature importance analysis conducted, and the optimisation processes implemented.

2.1. Materials

A compilation of 120 blast datasets sourced from the mining records in Jwaneng has been used to train and test the models proposed in this research. Figure 1 shows data collected from fragmentation images and seismographs. Split Desktop software was used to analyse the muckpile images captured with a digital camera, as shown in Figure 1a, and to generate the particle size distribution (PSD) curve for fragmentation. The resulting PSD curve is presented in Figure 1b. The quality of fragmentation at the mine is assessed based on 80% passing the 150 mm sieve size, less than 50 mm is considered undersize while above 150 mm is considered oversize. Ground vibration data in this study were recorded using seismographs placed at five different locations from the blast sources as shown in Figure 1c. The seismographs measured ground vibrations in three orthogonal components: radial, transverse, and vertical. Figure 1d shows a sample siesmograph report. The radial component captures vibrations in the direction of wave propagation from the blast source, the transverse component measures vibrations perpendicular to this direction within the horizontal plane, and the vertical component records vibrations in the up-and-down direction. Together, these three components provide a complete three-dimensional representation of ground motion, enabling accurate assessment of vibration intensity and directionality. The data used in the study span the period from 2017 to 2024, covering all blasting conditions and potential biases in these seven years. The parameters under consideration, along with their respective ranges, are outlined in Table 1. The data was preprocessed which included data cleaning and normalisation. The dataset was then divided into a training set (80%) and a testing set (20%).

2.2. Methods

According to the literature, there is a consensus among researchers that ANNs are effective for predicting blast-induced impacts. However, Yan et al. [42] highlighted their prolonged training times, while Jadav et al. [43] emphasised their susceptibility to local minima. To address these challenges and enhance prediction accuracy, this study proposes am ensemble ANN-RF model. Additionally, standalone ANN and RF models are utilised for comparison purposes. The model assessment was conducted using RMSE and R2 as performance indices. Sensitivity analysis was performed using the network weights of the best-performing model. The gradient descent method was used to optimise the outputs, and the Monte Carlo method was applied to find the optimal architecture of the best-performing model.

2.2.1. ANN

In this research, hyperparameters for an artificial neural network (ANN) model were optimised using the grid search technique. Hyperparameters selected for evaluation comprised learning rate values (0.001, 0.01, 0.1), configurations of hidden layers (20, 32, 64, 100, (50 50), (65 100), (25 75)), activation functions (rectified linear unit (RELU), hyperbolic tangent (tanh), sigmoid), batch sizes (32, 64), optimisers (Adam (adaptive moment estimation), stochastic gradient descent), and epoch counts (100, 200, 300). To facilitate adaptive ANN architectures capable of adjusting to varied hyperparameter combinations, a dynamic model generation function was created. The grid search was carried out utilising the GridSearchCV module from Scikit-learn, systematically testing all possible hyperparameter combinations. The dataset was partitioned into an 80% training set and a 20% testing set. For every hyperparameter set, the ANN was trained on the training subset and its predictive accuracy was validated against the testing subset to ensure robustness and generalisability. Model effectiveness was measured through loss metrics, and the hyperparameter combination that delivered the highest accuracy on the testing set was chosen as optimal. The loss trajectory for the optimal ANN configuration is presented in Figure 2. It shows stable convergence behaviour, confirming effective learning without overfitting, as indicated by the close alignment of validation loss with training loss, demonstrating strong generalisation capabilities to new data.

2.2.2. RF

In this research, hyperparameter tuning of the random forest (RF) algorithm was carried out employing a random grid search method to enhance its predictive performance. The hyperparameters considered included the number of estimators (50, 100, 200, 300), maximum tree depth (10, 20, 30, None), the minimum number of samples needed to split a node (2, 5, 10), the minimum number of samples required at each leaf node (1, 2, 4), and the bootstrap sampling approach (bootstrap: True, False). A randomised combination of these hyperparameters was utilised, which allowed for computationally efficient parameter exploration compared to exhaustive grid search methods.

The RandomizedSearchCV functionality from Scikit-learn was implemented to execute this randomised search across 100 iterations. The dataset was divided into an 80% training subset and a 20% testing subset. Each hyperparameter combination was used to train the RF model on the training subset, and its predictive accuracy was assessed using the testing subset. Model evaluation was conducted using root mean squared error (RMSE) and the coefficient of determination ($R^2$) as performance criteria. The optimal hyperparameter set was identified based on achieving the lowest RMSE and the highest $R^2$ value on the test data. This strategy effectively balanced computational efficiency and thorough hyperparameter optimisation. The learning curve of the optimal RF configuration is depicted in Figure 3, illustrating a distinct reduction in training and validation errors as the size of the training set increased. This trend confirms that augmenting the data volume enhances the predictive accuracy of the model. Additionally, the validation loss closely mirrors the decrease in training loss, signifying robust model generalisation.

2.2.3. ANN-RF

In this study, an ensemble approach combining an ANN and RF was implemented, utilising the optimal parameters identified during independent hyperparameter tuning of the ANN and RF models. The optimal parameters for the ANN were the learning rate (0.01), number of neurons (64), activation function (RELU), batch size (32), optimiser (Adam), and the number of epochs (200). The optimal parameters for the RF model were the number of estimators (100), maximum depth (20), minimum samples required to split a node (5), minimum samples required to be at a leaf node (2), and the bootstrap method for sampling (bootstrap: True). Using the best parameters from both methods, an ensemble ANN-RF model was constructed. The ANN and RF models were trained separately, with their predictions combined using a voting regressor to achieve optimised regression outcomes. This approach leverages the strengths of both algorithms, with the voting regressor integrating their outputs to improve overall predictive performance. The combined model was trained on the training dataset, comprising 80% of the data, and evaluated on the testing dataset, comprising the remaining 20%. Model performance was assessed using RMSE and the R2 as performance metrics. The ensemble model leveraged the complementary strengths of the ANN for capturing complex feature representations and the RF for robust decision-making. The flowchart illustrating the ANN-RF method is shown in Figure 4.

2.2.4. Feature Importance

In this study, a feature importance analysis was conducted to quantify the influence of each input parameter on the model outputs, namely ground vibration and rock fragmentation. The optimal RF model was employed due to its robustness, ability to capture nonlinear relationships, and inherent support for feature importance evaluation. Feature importance was computed using the mean decrease in impurity (MDI) method, which measures the average reduction in impurity, that is, the variance for regression introduced by each feature across all trees in the ensemble. The general formulation for computing the MDI-based feature importance $FI\left(f_j\right)$ of input feature $f_j$ is expressed in (1):

$$FI\left(f_j\right)={\displaystyle \frac{1}{T}}\sum _{t=1}^T\sum _{n\in N_t\left(f_j\right)}{\displaystyle \frac{N_n}{N}}·\Delta i\left(n\right)$$

(1)

In Equation (1), T denotes the total number of trees in the forest, and $N_t\left(f_j\right)$ represents the set of internal nodes in tree t where feature $f_j$ is used to split. $N_n$ is the number of training samples reaching node n, and N is the total number of samples. The term $\Delta i\left(n\right)$ refers to the impurity reduction achieved by the split at node n.

For regression problems, impurity is typically measured using variance. The impurity reduction $\Delta i\left(n\right)$ at a node can therefore be calculated as shown in (2):

$$\Delta i\left(n\right)=\mathrm{Var}\left(\mathrm{parent}\right)-\left({\displaystyle \frac{N_{\mathrm{left}}}{N_{\mathrm{parent}}}}·\mathrm{Var}\left(\mathrm{left}\right)+{\displaystyle \frac{N_{\mathrm{right}}}{N_{\mathrm{parent}}}}·\mathrm{Var}\left(\mathrm{right}\right)\right)$$

(2)

Here, $\mathrm{Var}\left(\mathrm{parent}\right)$ is the variance of the target variable at the parent node, while $\mathrm{Var}\left(\mathrm{left}\right)$ and $\mathrm{Var}\left(\mathrm{right}\right)$ are the variances of the left and right child nodes, respectively. $N_{\mathrm{parent}}$, $N_{\mathrm{left}}$, and $N_{\mathrm{right}}$ represent the number of samples at the parent, left, and right nodes. This equation quantifies how much the split at a given node reduces the variance in the target variable, thereby contributing to the overall importance of the feature used in that split.

To enhance interpretability, the computed feature importances were normalised such that their values sum to 1. This allows each importance value to be expressed as a percentage of the total, enabling a clearer comparison of the relative contribution of each feature. Formally, the normalised feature importance $F{I}^{*}\left(f_j\right)$ for feature $f_j$ is computed as:

$$F{I}^{*}\left(f_j\right)={\displaystyle \frac{FI\left(f_j\right)}{{\sum }_{k=1}^{K}FI\left(f_k\right)}}\times 100$$

(3)

where K is the total number of input features, and $FI\left(f_j\right)$ is the raw feature importance obtained from Equation (1). The normalised importances $F{I}^{*}\left(f_j\right)$ sum to 100, making them suitable for percentage-based sensitivity comparisons, as illustrated in Figure 6. This approach provides an intuitive and computationally efficient way to rank input variables based on their contribution to reducing prediction uncertainty. The resulting feature importances were used to interpret the model and guide further analysis.

2.2.5. Gradient Descent Optimisation

This study utilises gradient descent optimisation to determine the optimal values for fragmentation and ground vibration, along with their corresponding input parameters. The Monte Carlo method was applied to find the optimal model architecture, which was determined to be 10-60-45-2. The gradient descent method was employed to navigate the solution surface, maximising fragmentation while minimising ground vibration, thereby identifying the optimal input parameters. The specific technique employed for gradient descent is shown in (4):

$$p_i^{(k+1)}=p_i^{\left(k\right)}-\eta \left({\displaystyle \frac{\mathrm{ANN}-\mathrm{RF}(p_1^{\left(k\right)},\dots ,p_i^{\left(k\right)},\dots ,p_N^{\left(k\right)})-\mathrm{ANN}-\mathrm{RF}(p_1^{\left(k\right)},\dots ,p_i^{(k-1)},\dots ,p_N^{\left(k\right)})}{p_i^{\left(k\right)}-p_i^{(k-1)}}}\right)$$

(4)

In this equation, $p_i^{(k+1)}$ represents the updated parameter at iteration $k+1$, while $p_i^{\left(k\right)}$ is the current parameter value at iteration k. The term $p_i^{(k-1)}$ denotes the parameter value from the previous iteration. The learning rate is represented as $\eta$. The function ANN−RF represents the ANN-RF model.

3. Results and Discussion

To check the performance of all the predictive models RMSE and R2 were utilised as performance indices. These are expressed by Equations (5) and (6):

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{i=1}^N{(y-y\prime )}^2}$$

(5)

$$R^2={\left[{\displaystyle \frac{{\sum }_{i=1}^N(y-\overline{y})(y\prime -\overline{y}\prime )}{{\sum }_{i=1}^N{(y-\overline{y})}^2{\sum }_{i=1}^N{(y\prime -\overline{y}\prime )}^2}}\right]}^2$$

(6)

where y and $y^{\prime }$ are the measured and predicted values, respectively; $\overline{y}$ and ${\overline{y}}^{\prime }$ are the mean measured and mean predicted values, respectively.

Results from the chosen models (ANN, RF, and ANN-RF) for estimating fragmentation and ground vibration, along with their performance indices, are outlined in tab-training,tab-testing for the training and testing set, respectively. The high accuracy and low error of the models in Table 2 show that they have learned the data well, especially the ANN-RF model. Furthermore, the performance achieved on the test datasets, as presented in Table 3, suggests that the developed models are well generalised. The ANN-RF model demonstrates superior performance on the testing set as well, as indicated by the following R2 values: 0.956 for fragmentation, and 0.930 for ground vibration. Additionally, the RMSE values are 0.315 for fragmentation, and 0.380 for ground vibration. These results suggest that the ANN-RF model achieves a lower system error and higher accuracy in simultaneously modelling fragmentation and ground vibration than all the other implemented models.

Figure 5a,f depict the predicted versus actual values of fragmentation and ground vibration for the models used in this study within the testing datasets. Generally, all the models demonstrate a high predictive accuracy, but Figure 5e,f for ANN-RF model shows that the predicted values closely match the measured ones, following the ideal y:x line, which shows a strong correlation compared to other models. This confirms the effectiveness of the ANN-RF model in accurately predicting these outcomes. The superior performance of the ANN-RF hybrid model can be attributed to the complementary strengths of its components. The ANN (see Figure 5a,b) outperformed the RF (see Figure 5c,d) due to their ability to capture non-linear relationships and extract complex features through their layered architecture, which is particularly effective for datasets with intricate variable interactions like the dataset considered in this study. However, ANNs are prone to overfitting and lack interpretability, while RF excels in generalisation by handling structured patterns and mitigating overfitting through ensemble learning. The hybrid model leverages these advantages by using the ANN for feature extraction and the RF for robust decision-making, resulting in improved predictive accuracy. This enables the hybrid approach to address the limitations of each standalone model and achieve superior performance and robustness.

3.1. Feature Importance Analysis

Figure 6 illustrates the results of the sensitivity analysis conducted for both fragmentation and ground vibration. For fragmentation, powder factor (17%) emerges as the most influential input parameter. This finding underscores its critical role in determining the energy applied per unit volume of rock, which directly impacts the extent of breakage during blasting. Burden (15%) follows closely, reaffirming its importance in controlling the confinement and direction of explosive energy, both of which influence fragmentation outcomes. Other notably contributing factors include stemming and rock factor, indicating that both explosive confinement and rock mass characteristics significantly affect fragmentation. On the contrary, the blastability index, with a contribution of only 6%, was found to be the least influential parameter for fragmentation. This suggests that, within the context of this dataset, rock breakability characteristics play a secondary role compared to design-driven variables.

For ground vibration, the analysis identifies burden (15%) as the most impactful parameter. This highlights the burden’s key function in moderating explosive energy dissipation, where larger burdens can absorb more energy, thereby reducing the intensity of vibrations transmitted through the ground. Hole diameter (6%) is observed to have the least effect on ground vibration, indicating that either its variation is limited across the blast designs, or that its contribution is indirect and less significant compared to other parameters. Parameters such as maximum charge per delay, powder factor, and stemming also demonstrated moderate influence, suggesting their relevance across both performance metrics.

These results demonstrate that certain parameters, particularly powder factor and burden are critically influential in controlling both fragmentation and ground vibration, making them essential targets in blast design optimisation. In contrast, parameters such as blastability index and hole diameter may be of lesser priority, particularly under consistent geological conditions. These observations are consistent with findings from previous studies [39,44,45,46,47,48], which similarly report the dominant effects of explosive energy distribution and confinement on blast-induced outcomes.

Figure 6. Sensitivity analysis results

Figure 6. Sensitivity analysis results

3.2. Analysis of the Optimisation Results

Figure 7a,b display the results of the optimisation procedure. In each figure, the red dots represent the randomly selected initial points, with seven such points chosen. Applying the gradient descent method to each initial point led to convergence at the optimal solution, depicted by the blue dot. The optimised values for fragmentation and ground vibration are approximately 84% and 0 mm/s, respectively. Table 4 shows the optimised input parameters. Variations in optimised parameters are influenced by the model’s complexity, random initialisation, and interpolation errors. The model’s inherent complexity, with numerous interdependent variables and non-linear relationships, makes it sensitive to initial conditions, causing fluctuations in outcomes. Additionally, the random initialisation of network weights results in different starting points for each training run, leading to varied convergence paths and optimised parameters. Furthermore, interpolation errors arise from the grid resolution used to identify optimal points; the grid is not sufficiently fine, hence it misses more precise locations, resulting in inaccuracies in the final parameters.

From the sensitivity analysis, blastability index and hole diameter are the least influential input parameters on fragmentation and ground vibration, respectively. Because of this, we compare these parameters versus all the other input parameters kept constant, to observe the variations in the values of the output parameter in a two dimensional plot. We assign blastability index and hole diameter to be the x-axis, while ground vibration is assigned as the y-axis.

In Figure 7a, when the powder factor increases while the blastability index is fixed, rock fragmentation value increases towards an optimal value then decreases with fluctuations. In Figure 7c, increasing the hole diameter while keeping the burden fixed results in ground vibrations initially increasing then decreases to an optimal value before increasing again.

The solution space also confirms the results of the feature importance analysis. As can be seen in Figure 7b, generated by varying blastability index which is the least influential input parameter on fragmentation and fixing all the other parameters, and taking a slice through the lowest, optimised point in the solution space. Powder factor is the most influential input parameter on fragmentation, and exhibits significant fluctuations over the blastability index range compared to other parameters. Figure 7d, shows that burden has the most significant variations compared to other parameters, suggesting a higher sensitivity to ground vibration.

This multi-output optimisation of fragmentation and ground vibration has significant practical applications. It allows for the simultaneous consideration of these critical factors, enabling a balanced approach to blasting that enhances operational efficiency while minimising environmental impacts. By optimising fragmentation, mining operations can achieve the desired material size, which improves the efficiency of downstream processes like loading, hauling, and milling. At the same time, controlling ground vibration helps to reduce the risk of damage to nearby structures, ensures compliance with regulatory standards, and mitigates the environmental impact of blasting activities. This approach leads to more efficient and safer mining operations, ultimately contributing to cost savings, improved productivity, and reduced environmental footprint.

4. Conclusions

This study demonstrates the successful application of a hybrid machine learning technique to simultaneously predict and optimise blast-induced rock fragmentation and ground vibrations. The method was applied at the Debswana Jwaneng Diamond open-pit mine, utilising 10 input parameters. The evaluated models included ANN, RF, and the ensemble ANN-RF. Among these, the ANN-RF model proved to be the most effective, achieving the highest predictive accuracy on the test set, with $R^2$ values of 0.956 for fragmentation, 0.930 for ground vibration. Correspondingly, the RMSE values were 0.315 for fragmentation, 0.380 for ground vibration. The architecture of the ANN-RF model was optimised using the Monte Carlo method, resulting in an architecture with two hidden layers containing 60 and 45 neurons, respectively. The model was further refined using the gradient descent method to simultaneously minimise ground vibration (approximately 0 $mm/s$) while maximising fragmentation (84%). Feature importance analysis revealed that powder factor (17%) was the most influential parameter for fragmentation, whereas the blastability index (6%) had the least impact. For ground vibration, burden (15%) emerged as the most significant factor, while hole diameter (6%) was the least influential. While this study focuses on the Jwaneng Diamond Mine, the methods used are scalable and adaptable to other mining operations. With site-specific data and model retraining, the approach can be applied to mines with varying topographical, geological, and operational conditions.