Speeding up Conditional Simulation: Using Sequential Gaussian Simulation with Residual Substitution

...

data values and every realisation at the data locations:

where is the k-th non-conditional realisation of the Gaussian random field (x). Estimation of the residuals over the domain of interest by simple kriging (sk), using the covariance model . Because the weigths are the same in all the realisations, the residual estimates are simultaneously obtained in a single kriging run. Addition of the residual estimates and the corresponding non-conditional realisation to generate the conditional realisation over the domain:

Back-transformation of the conditional Gaussian realisations to the original values.

The non-conditional Gaussian simulation step is detailed below:

Define a single visiting sequence, which can be achieved using a low discrepancy sequence [7], a regular sequence or just a random sequence with a multiple grid approach [13]. With the LU algorithm, generate several () realisations for the first thousand nodes of the visiting sequence. These simulated values, which will be used as conditioning data for the subsequent nodes, account for the covariance model and are by construction independent from one realisation to another. Continue with the sequential approach and generate non-conditional realisations in each visited node :

where stands for the vector of simple kriging estimates for the realisations given the previously simulated nodes, for the simple kriging standard deviation, and for an independent Gaussian random vector. This way, all the realisations can be generated simultaneously.

As an example, for a domain of 100 × 100 nodes, one hundred non-conditional realisations are generated using an isotropic spherical variogram of range 20. Figure 1 (left) presents the probability intervals of the simulated values at each node as a function of the sequence order: it is remarkable that the intervals are almost the same for all the nodes. To give another look at this feature, the conditional variance is calculated at each node and plotted as a function of the sequence order (Figure 1, right) without exhibiting any trend or artificial pattern.

Figure 1: Probability intervals (left) and conditional variances (right) as a function of the visiting order

APPLICATION TO A MINING DATASET

Presentation of the Case Study

The area under study is part of the Río Blanco – Los Bronces porphyry copper deposit [14], a breccia complex located in the Chilean central Andes. A set of 2376 diamond drill hole samples, located in a 400 × 600 × 130 m3 volume, are available with information on total copper grades. Figure 2 presents the available data coloured by copper content.

Figure 2: Location map of copper grade data

Simulation Approaches

Three different algorithms are compared by simulating the copper grades over a 2D regular grid of 390 × 600 nodes with a 1 × 1 m spacing: Sequential Gaussian simulation (SGSIM): this approach is performed using the SGSIM program of the GSLIB package [3], using multiple grids and migration of data to nodes. Sequential Gaussian simulation with residual substitution (SGSIM-RS): a single path is used with multiple grids and a uniformly random ordering of the nodes of each grid. Turning bands simulation (TBSIM): this approach is performed in the ISATIS software, using 1000 turning lines.

The same search radius (250 m) and number of conditioning data (16) are considered for each method. Simple kriging is used for non-conditional simulation and conditioning step in SGSIM-RS, and in the conditioning steps of TBSIM and SGSIM. There is no loss of data in the migration to the nodes in SGSIM and SGSIM-RS, so the effective datasets are the same in each method. The same normal score transformation, back-transformation and isotropic variogram model (Table 1) are used in each method to avoid differences due to implementation.

Table 1: Isotropic variogram model for transformed copper grades

Structure

Range (m)

Sill contribution

Nugget

-

0.12

Spherical

112

0.7

Exponential

416

0.18

Comparison Between Methods

In the following subsections, the resulting realisations are compared in several ways. First of all, the visual inspection of several realisations of SGSIM-RS does not indicate the presence of any artefact or strange pattern. Even more, it is impossible to distinguish if the realisation comes from SGSIM, SGSIM-RS or TBSIM (Figure 3).

Figure 3: Examples of realisations

Basic statistics: Figure 4 (left) shows the distributions of the average copper grade per realisation for each algorithm. These distributions are close to each other, although that of SGSIM-RS presents slightly higher values. The variance per realisation is also indicated in Figure 4 (right): SGSIM-RS and SGSIM have similar distributions, while TBSIM shows a slightly wider range of variances.

Figure 4: Distribution of the average (left) and variance (right) of simulated copper grades per realisation (right)

Grade-tonnage curves: the tonnage and grades above a set of cutoff grades are calculated on each realisation. The expected curves are presented in Figure 5 (left), showing virtually the same values for all the algorithms. In counterpart, Figure 5 (right) presents the widths of the 80% confidence intervals for the grades and tonnages by cutoff. It is seen that the turning bands algorithm presents a higher variability in grades for almost every cutoff, whereas SGSIM-RS and SGSIM exhibit a similar behaviour for this measure.

Figure 5: Expected grade-tonnage curves (left) and widths of 80% confidence intervals for grade-tonnage curves (right)

Local expectation and uncertainty measures: Figure 6 (left) shows the distributions of the conditional expectation (mean of realisations) calculated at each node. The curves are almost identical for the three algorithms. The conditional variance distributions (Figure 6, right) are close for SGSIM and SGSIM-RS, whereas TBSIM shows a 5% of the nodes with higher values.

Figure 7 displays the conditional coefficient of variation. The high and low zones are located