1 Introduction

In mineral resource extraction, the main objective is meeting production targets in terms of ore tonnage, mineral content and grades. This contributes to the optimal utilization of the mill with regard to throughput and metal recovery, while minimizing costs. Short-term production scheduling and operations control aim to maximize the value of the mining operation by considering the maximum amount of information available at the time. The decision-making is based on models of the spatial distribution of attributes of interest within the ore body.

Typically, grade control models predict attributes relevant to processing and recovery operations at the scale of a selective mining unit (SMU).

The available information consists of exploration holes augmented by chip or channel samples from the grade control process. This database is not only spatially sparse, but also hampered by time delay, as samples are typically analysed offline in a laboratory, and the quality is lower. It can take several days before results are available; the SMU may have already been mined, with decisions about its destination/processing based on an outdated model. This is especially problematic in highly complex deposits.

In many operations, online sensors provide a large amount of georeferenced data during production monitoring. The ability to rapidly incorporate these data into grade control models provides the opportunity to continuously improve the grade control and resource models while mining the deposit. To incorporate production data obtained during the mining process in real time, recent studies have developed a method for backward integration of online production information to improve the accuracy in grade control models (Benndorf 2015). This integration task applies a methodology akin to an ensemble Kalman filter. Methods have been implemented and tested in bulk mining operations in open pit mining settings (Wambeke and Benndorf 2016; Yüksel et al. 2017) and have shown improvements in prediction accuracy. The potential of these methods is obvious, since they are able to assimilate direct and indirect information, improving the resolution at a local model scale.

The updating step for a grade control model is optimal under the conditions of (i) Gaussianity in the non-updated model (prior ensemble), (ii) linearity of the observation operator and the observations and (iii) Gaussianity in the additive observation error. When these conditions are not satisfied, the application of the ensemble Kalman filter in the analysis step is suboptimal, but it can still be satisfactory in a non-linear, non-Gaussian environment (Wikle and Berliner 2007; Amezcua and Van Leeuwen 2014; Carrassi et al. 2018). Hence, the variables under assimilation are positive (can only take values equal to or higher than zero). Thus, zero probabilities should be assigned to negative values of the physical observations, since these cannot be a negative value. To ensure this condition, a Gaussian anarmophosis transformation is commonly used in data assimilation (Bertino et al. 2007; Wambeke and Benndorf 2017; Yüksel et al. 2017; Carrassi et al. 2018). Kumar and Srinivasan (2019) provided a version of an ensemble Kalman filter that can deal with non-Gaussian and non-linear operators by transforming the model parameters into indicator variables.

The focus of this paper is the application of the ensemble-based grade model updating framework in underground mining settings. Here, a generalized mining process is considered, consisting of a cycle that involves drilling, blasting, loading, scaling and changing the support. The updating algorithm has been extended to adapt to the underground mining environment. This refers in particular to the ability to incorporate multiple sensor inputs at different stages of the production cycle, which represent different supports. It has also been tested for three different grade control scenarios. The first scenario shows no conditioning information. In the second scenario, the realizations are conditioned to the drift areas surrounding the block to be extracted, and in the third, borehole information is considered in addition to the data from the second scenario.

Each SMU unit is blasted and loaded onto the conveyor belt, where sensors are installed to analyse the material. Note that, in contrast to the bulk mining applications published so far, in a selective mining setting, this material may represent a blend of ore and dilution. Additionally, before blasting, the SMU may be classified using an online image sensor delivering a source of two-dimensional mine face information. The information obtained from these two types of sensors represents the mean of the whole SMU blasted or part of it. The grade control model, on which basis the operational decisions are performed, has SMU support, while the underlying grid model is typically on a finer scale support. The integration of information needs to account for this change in support. The grid support is converted into SMU support by re-blocking.

This paper is structured as follows: In Sect. 2, the methodology is provided for the updating grade control models based on the online sensor data. The three scenarios are discussed in Sect. 3, followed by the application and results of the updating for the three scenarios in Sect. 4. A discussion and conclusion is provided in Sect. 5.

2 Updating

Data integration can be achieved by an inverse modelling approach (Tarantola 2005). An inverse problem aims to determine the unknown model parameters from the observed state data (Zhou et al. 2014). The basic inverse model is given as

$$\begin{aligned} z = {\mathcal {A}}^{-1}(d) ,\, \end{aligned}$$
(1)

where \({\mathcal {A}}\) is the forward observation operator, which maps the spatial attributes (z) into the real-valued space of sensor observations (d). The operator \({\mathcal {A}}\) can be linear or non-linear.

One inverse modelling approach is the ensemble Kalman filter (Evensen 1994; Bertino et al. 2002; Dubrule 2018; Carrassi et al. 2018), which sequentially incorporates information observed into the considered model. This assimilation method provides an optimal framework for linear operators and Gaussian error distributions (Wikle and Berliner 2007). For that approach, a Markovian assumption is applied to the process. That is, the state at time t, when conditioned on all previous states, depends only on the state at time \(t-1\).

The random function representing the variable of interest at time step t will be denoted by \(\mathbf{Z }_{t}(\cdot )\), being the matrix-valued state variable \(\mathbf{Z }_{t}\). Realizations at \(\mathbf{Z }_{t}\) are updated to \(\mathbf{Z }_{t+1}\) based on sensor readings \(\mathbf{D }_{t+1}\) and evolve over time as

$$\begin{aligned} \mathbf{Z }_{t+1} = \mathbf{Z }_{t} + \mathbf{W }_{t+1}[\mathbf{D }_{t+1} - {\mathcal {A}}(\mathbf{Z }_{t})] .\, \end{aligned}$$
(2)

