The first step in our methodology involved the collection and preprocessing of spatial data (Bai and Tahmasebi, 2021) relevant to experimental variogram modelling (Pesquer et al., 2011). For this study, we obtained data from four orebodies characterized by similar genetic aspects (geological background). These orebodies were chosen to ensure consistency in the spatial characteristics of the data (Li Z. et al., 2018), facilitating meaningful comparisons in our analysis. The spatial data included measurements of ore grades at various locations within each orebody (Liu et al., 2022). These measurements were used to compute experimental variograms (Pardo-Igúzquiza and Dowd, 2001), which quantify the spatial dependence between pairs of data points (Fouedjio, 2016). To warrant the data reliability and accuracy, we performed careful quality control processes, including outlier detection and data cleaning. The used dataset in this research includes a medium-sized database containing 243,808 composite samples from four orebodies, all sharing the same geological characteristics (McCormick and Heaven, 2023) as the misinformed orebody, extracted from 477 drillholes. The assays encompass 16 variables, including sample coordinates (northing, easting, and elevation), ore grades and sample length. Sampling was performed at both regular and irregular intervals, with composite data at a 5 m sampling interval (Figure 2).

Once the spatial data were collected and preprocessed, we proceeded to model experimental variograms (Pardo-Igúzquiza et al., 2013) for each orebody. Experimental variograms were computed using the traditional method of pairwise differences, where the variance of the differences between data points, at different distances is calculated (Rivoirard, 2007). This procedure grants helpful perceptions into the data spatial structure and variability, which are fundamental for successive predictive modeling (Niu et al., 2024). It’s essentially half the anticipated square deviation among pair off random functions (Lui et al., 2022), Z x and Z x + h situated at a certain space and bearing vector (with anisotropy factored in), or else what’s termed as lag h Equation 1: γ h = 1 2 E Z x − Z x + h 2 (1)

To compute experimental variograms, we followed standard procedures outlined in the geostatistics literature (Atkinson and Lloyd, 2007). Specifically, we calculated the semivariance between pairs of data points at various lag distances, using a predefined lag tolerance to ensure enough data pairs for reliable estimation (Afeni et al., 2021). The resulting experimental variograms were then plotted and analyzed to identify spatial trends and patterns. Instead of continuous variables, the “experimental semi-variance” stays described as quasi of the mean square off variation among quantities that are a certain lag h apart. This means that the γ h (experimental variogram) could be derived for α = 1 , 2 , 3 , … . . , N h couples of samples Z x α , Z x α + h at positions x α , x α + h split with a constant lag h Equation 2: γ h = 1 2 N h ∑ α = 1 N h z x α − z x α + h 2 (2)

A mathematical model can be applied to the variogram, and its coefficients can find the best weights for spatial prediction through Kriging. The model must be conditionally negative semi-definite, as emphasized by (Atkinson and Lloyd, 2007). Typically, the model is selected from a set of approved or valid models that meet this criterion, as discussed in a review by (Li Z. et al., 2018) of commonly used valid models as spherical model Equation 3 ( C 0 : Nugget effect, C : Sill, a : Lag, h : Distance between samples). γ h = C 0 + C .   3 2   h   a − 1 2   h   3 a 3   ;   0 <   h   ≤ a C 0 + C   f o r     h   > a   w i t h   C , a > 0 (3)

With the experimental variograms computed, we proceeded to develop neural network regressors (Li X. et al., 2018) for predicting variogram based on input data using conceived MATLAB scripts. Neural networks are impressive machine learning models (Friedman, 2001; Manouchehrian et al., 2012) able of catching complex relationships in data (LeCun et al., 2015), through interconnected layers of neurons. In our case, we used feedforward neural networks (Nwaila et al., 2024), which involve of an interconnected layers divided into three input layers, one or many hidden layers, and one output layeryperparameters to evalua Figure 1.