Here, \(\mathbf{W }_{t+1}\) is a weight matrix that balances new observations obtained against the predicted model. The difference between model-based predictions (\({\mathcal {A}}(\mathbf{Z }_{t})\)) and observations (\(\mathbf{D }_{t+1}\)) is the innovator term. Observations is a replicated vector of observations \(\mathbf{d }_{t+1}\) into the matrix \(\mathbf{D }_{t+1} = [\mathbf{d }_1,\dots ,\mathbf{d }_i]_{t+1}\). The rows represent the number of observations and the columns the number of realizations. The weight matrix is given by

$$\begin{aligned} \mathbf{W }_{t+1} = \mathbf{C }_{t,\mathbf{zz} }\mathbf{A }_{t+1}^\text {T}\big (\mathbf{A }_{t+1}\mathbf{C }_{\text {t,}{} \mathbf{zz} }\mathbf{A }_{t+1}^\text {T} + \mathbf{R }\big )^{-1} ,\, \end{aligned}$$
(3)

where the matrix \(\mathbf{R }\) represents the precision of the sensor measurements. Since measurement errors are assumed to be independent (uncorrelated), the matrix R is taken to be the scalar matrix \(\mathbf{R } = \sigma ^2 \mathbf{I }\), where \(\sigma ^2\) is the variance in the measurement error and I is the identity matrix of size \(n \times n\). The matrix \(\mathbf{A }\) is the linear approximation of the operator \({\mathcal {A}}\) (Wambeke and Benndorf 2017), and the term \(\mathbf{C }_{t+1,\mathbf{zz} }\) is the updated covariance matrix at time \(t+1\).

In further time steps, the covariance evolves as

$$\begin{aligned} \mathbf{C }_{t+1, \mathbf{zz} }&= (\mathbf{I } - \mathbf{W }_t\mathbf{A }_t)\mathbf{C }_{t,\mathbf{zz} } ,\, \end{aligned}$$
(4)

where \(\mathbf{C }_{t,\mathbf{zz} } \) is the covariance at the previous time step.

The conditional mean and covariance of the new observations are propagated to the simulated realizations of the model. By adding new observations, the covariance is reduced and the stationary assumption may no longer be satisfied. This is estimated in Eq. (4), where the covariance now depends on the previous time covariance (\(\mathbf{C }_{t,\mathbf{zz} }\)) and the Kalman gain (\(\mathbf{W }_t\)). The sequential updating process reduces the variability in the covariance of the initial random field represented after several iterations. This can provoke ensemble collapse, which happens when the same set of realizations is used to update the model during a certain time step (Leeuwenburgh et al. 2005; Baker 2007; Petrie and Dance 2010), and the evolved background error covariances are systematically underestimated after each assimilation cycle (Stroud et al. 2007).

3 Implementation Details for Underground Mining

This section describes the information that can be obtained during each operational phase in an underground extraction process.

Figure 1 shows the considered SMU extraction process in the mining area. The first step of the extraction process is in the pre-blasting phase, when the SMU to be extracted is exposed and has to be characterized. This is shown in Fig. 1(a), where Fig. 1(a1) is the sensor that characterizes the mine face, Fig. 1(a2). Channel samples [Fig. 1(a3)] may also be obtained out of the mine face during grade control. The composition of the element grade and the vein thickness within the mine face can be obtained based on different sensors imaging the exposed face, which is considered point-wise information. Therefore, there is no change in support, since the exact locations of the data points are known, and the observation support matches the support of the resource model considered.

To improve the prediction of the SMUs grade using the updating approach, the information is fed back into the grade control model before blasting in real time. The second step of the extraction process chain consists of blasting the SMU [Fig. 1(b)] that has been selected in the extraction sequence. Ore is drilled, blasted [Fig. 1(b6)], loaded and transported. The transport may be by means of a conveyor belt [Fig. 1(c8)], which leads to a primary crusher and subsequently into a storage bin. During this phase, the material can be characterized by cross-belt sensors [Fig. 1(c7)] installed along the conveyor line. The grade content is measured and represents a mean value of the whole SMU extracted. The grade content prediction of extracted SMU is usually performed as the mean value of the grid nodes within the SMU in the grade control model and corresponds to the operation modelled in Eq. (8). This corresponds to Fig. 1(b) and (c). When mining thin veins, the material scanned includes dilution from waste rock. The diluted grade and mineral content results from the ore/waste ratio within the blasted SMU and the grade or mineral content within the vein.

Fig. 1
figure 1

General extraction process where a represents the mine face exposed, b represents the scheduled blasted SMU, c represents the extraction and logistic of the blasted material and d shows the new mine face exposed

Full size image

3.1 Sensors

Two different types of sensors are considered: sensors characterizing the exposed mine face and sensors installed in the conveyor belt that characterize the blasted material. The sensors that characterize the mine face provide georeferenced information within the area that is being characterized and enable an accurate estimation of the vein thickness and thus the volume of ore to be extracted.

Desta and Buxton (2017) focused on the application of sensor technologies for raw material characterization in underground mining settings. Their application is made on the same test case as considered here. However, after blasting, the material that is being characterized on the conveyor belt cannot be georeferenced to a higher resolution than a SMU support. Thus, fuzziness in material tracking due to muck pile management and loading is assumed for the sensor data on the conveyor belt. To accommodate for this, different errors have been associated to the different sensors, accounting for the accuracy of the measurements.

Note that for this study, without loss of generality, the blasting of the SMU is assumed to be performed perfectly (as planned). In reality, information about over- or under-breakage will be considered during the operation using three-dimensional laser profiling data from the blasted SMU. These profiling sensors provide information on the surface of the SMU before and after blasting (Desta and Buxton 2017). With this information, the volume of SMU extracted is estimated.