The architecture of the neural network regressors (Heaton, 2018) was carefully designed to optimize predictive performance while minimizing computational complexity (Lozano et al., 2011; Adeniran et al., 2019; Lundberg et al., 2020). We tested with many configurations, counting different numbers of hidden layers, activation functions, neurons per layer, and regularization techniques (Kim et al., 2023). These configurations were chosen based on practical evidence and domain expertise to confirm the neural networks efficiency in experimental variogram modelling.

To fine-tune the hyperparameters (Ystroem et al., 2023) of the neural network regressor, we engaged Bayesian optimization (Zhang et al., 2020), an overwhelming optimization method that powers probabilistic models to guide the quest for optimal hyperparameters (Asante-Okyere et al., 2022). Bayesian optimization runs iteratively, operating previous valuations to update its probabilistic model and select the next set of hyperparameters to evaluate (Xie et al., 2022; Rong et al., 2023).

In our implementation of Bayesian optimization (Shahriari et al., 2016), we used Gaussian process regression (Arabpour et al., 2019; Phelps and Cronkite-Ratcliff, 2023) to model the objective function, which in this case was the performance of the neural network regressor in predicting experimental variogram. We defined appropriate acquisition functions (Zhang et al., 2020; Hallam et al., 2022), such as expected improvement or probability of improvement, to guide the search for optimal hyperparameters efficiently.

To evaluate the performance of the Bayesian-optimized and the others neural network regressors, we conducted rigorous validation experiments using a holdout dataset (Adeniran et al., 2019). The dataset was randomly split into training and validation groups, ensuring that each set comprised a representative data sample (Guo et al., 2022).

The neural network regressors were trained on the training set using the best hyperparameters achieved among Bayesian optimization (Zhang et al., 2021). The accomplished models were afterward, evaluated on the validation set, using proper performance metrics (Wu and Zhou, 1993), like: Mean squared error (MSE), Coefficient of determination (R²), Mean absolute error (MAE).

Finally, we conducted a comparative analysis (Soltanmohammadi and Faroughi, 2023) to assess the performance of the Bayesian-optimized neural network regressor against non-optimized configurations (Pavlov et al., 2024). We compared the predictive accuracy of Bayesian-optimized neural network regressor with wide, bilayer, and tri-layer networks (Lauzon and Marcotte, 2022; Kim et al., 2023).

The comparative assessment engaged quantitative valuation of the performance metrics (Li Z. et al., 2018; Liu et al., 2022), as well as qualitative assessment of the predictive models accuracy and robustness (Dutta et al., 2010). By comparison of different models performance, we aimed to prove the lead of Bayesian-optimized neural network regressor (Fasnacht et al., 2020; Houshmand et al., 2022; Costa et al., 2023) in experimental variogram modelling.

Statistical assessment of the composites (Pardo-Igúzquiza et al., 2013) indicated a notably low difference in ore grades, with a mean of 0.44% and a standard deviation of 0.76%. However, the coefficient of variation exceeded one. To ease training of the artificial neural network (ANN) (Hu and Shu, 2015), the data was normalized using log transformation. Figure 3 illustrates histograms of normalized and clustered data for grades of four deposits. Chart inspection of the histograms indicates that the data predominantly consists of medium-grade values, with only a small percentage of very high-grade across all deposits Figure 3.

Following data analysis, we investigated spatial continuity by creating variogram models (Phelps and Cronkite-Ratcliff, 2023). Both omni-directional and directional variograms and are crucial in spatial analyses (Shi and Wang, 2021). However, in our case, we focused on constructing a downhole variogram for each orebody, using a conceived macro in Datamine Studio RM and Supervisor. The dominant direction for each orebody was found, and it was found that the four obtained directions were nearly identical due to the shared genetic context. The spatial structures (Mueller et al., 2020; Souza et al., 2023), as depicted in Figure 4, showed substantial impacts from the nugget effect, implying challenging conditions for variogram modelling (Das et al., 2020).

The obtained variogram models offered improved insight into the deposit, helping in model fitting (de Carvalho and da Costa, 2021). Figure 4 displays the downhole variogram models fitted using a spherical model. A significant portion of the spatial irregularity arising from the nugget effect suggests a medium spatial correlation structure across the study area, as indicated by the variogram plot, demonstrating good spatial correlation (Sharifzadeh Lari et al., 2021).