3.2 Methodology for Underground Mine Settings

The underground mineral deposit is characterized by different attributes. In the present work these are the arsenic grade and thickness of the vein over each grid location. The representation of any of the attributes in this paper has been defined by the univariate spatial random function previously described \(\mathbf{Z }(\cdot )\). Each variable may be represented in different supports over the whole mineral extraction process.

Figure 2 represents a flowchart of the data assimilation process for grade control model updating. This is divided into two main supports. The grid support is a fine scale, one that matches with the initially obtained data support as chip samples, channel sampling, boreholes or sensor measurements on the mine face. The resource model is simulated on a fine-scaled grid comprising n grid nodes. The associated SMU model is obtained through re-blocking \(\mathbf{Z }(\mathbf{x }^F)\) into \(\mathbf{B }(\mathbf{x }^B)\). The phases where the assimilation algorithm is implemented are mine face mapping and bulk sampling. The first one is the phase where the sensor characterizes the mine face, as in Fig. 1(1). The bulk sampling corresponds to the phase where the sensors analyse the information transported by the conveyor belt, as shown in Fig. 1(8). In real applications, the analysis of this material may have an influence on blending or time delay. In the present paper, the sequence of sampling is simplified, as the data are available at every time step, and blending has no influence on the process. Thus, the aim of this paper is to provide an analysis of updating the model with information of different supports.

Fig. 2
figure 2

Flowchart describing the data assimilation process for control model updating

Full size image

Moreover, the two simulated models are informed respectively on two sets of points: \(S^F_0 = \{\mathbf{x }^F_1, \dots , \mathbf{x }^F_{n}\}\), where \(\mathbf{x }^F\) refers to any grid contained in \(S^F_0\); and \(S^B_0 = \{\mathbf{x }^B_1, \dots , \mathbf{x }^B_{m}\}\), where \(\mathbf{x }^B\) is any node contained in \(S^B_0\).

As discussed in Wikle and Berliner (2007), Amezcua and Van Leeuwen (2014) and Benndorf (2015), the sequential updating methodology is optimal when all variables involved are Gaussian. For most variables, it is necessary to perform a mapping to a Gaussian space prior to updating. The most commonly used transformation is the normal score transform, where a transformation to normal scores is achieved via quantile matching (Simon and Bertino 2009; Carrassi et al. 2018).

Two different transformations (\(\Phi \)) have been obtained for each support as

$$\begin{aligned} \Phi _{n}(\mathbf{U }_{\alpha }(\mathbf{x }_{\alpha }))&= \mathbf{Z }_{0}(\mathbf{x }_{\alpha }); \quad&\mathbf{x }_\alpha \in S_{0}^{F} ,\, \end{aligned}$$
(5)
$$\begin{aligned} {\Phi }_{m}(\mathbf{Q }^i(\mathbf{x }^{B}))&= \mathbf{B }^i(\mathbf{x }^{B})+\mathbf{E }_k^i(\mathbf{x }^{B}) ; \quad&\mathbf{x }^{B} \in S_{0}^{B} ,\, \end{aligned}$$
(6)

where i denotes each realization of the whole ensemble of each simulated random field. The Gaussian-transformed random fields for the grid-support and the SMU-support simulated models are represented by \(\mathbf{U }_{\alpha }(\mathbf{x }_{\alpha })\) and \(\mathbf{Q }^i(\mathbf{x }^{B})\). Equation (5) shows the normal score transformation for the initial data point support (\(\Phi _\text {n}\)).

This was obtained from the initial information time 0 at the data locations \(\mathbf{x }_\alpha \). After the application of this transformation, a simulation of the random field in a grid support is implemented by sequential Gaussian simulation.

Equation (6) shows the normal scores of the transformation for a SMU support (\({\Phi }_{m}\)). This was obtained after adding a random noise \(\mathbf{E }_k^i(\mathbf{x }^{B})\) to each SMU of each realization of the re-blocked model \(\mathbf{B }^i(\mathbf{x }^{B})\). This error is associated with the bulk sampling instrumental measurement, and it is simulated assuming a mean value equal to zero and uncorrelated covariance.

Figure 2 shows two different sources of observations: those made by mine face mapping and those made by bulk sampling. Mine face mapping obtains grid support observations, while bulk sampling obtains SMU support observations.

The mining process proposed here makes a two-step characterization of the SMUs studied. A time step index \(t+\frac{1}{2}\) refers to updates based on mine face mapping, and a time step index \(t+1\) refers to updates based on bulk information. Both measurements always map the same SMU, but each case refers to a different support. The set of grid nodes sampled for each time step is defined by \(S^F_{t+\frac{1}{2}}\) and \(S^B_{t+1}\). These are subsets from \(S_{0}^{F}\) and \(S_{0}^{B}\), respectively. By this assumption, one SMU update consists of two measurements and two updates. The reason that this double update is implemented over one SMU is to update the information regarding thickness and grade content before the blasting. A sequential extension of this assumption to the assimilation process is implemented, alternating the two types of observations.

With the proposed methodology, there exist two different data assimilation windows, depending on the observations that are being obtained: the window grid-grid that updates based on grid observations, and the window SMU-grid that updates based on SMU observations. Nevertheless, both windows update the grid model.

One of the main differences between the two assimilation windows is the forward operator that maps the observations of each forward model

$$\begin{aligned} {\mathcal {A}}_{t+\frac{1}{2},k}\big (\big \{\mathbf{Z }_{t}^i(\mathbf{x }^{F}) | \mathbf{x }^{F} \in S_{0}^{F} \big \}\big )&= \left\{ \mathbf{Z }_{t}^i(\mathbf{x }^{F}) |\mathbf{x }^{F} \in S^F_{t+\frac{1}{2}} \right\}&= \mathbf{Z }_{t+\frac{1}{2}}^i\big (\mathbf{x }_k^{F}\big ) ,\, \end{aligned}$$
(7)
$$\begin{aligned} {\mathcal {A}}_{t+1,b}\big (\big \{\mathbf{Z }_{t}^i(\mathbf{x }^{F}) | \mathbf{x }^{F} \in S_{0}^{F} \big \}\big )&= \big (\#S_{t+1}^B\big )^{-1} \sum _{\mathbf{x }^{F} \in S_{t+1}^B} \mathbf{Z }_t\big (\mathbf{x }^{F}\big )&= \mathbf{B }_{t+1}^i\big (\mathbf{x }_b^{B}\big ) .\, \end{aligned}$$
(8)

Equation (7) is the forward operator of the mine face observations. This operator maps the observations onto the resource model with a correspondence in the support between the observations and the background model. Equation (8) is the forward operator for observations with SMU support. This provides an average value of the grids contained in the volume of the SMU observed from the resource model simulated. Both cases are mapped from the observation locations at each time defined by the sets \(S_{t+\frac{1}{2}}^F\) and \(S_{t+1}^B\). These are the subsets of the grid and SMU sets of locations. These constitute the locations of the observations made in each time step. As stated above, an optimal implementation of the algorithm requires linearity of the observation operator, and observations must be ensured. In the present work, Eqs. (7) and (8) are assumed to be linear. However, these can also be non-linear (e.g., when the property is being measured over or under a threshold or when the measurement error is non-Gaussian). Implementing the sequential updating algorithm under suboptimal conditions can still be satisfactory in a non-linear, non-Gaussian and high-dimensional setting (Wikle and Berliner 2007; Amezcua and Van Leeuwen 2014; Carrassi et al. 2018).

After applying the forward operator and mapping the predicted observations into the forward model, these are transformed by the defined transformations as

$$\begin{aligned} \mathbf{U }^i_{t+\frac{1}{2}}\big (\mathbf{x }_k^F\big )&= {\Phi }^{-1}_{n}\big (\mathbf{Z }^i_{t+\frac{1}{2}}\big (\mathbf{x }_k^F\big )\big );&\mathbf{U }_{t+\frac{1}{2}}\big (\mathbf{x }_k^F\big ) \in {\mathbb {R}}^{b \times I} ,\, \end{aligned}$$
(9)
$$\begin{aligned} \mathbf{Y }^i_{t+1}\big (\mathbf{x }_b^B\big )&= {\Phi }^{-1}_m\big (\mathbf{B }^i_{t+1}\big (\mathbf{x }_b^B\big )\big );&\mathbf{Y }_{t+1}\big (\mathbf{x }_b^B\big )\in {\mathbb {R}}^{k \times I} ,\, \end{aligned}$$
(10)

where k and b are the number of observations for each support.

In this way, the predicted observations are mapped into the Gaussian space. In order to compute the estimated covariances, the mine face mapping must also account for the instrumental measurement error. This is denoted by \(\mathbf{E }_b(\mathbf{x }_k^F)\), and is obtained by randomly drawn values from a normal distribution with mean zero. Since measurements are independent of location, this covariance is assumed to be diagonal. The transformation of the grids is also applied, obtaining a new variable \(\mathbf{T }_{t+1}(\mathbf{x }_k^F)\) as

$$\begin{aligned} \mathbf{T }_{t+1}\big (\mathbf{x }_k^F\big )&= {\Phi }^{-1}_{n}\big (\mathbf{Z }_{t+1}\big (\mathbf{x }_k^F\big ) + \mathbf{E }_b\big (\mathbf{x }_k^F\big )\big );&\mathbf{T }_{t+1}\big (\mathbf{x }_k^F\big ) \in {\mathbb {R}}^{k \times I} .\, \end{aligned}$$
(11)

After transforming the grids and SMU nodes, the normal score transformation is applied to the actual observations \(\mathbf{d }_{t+\frac{1}{2}}\) and \(\mathbf{d }_{t+1}\) as

$$\begin{aligned} \mathbf{s }_{t+\frac{1}{2}}(\mathbf{x }_k)&= {\Phi }^{-1}_{n}\big (\mathbf{d }_{t+\frac{1}{2}}(\mathbf{x }_k)\big );&\mathbf{s }_{t+\frac{1}{2}}(\mathbf{x }_k) \in {\mathbb {R}}^{k \times I} ,\, \end{aligned}$$
(12)
$$\begin{aligned} \mathbf{s }_{t+1}(\mathbf{x }_b)&= {\Phi }^{-1}_{m}\big (\mathbf{d }_{t+1}(\mathbf{x }_b) \big ) ;&\mathbf{s }_{t+1}(\mathbf{x }_b) \in {\mathbb {R}}^{b \times I} ,\, \end{aligned}$$
(13)

where \(\mathbf{s }_{t+\frac{1}{2}}(\mathbf{x }_k) \) and \(\mathbf{s }_{t+1}(\mathbf{x }_b)\) are replicated into the matrix \(\mathbf{S }_{t+\frac{1}{2}}(\mathbf{x }_k) \) and \(\mathbf{S }_{t+1}(\mathbf{x }_b)\) as \(\mathbf{D }_{t+1}\) in the previous section.

After setting the different relations between the observations and the forward predicted model, the update methodology is implemented. The weights stated in Eq. (3) are approximated as

$$\begin{aligned} \mathbf{W }_{t+\frac{1}{2},k}&= \hat{\mathbf{C }}_{t+\frac{1}{2},\text {UU}} \hat{\mathbf{C }}_{t+\frac{1}{2},\text {TT}}^{-1};&\mathbf{W }_{t+\frac{1}{2},k} \in {\mathbb {R}}^{n\times k} ,\, \end{aligned}$$
(14)
$$\begin{aligned} \mathbf{W }_{t+1,b}&= \hat{\mathbf{C }}_{t,\text {UY}} \hat{\mathbf{C }}_{t,\text {QQ}}^{-1} ;&\mathbf{W }_{t+\frac{1}{2},b} \in {\mathbb {R}}^{n\times b} .\, \end{aligned}$$
(15)

Here, \(\hat{\mathbf{C }}\) are the approximated auto- and cross-covariances of the model and the predicted observations for different supports. However, the variables selected to approximate these covariances in the different assimilation windows are different.

For the grid-grid window, these covariances are approximated as

$$\begin{aligned} \begin{aligned} \hat{\mathbf{C }}_{t+\frac{1}{2},\mathbf{UU} }&= E\big [\mathbf{U }^i_{t+\frac{1}{2}}\big (\mathbf{x }^{F}\big ) - E\big [\mathbf{U }_{t+\frac{1}{2}}\big (\mathbf{x }^{F}\big )\big ]\big ][\mathbf{U }^i_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big ) - E\big [\mathbf{U }_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big )\big ]\big ]^{\mathrm{T}} ,\,\\ \hat{\mathbf{C }}_{t+\frac{1}{2},\mathbf{TT} }&= E\big [\mathbf{T }^i_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big ) - E\big [\mathbf{T }_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big )\big ]\big ]\big [\mathbf{T }^i_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big ) - E\big [\mathbf{T }_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big )\big ]\big ]^{\mathrm{T}} ,\, \end{aligned} \end{aligned}$$
(16)