Experimental variogram modelling is a crucial aspect of geostatistics (Allotey and Harel, 2023), providing insights into the spatial dependence (Fouedjio, 2016) of ore grades within mining environments. In this study, we evaluated four different neural network models for their effectiveness in predicting experimental variograms based on spatial data from multiple orebodies. Here, we provide a complete comparison of assessed models founded on many characteristics and performance metrics Table 1.

The wide neural network (Model 1) (He et al., 2015), with a single fully connected layer comprising 100 neurons and ReLU activation, exhibits moderate performance in experimental variogram modelling. On the validation dataset, it gets an R² of 0.3294 and an RMSE of 0.1318, then proving an adequate data fit. Nevertheless, it is performing on the test dataset is relatively lower, with an RMSE of 0.1461 and a negative R² of −0.7646, advocating overfitting or inadequacy in capturing the principal spatial relationships. Additionally, the model’s Mean Absolute Percentage Error for both validation (52.1646%) and test (28.5176%) datasets show an important discrepancy between predicted and actual values. The absence of regularization in this model may contribute to its susceptibility to overfitting, particularly given the limited architecture complexity Figure 5A.

The bilayered neural network (Model 2), featuring two fully connected layers with 10 neurons each and ReLU activation (He et al., 2015), demonstrates slightly inferior performance compared to the wide neural network. Though it reaches the same RMSE on the validation dataset (0.1492), its R² value is remarkably lower (0.1400), suggesting weaker predictive capability. On the test dataset, however, the bilayered network outperforms the wide network with a lower RMSE (0.1323) and a less negative R-squared value (−0.4468). This suggests that the bilayered architecture may generalize better to unseen data despite its simpler structure. The MAPE values for both validation (42.9603%) and test (22.9784%) datasets remain high, indicating notable prediction errors (Kim et al., 2023) Figure 5B.

The trilayered neural network (Model 3), featuring three fully connected layers with 10 neurons each and ReLU activation, exhibits the weakest performance among the neural network models evaluated. It attains the greatest RMSE on both validation (0.1945) and test (0.1553) datasets, revealing the smallest accurate predictions. The negative R² values on both datasets (−0.4603 on validation, −0.9939 on test) further signify poor model fit. Additionally, the high MAPE values for both validation (45.0308%) and test (27.2294%) datasets highlight substantial discrepancies between predicted and actual values. The trilayered architecture’s increased complexity does not translate to improved performance, suggesting potential issues with model capacity or training convergence (Figure 5C).

The custom neural network (Model 4), optimized through Bayesian optimization (Figure 5D), emerges as the top-performing model for experimental variogram modelling. It achieves the least RMSE on validation and test datasets, 0.0749 and 0.0898, indicating superior predictive accuracy. With three fully connected layers, sigmoid activation, and optimized layer sizes (267, 14, and 3 neurons), Likewise, the model presents a high R² on the validation dataset (0.7836), signifying a strong fit to the data. On the test dataset, although the R² of 0.3335 is lower, it stays positive, implying satisfactory model performance. The MAPE values for both validation (18.2583%) and test (16.9594%) datasets are significantly lower than those of other models, indicating improved prediction accuracy and reduced errors. The absence of data standardization in this model suggests that it effectively handles the input data without requiring normalization, further simplifying the modelling process Figure 6.

In summary, the custom neural network (Model 4) outperforms the wide, bilayered, and trilayered neural network models in experimental variogram modelling. Its greater performance is ascribed to the hyperparameters optimization through Bayesian optimization, resulting in an architecture that effectively captures the underlying spatial patterns in the data. Compared to the other models, the custom neural network demonstrates higher predictive accuracy, stronger model fit, and reduced prediction errors, making it the preferred choice for experimental variogram modelling in mining applications.

Whereas, the study investigates the use of optimized artificial neural networks models, through Bayesian optimization, for experimental variogram modelling in geostatistics, which discussion can cover an investigation of the results, inferences, limits, and future directions of the research.