where \(\hat{\mathbf{C }}_{t+\frac{1}{2},\mathbf{UU} } \in {\mathbb {R}}^{n \times k}\) and \(\hat{\mathbf{C }}_{t+\frac{1}{2},\mathbf{TT} } \in {\mathbb {R}}^{k \times k}\), and i denotes each realization of the whole ensemble of each simulated random field.

On the other hand, for the SMU-grid window, these covariances are approximated as

$$\begin{aligned} \begin{aligned} \hat{\mathbf{C }}_{t+1,\mathbf{UY} }&= E\big [\mathbf{U }^i_{t+1}\big (\mathbf{x }^{F}\big ) - E\big [\mathbf{U }_{t+1}\big (\mathbf{x }^{F}\big )\big ]\big ]\big [\mathbf{Y }^i_{t+1}\big (\mathbf{x }_b^{B}\big ) - E\big [\mathbf{Y }_{t+1}\big (\mathbf{x }_b^{B}\big )\big ]\big ]^{\mathrm{T}} ,\,\\ \hat{\mathbf{C }}_{t+1,\mathbf{QQ} }&= E\big [\mathbf{Q }^i_{t+1}\big (\mathbf{x }_b^{B}\big ) - E\big [\mathbf{Q }_{t+1}\big (\mathbf{x }_b^{B}\big )\big ]\big ] \big [\mathbf{Q }^i_{t+1}\big (\mathbf{x }_b^{B}\big ) - E\big [\mathbf{Q }_{t+1}\big (\mathbf{x }_b^{B}\big )\big ]\big ]^{\mathrm{T}},\, \end{aligned} \end{aligned}$$
(17)

where \(\hat{\mathbf{C }}_{t+1,\mathbf{UY} } \in {\mathbb {R}}^{n \times b}\) and \(\hat{\mathbf{C }}_{t+1,\mathbf{QQ} } \in {\mathbb {R}}^{b \times b}\).

After obtaining the Kalman gains for each window, the updating algorithm is implemented

$$\begin{aligned} \mathbf{U }_{t+\frac{1}{2}}\big (\mathbf{x }^{F}\big )&= \mathbf{U }_{{t}}\big (\mathbf{x }^{F}\big ) + \mathbf{W }_{t+\frac{1}{2},k}\big [\mathbf{S }_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big ) - \mathbf{U }_{t+\frac{1}{2}}\big (\mathbf{x }_k^{F}\big )\big ],\, \end{aligned}$$
(18)
$$\begin{aligned} \mathbf{U }_{t+1}\big (\mathbf{x }^{F}\big )&= \mathbf{U }_{{t}}\big (\mathbf{x }^{F}\big ) + \mathbf{W }_{t+1,b}\big [\mathbf{S }_{t+1}\big (\mathbf{x }_b^{B}\big ) - \mathbf{Y }_{t+1}\big (\mathbf{x }_b^{B}\big )\big ] .\, \end{aligned}$$
(19)

So far, the implementation details for the method for underground mining settings have been presented. However, the results must further be evaluated in order to quantify the improvement in the SMU prediction. For that purpose, the mean absolute error (MAE) is used. This is done for the grid and the SMU models (\(\text {MAE}^{F}\) and \(\text {MAE}^{B}\)). It is defined as the mean of the absolute value of the difference between the true value (\(\mathbf{Z }^*(\mathbf{x }^{F})\)) and the model-based prediction (\(\mathbf{Z }_{t+1}(\mathbf{x }^{F})\)) for the grid model and for the SMU model (\(\mathbf{B }^*(\mathbf{x }^{B})\) and \(\mathbf{B }_{t+1}(\mathbf{x }^{B})\)).

$$\begin{aligned} \text {MAE}^{F}&= (I-1)^{-1} \sum _{i=1}^{I} \Vert \mathbf{Z }_{t+1}^i(\mathbf{x }^{F}) - \mathbf{Z }^*(\mathbf{x }^{F}) \Vert ;&\quad&\mathbf{x }^{F} \in S_0^{F} ,\, \end{aligned}$$
(20)
$$\begin{aligned} \text {MAE}^{B}&= (I-1)^{-1} \sum _{i=1}^{I} \Vert \mathbf{B }_{t+1}^i(\mathbf{x }^{B}) - \mathbf{B }^*(\mathbf{x }^{B}) \Vert ;&\quad&\mathbf{x }^{B} \in S_0^{B} .\, \end{aligned}$$
(21)

3.3 Diluted Grade

In order to incorporate sensor measurements from diluted ore into the prediction model, an approach is implemented that considers information related to different grades and vein thickness in the grade control model. After blasting, the material is carried on conveyor belts, where sensors are able to measure diluted grade. Equation (22) provides the calculated value for the recovery model. This combines the final percent dilution related to the resource model and the mining reserve estimate related to the grade control model

$$\begin{aligned} \text {Diluted Grade} = \frac{V_O \times O\%}{V_O + V_B } .\, \end{aligned}$$
(22)

Volume \(V_B\) of the SMU as blasted can be obtained from a laser profiling of the mining SMU, and \(V_O\) is the volume of the vein contained in the extracted SMU. In the case of perfect blasting, this volume equals the area of the SMU multiplied by the mean thickness value. This information is available from the already updated model based on the face image. Thus, the diluted grade sensed on the conveyor belt provides indirect information about the SMU grade.

4 Experiment

In this section, a description of the experiments performed is provided to analyse the algorithm implemented. The following aspects will be discussed: (i) the Reiche-Zeche test case and the generated ground truth model, and (ii) the available ground truth model and sampling scenarios, with two different sampling regimes providing operational scenarios resulting in grade control data; (iii) conditional simulation of prior grade control models for each of the scenarios and an unconditional simulation scenario, (iv) generation of artificial sensor data representing a predefined extraction process by means of the ground truth models, (v) updating of prior models to generate a posterior grade control model utilizing the developed algorithm, and (vi) comparison of posterior models to the ground truth model and assessment of method performance.

4.1 Reiche-Zeche Test Case

The Reiche-Zeche mine is located in Freiberg, Germany. The Freiberg region is the oldest mining district in the eastern part of Erzgebirge. It was mined for Ag and for Cu, Pb and As. The Freiberg ore vein network is characterized by two (NNE-SSW to N-S and E-W to ENE-WSW) shear systems, and spatially associated fissure veins. In general, ores in the Freiberg mining district are associated with a system of dykes. The Freiberg polymetallic sulphide deposit was formed by two hydrothermal mineralization events of late Variscan and post-Variscan age. The central part of the mine comprises predominantly late Variscan mineralization and is rich in sulphur, iron, lead, zinc and copper. The ore minerals typical for this area are arsenopyrite, pyrite, sphalerite, chalcopyrite and agglomerates of galena–quartz–carbonate gangue (Stockmann et al. 2013).

The test site is at the Wilhelm Nord part of the Wilhelm Stehender orebody, located in the third level of the mine, at 200 m depth. The solid part of the orebody is a continuous vein that has a dip of around 50\(^{\circ }\) ENE-WSW and thickness that varies from 0 to 1.4 m.

In order to implement and benchmark the data assimilation process within a fully known environment, a synthetic data set based on real information was created based on the result of a sampling campaign carried out between 1985 and 1989.

4.1.1 Available Database

The main data set that was used consists of 114 channel samples across the mineralized zone distributed along two different drifts. The samples were taken at 1-m intervals and analysed in a laboratory. Each channel sample provides information about both the grade and the thickness of the vein. This geological feature was modelled as the resource model. Thickness is an indicator of the ore content and gangue for each SMU. The ore grade represents the available grade control data.

The grade that is considered for demonstrating the assimilation process is arsenic (As), which presents a maximum value of 18% (Fig. 3). The remaining grades are of no significance for this study and are considered waste. The size of the resource block is 60 m \(\times \) 60 m. The panel being extracted has a size of 30 m \(\times \) 30 m and corresponds to the lower left corner of the SMU.

4.1.2 Ground Truth Model Estimation

The ground truth model generated from the available data set was modelled as a realization by conditional sequential Gaussian simulation.

Fig. 3
figure 3

The histograms show the distribution of arsenic (left) and thickness (right) for the 114 samples at \(\text {x}_\alpha \)

Full size image

To provide a ground truth model capturing both the vein thickness and grade content distribution, the attributes grade percentage content, elevation of the hanging wall and the vein thickness were simulated. This study was performed in two dimensions. The ground truth model was generated on a grid with dimensions of 0.2 m \(\times \) 0.2 m for the considered resource block of 60 m \(\times \) 60 m.

4.1.3 Sampling

To analyse the performance of the sequential updating algorithm applied to different grade control strategies, different scenarios constrained to different levels of information were evaluated.

  1. 1.

    Sampling 1: The first scenario mimics a sampling campaign, where grade control data are available at the contour (mine face exposed) of the block considered (Fig. 4). In a two-dimensional case, these data are lines in the contour of the whole block and are represented in black as shown in Fig. 4. By means of available sensor technology, it is possible to characterize the mine face exposed in the SMU. The prior model generated based on this sampling campaign is denoted as a mine face scenario.

  2. 2.

    Sampling 2: The second scenario combines information from the mine face scenario with data from three boreholes that are oriented towards the mining direction, and explores the centre of the mining panel. Information about the mineral content is considered to be available on a 0.2 m borehole interval. The right-hand side of Fig. 4 shows the information that constrains this prior model as black lines. This scenario is referred to as a borehole scenario.

Fig. 4
figure 4

Representation of the information obtained in the sampling campaigns for a whole panel. The right part represents the mine face sample information and the left part the borehole sample information

Full size image

4.1.4 Generation of the Prior Model Conditional Simulation of the Key Variables

Having a ground truth model and defining the sampling scenarios considered, three different prior grade control models are simulated. The first model is unconditional on any information sampled, meaning that the random field is simulated with only the information of the covariance function. The second model is conditional on the mine face data, while for the third model the borehole data provide the conditioning information.

The ground truth model and the models based on the three scenarios are simulated using the same covariance function, derived from sampling the exposed mining face. Sequential Gaussian simulation at grid size scale is performed for the arsenic content and thickness. A total of 300 realizations are generated for each of the three desired models. Figure 5 shows the absolute error of the mean (mean value of all realizations versus reality), the mean grade and standard deviation of the of the 300 realizations for the arsenic grade. This is depicted for the three sampling scenarios.

Fig. 5
figure 5

Map in point scale (each cell is 0.2 m \(\times \) 0.2 m) of As grade for error, expected value and standard deviation (SD). From up to down, unconditional, conditioned on mine face information and conditioned on borehole information models

Full size image

The grid size settings considered for the conditional simulation are kept the same as for the ground truth model.

4.2 Algorithm Settings

The mining scenario was set for a SMU of dimensions 2 m \(\times \) 2 m \(\times \) 2 m. The algorithm was implemented in a two-dimensional test case. However, due to the vein thickness considered, this gives a third dimension to the assimilation procedure. The thickness presents a maximum of 1.4 m, which means that the vein is fully contained within each SMU. The block simulated has an area of 60 m \(\times \) 60 m. A re-blocking was performed as indicated in Eq. (8) on the simulated prior models but also in the ground truth model. For purpose of illustration, the area was divided into four quarters (panels), and the assimilation sequence was performed in the lower left panel.

The algorithm was applied to a total of 225 SMUs, representing 15 drifts with 15 SMUs per drift. The updating sequence simulated for the first drift was defined from the lower left SMU to the lower right SMU. This direction is also applied to all subsequent drifts. The final drift starts in the upper left corner and ends in the upper right corner. To update locally, a localization function was implemented to exponentially decrease the correlations with the distance to the extracted SMU (Leeuwenburgh et al. 2005; Baker 2007; Petrie and Dance 2010; Carrassi et al. 2018). This function has a value of 1 in the centre and a radius of 8 m from the centre of the SMUs being updated. Beyond this range, grid nodes are not updated. The application of this function guarantees no abrupt change in domain between the updated (with value of 1) and non-updated (with value of 0) areas.

The synthetic experiment mimics reality by acquiring information from the production process via the observations provided by sensor technologies. The updating process accounts for sensor errors. This is shown in Eq. (3), represented by the term R. The relative error assigned to the mine face characterization is 0.05, while the relative error assigned to the conveyor belt sensor analysis is 0.1. The loss of material during loading and blasting was considered by giving a higher error to the conveyor belt sensor analysis.

4.3 Results of the Experiment

In this section, the results of the previously defined experiments are presented. Different visualizations characterize the statistical moments and show the improvements achieved by the data assimilation implementation. First, two-dimensional overview maps are used to depict the model standard deviation, the true error (residual) and the mean value of the set of realizations. Box plots represent the evolution of the distribution of each SMU for every drift. Absolute mean error analysis is applied to quantify the improvement in the data assimilation process performed.

While Fig. 5 shows the different initial scenarios of the arsenic grade in a detailed grid of 0.2 m \(\times \) 0.2 m, Fig. 6 shows the mean and error for all realizations of arsenic grade in a point and SMU support of 2 m \(\times \) 2 m. The first two maps on the left represent the arsenic grade maps for the unconditional case after being updated at times \(t=21\) and \(t=141\). The first two maps on the right correspond to the error associated with the corresponding cases. The four remaining maps represent the same in a SMU support. When the information observed during the production process is incorporated into the model, areas with extreme values are appreciable. The maps change notably by incorporating the new information into the model. The resolution of areas next to the assimilated SMU improves when the observed information is incorporated.

Fig. 6
figure 6

Map in point and SMU scale of unconditional arsenic grade in the initial space. The black square indicates which SMU is being assimilated at time \(t=21\) and \(t=141\)

Full size image

Figure 7 depicts the histograms for the 21st and 141st SMUs at different updated times for the arsenic variable in the Gaussian-transformed space. The underlying histogram represents the initial ensemble. The variance in the distribution of the updated ensembles decreases when the SMU is closer to the observation. Additionally, the mean value of the distribution is closer to the observed value.

Fig. 7
figure 7

Histograms of the realizations of the 21st and 141st SMUs assimilated in the Gaussian-transformed space at times previous to be extracted

Full size image

Figure 8 represents the box plots of SMU mean of arsenic grade values for each SMU in the ninth drift unconditioned on sampling data. The red point plotted in each SMU corresponds to the true SMU arsenic grades. The upper plot shows the model before the assimilation of any SMU, and the lower plot shows the prediction after assimilating SMU 140. There is a marked reduction in the variability of the arsenic SMU means after assimilation. The SMU means for SMUs 136 to 141 are now very close to their respective true means.

Fig. 8
figure 8

Box plot representation of the SMU before assimilation (above) and after assimilation of SMU 141 (below). Both correspond to the ninth drift and for the unconditional scenario of the arsenic variable in the raw space

Full size image

Figure 9 shows the analogous representation from Fig. 8 for the borehole scenario. In this case, one can observe SMU 143, which is conditioned on the borehole intersecting that SMU. SMU 143 appears to have a smaller variance compared to those in Fig. 8. Figure 9 also depicts a lower variance for all the distributions at the initial time. After updating at time \(t=143\), similar distributions are displayed for SMUs from 144 to 150 in both scenarios. The initial information is no longer appreciable between the updated scenarios.

Fig. 9
figure 9

Box plot representation of SMU 140 before the assimilation (above) and after the assimilation (below). Both correspond to the ninth drift and for the borehole scenario of the arsenic variable in the raw space

Full size image

Figure 10 represents the implementation of Eq. (10). This is a SMU-level map that combines the vein thickness volume and the arsenic grade. This representation is an evaluation of the dilution effect at each SMU.

Fig. 10
figure 10

SMU-map of the diluted grade representation at time \(t=21\) and \(t=141\)

Full size image

Figure 11 shows the MAE of a whole panel based on one further SMU to be mined. This is shown before and after the assimilation process. From left to right, the three different scenarios are considered: the borehole conditioned scenario, the mine face conditioned scenario and the unconditional simulation scenario.

The analysis focused on the non-assimilated models shows that the error for the unconditional model is higher than that for the models conditional on borehole and mine face information. One can also see that the error decreases with the amount of conditioning information available. Therefore, the borehole scenario has a smaller error than the mine face scenario.

On the other hand, the analysis focused on the models after applying the assimilation algorithms in Fig. 11 shows minimal differences in mean absolute errors. The simulated models updated by observations obtained during the production phase have lower error. Also, the error decreases to almost the same levels. This is because the amount of information assimilated from various data sources is similar. The information assimilated reduces the uncertainty presented by the initial simulated models even, though the support for the observations assimilated and the model are different.

As one can see, the uncertainties of all the SMUs of a panel for the simulated assimilated models are lower than those for the non-assimilated counterparts. This result validates the applicability of the assimilation method for the panel. Conditioning the simulated models on the production information overwrites initial conditioning information for each degree (borehole, mine face, unconditional) that was previously used to generate the resource model.

Fig. 11
figure 11

MAE for each of the three different simulations before and after the updating process evaluated in all the SMUs that are one SMU further from the extracted SMU

Full size image

Figure 12 presents the MAE [as in Eqs. (20) and (21)] of each drift for the whole panel. Each X step represents different drifts. The SMUs covered by this drift are represented on the axes. The upper plot represents the mean value of all the SMUs contained by each drift that will be mined next. This means that the averaged drift has been updated with updated SMUs by the contiguous SMU.

The middle plot shows the MAE per drift for all the updated SMUs considering the information from SMUs two steps contiguous. The lower plot shows the error for all the averaged drift updated with SMUs that are three steps contiguous. Figure 12 represents the uncertainty level for the decision-making in SMU extraction and processing.

In order to put the MAE in perspective, each graph from Fig. 12 also shows the error for the three different scenarios not considering any assimilation process but showing the corresponding conditioning levels. The error strongly decreases between pairs of models before and after the assimilation. The difference is more abrupt for models that have less conditioning information. One aspect to highlight is that for most ranges evaluated, each model derived from the assimilation has a lower error per drift than the model with a higher amount of conditioned information (borehole scenario). This is applicable to all the models except the first drift of the unconditional model against the borehole and mine face scenario. This indicates that the influence of the information assimilated is greater during extraction relative to the a priori available conditioning information based on traditional grade control.

Fig. 12
figure 12

Representation of the MAE per drift for the models applying the assimilation algorithm and without the assimilation algorithm. The first drift is covered from SMUs 1 to 15

Full size image

5 Conclusions

Ensemble-based sequential model updating has been implemented for real-time updating of an underground grade control model. This work provides a case study in a synthetic environment based on real information sampled from the Reiche-Zeche underground mine. The methodology was implemented for geochemical variables including arsenic grade and thickness of the ore vein, also considering the dilution during the extraction process.

An initial ground truth model was developed based on the exploration information. The model was sampled for creating three scenarios with different conditioning information levels prior to extraction. Based on this initial information, three initial grade control models were generated, employing a sequential Gaussian simulation. The information that feeds the mining production process was simulated considering different kinds of sensors suitable for characterizing the material in different mine phases, namely before blasting at the mining face and during transportation after the blending phase.

The algorithm updates a grid support model based on grid information obtained from the mine face before blasting, and SMU support information that is obtained after blasting the SMU unit within a drift. In order to ensure an optimal condition for the Gaussian-based updating algorithm, data were mapped onto a Gaussian space by an adjusted Gaussian anamorphosis strategy. The updating procedure considers the sample covariance of the observations, the sensor measurement error and the model uncertainty captured in the simulated model. The change in support between the observations and the simulated model was mapped by the forward simulator, which thus ensures that the estimated sample covariance accounts for the change in support (over a constant density of the material). This study shows the importance of updating the resource model by SMU support and SMU face characterization when the initial models were conditioned on different levels of information. This was analysed for 15 different drifts of a panel and for the whole panel. The results also show that the updated models are more informed than those that are not updated but are conditioned on the exploration borehole and mine face information. Analysis of the updated SMU errors at one, two and three SMUs from the currently extracted SMU shows that the errors decrease to comparable levels after assimilation, irrespective of the conditional sampling information on which the model was conditioned. Also, the error decrement is higher for the unconditional model and the updated unconditional model than for those initial models conditioned by a high level of information from borehole and mine face characterization.

In conclusion, the models obtained by implementing the updating algorithm minimize the influence of the exploration information at a short-term decision-making level, such as SMU, drift or panel support. One may consider revisiting the density of grade control sampling within a panel when integrating online production data. Future studies will focus on the reconciliation of the information in a multivariate case considering the relationship among mineral attributes as compositional variables.