- Email: [email protected]
Comparison of indicator kriging, conditional indicator simulation and multiple-point statistics used to model slate deposits
Available online at www.sciencedirect.com
Engineering Geology 98 (2008) 50 – 59 www.elsevier.com/locate/enggeo
Comparison of indicator kriging, conditional indicator simulation and multiple-point statistics used to model slate deposits F.G. Bastante a , C. Ordóñez a,⁎, J. Taboada a , J.M. Matías b a
Department of Natural Resources and Environmental Engineering, University of Vigo, Vigo, Spain b Department of Statistics, University of Vigo, Vigo, Spain
Received 17 September 2007; received in revised form 25 December 2007; accepted 26 January 2008 Available online 8 February 2008
Abstract The resources in an ornamental slate deposit can be estimated using geostatistical estimation techniques applied to information collected from drill cores. The result, however, is a smooth approximation that fails to take account of the natural variability in mineralization, which is fundamental to proper design and evaluation of the financial viability of a mining deposit. Geostatistical simulation techniques are more useful in this respect, as they reflect different realizations of the reality and better reflect natural mineral dispersion. In this work, we evaluate the resources of a slate quarry by comparing the results obtained using two geostatistical techniques—indicator kriging (ik) and sequential indicator simulation (sisim)—with the results obtained using the single normal equation simulation (snesim) technique based on multiple-point statistics (mps), analyzing their usefulness in evaluating the financial risk derive from uncertainty in regard to knowledge of the deposit. Our results indicate that although the multiple-point statistics approach produces models that are closer to reality than the models produced by the geostatistical techniques, the simulation relies in part on information obtained via indicator kriging. © 2008 Elsevier B.V. All rights reserved. Keywords: Geostatistics; Indicator kriging; Conditional simulation; Multiple-point statistics; Mining resource evaluation
1. Introduction The mining of ornamental rocks, and especially slate for roofing, has undergone spectacular development during the last 30 years in Spain, a fact which has placed this country first among the world's slate producers. This has brought about a considerable improvement in both productive processes and the machinery used, but there are a number of technical factors that have not yet been satisfactorily resolved. One example is the exploration and evaluation of slate deposits for mining purposes. Slate is a fine-grained, homogeneous, metamorphic rock derived from an original shale-type sedimentary rock composed of clay through low-grade regional metamorphism. It is mainly
⁎ Corresponding author. E.T.S.I.MINAS, Universidad de Vigo, Campus Universitario, 36310 Vigo, Spain. Tel.: +34 986814052; fax: +34 986812201. E-mail address: [email protected] (C. Ordóñez). 0013-7952/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.enggeo.2008.01.006
composed of quartz and muscovite or illite, often along with biotite, chlorite, hematite and pyrite, and less frequently, apatite, graphite, kaolin, magnetite, tourmaline, zircon and feldspar. As methods for slate quarrying, we can distinguish clearly between the more frequently utilized surface mining method, and the underground method that is gradually becoming more popular (Taboada et al., 2006a). Slate is commonly used as a roofing and flooring material and for electrical insulation. Geological investigation of a slate deposit is mainly based on continuous drill cores, since surface alterations impede surface sampling of fresh rock. The cores are used to study the sedimentary, metamorphic and structural characteristics that define the orebody. The orebody typically contains a mass of useful slate (which is defined, as we shall see subsequently, as the equivalent of mineral grade) and areas of spoil, the latter mainly a consequence of structural discontinuities in the orebody. In slate mining as in the rest of mining, however, all knowledge is inevitably linked to uncertainty, and this is represented as risk when a mining deposit is being assessed. In order to obtain the necessary decision-making
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
elements for each of the developmental phases of the mining project, the mining company has to quantitatively calculate a value for each parameter that intervenes in the mining evaluation along with the uncertainty associated with these values. This uncertainty has a significant bearing on the mining company's final decision as to whether to assume the risk of allocating human and financial resources to the development of a mining project. The problem of grade uncertainty and risk effects in open-pit design is addressed by Dimitrakopoulos et al. (2002), who use generalized sequential Gaussian simulation and direct block simulation techniques to model and evaluate the financial risk implied by a typical, disseminated, low-grade, epithermal, quartz breccia-type gold deposit. In this article, we tackle the same problem for slate deposits, although taking a somewhat different approach. We estimate useful slate reserves using indicator kriging and compare the results with those obtained using conditional simulation techniques and multiple-point statistics (Guardiano and Srivastava, 1993). 2. Roofing slate properties: Useful slate and saleable slate Only part of the slate in a quarry can be exploited as roofing slate, as slate for this purpose must have certain physical, chemical and aesthetic properties and must not be affected by structures that interrupt the cleavage direction. This slate is referred to as useful slate. From a general perspective, the usefulness of slate for roofing is determined from samples extracted from quarry fronts which are laboratory tested for compliance with the Spanish UNE standard 12326-1: 2005 (AENOR, 2005). These tests analyze slate samples in terms of a range of physical and chemical properties such as fissility, impermeability, carbonate content, colour, alterability, inclusions, homogeneity of grain size, and roughness. (See Taboada et al., 2006b) for a more detailed description of these properties). The test results combined with data from continuous drill cores enable experts to determine which parts of the drill cores represent useful slate. Thus, the expert divides the cores into 1-m long stretches and assigns them to Class 1 if classified as useful slate and Class 0 if classified as non-useful slate. The choice of 1-m lengths is justified by the fact that mining requires a minimum block thickness of 50 cm; consequently, as thickness is given by the normal direction to the cleavage plane that dips 60°, the support can be considered to be representative of the minimum block size in that direction. As will be demonstrated in Section 5.1 below, using indicator kriging will enable us to estimate both the percentage of stretches belonging to Class 1 within the blocks into which the quarry will be discretized and the percentage of useful slate contained in this support (in other words, the mineral grade for the block). In this way gross slate reserves can be estimated and mining can be designed optimally (Bastante et al., 2005). Of the useful slate indicated as such by the drill core sampling, only part can be recovered from the quarry in the form of blocks. Only a small proportion of this recovered slate is saleable, moreover, and only after the slate has undergone various mechanical operations in the slate transformation processing plant. This slate is referred to as finished or saleable slate.
51
3. Methodologies used to evaluate slate resources In our research, we used different techniques to evaluate useful slate reserves from the information provided by cores made in a slate quarry. Below we summarize the theory underlying each of the techniques and subsequently compare and discuss the results obtained for each method. 3.1. Indicator kriging The basic approach taken by predictive statistics is to model the uncertainty associated with an unsampled value, z(u), with u as the coordinate location vector, as a random variable (RV), Z(u), the probability distribution of which characterizes the uncertainty about z(u) (Deutsch and Journel, 1998). Categorical variables such as rock types can be effectively modelled by RVs. Let Z(u) be a categorical RV that can take either of 2 outcome values k = {k1, k2}. In the particular case to be analyzed here, z(u) is the variable stretch of a metre of slate centered on u; Z(u) is the random function that models z(u); and (k1, k2) corresponds to the categories (Class 1, Class 0) or useful/non-useful slate. The conditional probability density function is: f ðu; kjðnÞÞ ¼ProbfZ ðuÞ a category kjðnÞg (n) consisting of n neighbouring data values z(uα), α = 1,…, n. The essence of the indicator approach is the binomial coding of Z(u) into either 1 or 0 depending upon its relationship to a categorical class. By considering binary indicator transforms of Z(u) defined as: 1; if Z ðuÞ a k I ðu; k Þ ¼ 0; otherwise the probability that a category k prevails at u can be expressed as: EfI ðu; kjðnÞÞg ¼ Probfðu; k Þ ¼ 1jðnÞg ¼ f ðu; kjðnÞÞ: Applying the kriging algorithm to indicator data we obtain the least-squares estimate of its conditional expectation, f ⁎(u;k| (n). Note that the bivariate (two-point) distribution of 2 RV's Z (u1), Z(u2), appears as the non-centered covariance of the indicator variables: f ðu1 ; u2 ; k1 ; k2 Þ ¼ProbfZ ðu1 Þ ¼ k1 ; Z ðu2 Þ ¼ k2 g ¼ E fI ðu1 ; k1 ÞI ðu2 ; k2 Þg and the indicator variogram function measures how often two locations that are a vector h apart belong to different categories. Indicators are useful for characterizing the spatial variability of categorical variables as the ranges and shapes of the directional indicator semivariograms that reflect geometric patterns (Goovaerts, 1997). 3.2. Sequential indicator simulation Estimating reserves by means of kriging produces a smoothed representation of the reality, and this typically underestimates real data dispersion (Journel, 1974). Since geostatistics considers the
52
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
variable that is the object of the study to be a realization of a random function, it is possible to simulate an infinite number of realizations of random functions to represent the variability and spatial correlation characteristics observed in the experimental data (Journel and Huijbregts, 1978). The conditional simulation techniques reproduce the first two experimental moments (mean and covariance, or variogram) for the real data; in other words, they reproduce a measure of real dispersion in accordance with the information available, with the variogram acting as a two-point measure of geological heterogeneity or continuity. Therefore, whereas kriging provides the best unbiased estimator for each point, the conditional simulation produces values that reproduce real fluctuations. In the sequential indicator simulation technique, the simulated value for each point (1 or 0) is obtained by applying the Monte Carlo simulation technique to the distribution function for the categorical variable, estimated by means of indicator kriging f ⁎(u;k|(n). Calculating the mean value of the simulated data for the support of interest we obtain an estimate of the proportion of class ki in this support. 3.3. Multiple-point statistics Complex geometries, such as those reflecting slate deposits, cannot be modeled accurately using only traditional two-point statistics such as the variogram. It is possible to produce images with very different structures which all share the same two-point statistics (Strebelle, 2000), but two-point statistics are not enough to fully characterize low-entropy structures. This is why higher-order statistics have to be used to model those structures. Reproduction of such random geometries requires parameterization of specific shapes or consideration of joint categorical variability at 3 or more points at a time. This is why specific geometries are poorly reproduced by traditional pixel-based algorithms—such as indicator simulation technique—which manage only to reproduce proportions and two-point variograms. In multiple-point statistics (Strebelle, 2000), the value of a categorical variable Z(u)—or a continuous variable discretized in categories—belonging to a specific domain is determined from information on the spatial structure of the values for the variable, obtained from a training image of that domain. A data dn event of size n centred at location u is constituted by: • the data geometry defined by n vectors {hα, α = 1,… n} • the n data values {z(u + hα), α = 1,… n}. If only the geometric configuration is considered for a data event dn, then we refer to a data template τn. The aim of multiple-point statistics is to calculate the probability that a variable has a specific value (or attribute) at a point conditioned by a data event (geometric configuration and value) occurring in the region of this point. ProbðZ ðuÞ ¼ kjdn Þ:
ð1Þ
Noteworthy is the fact that whereas simulation techniques based on geostatistics use two-point statistics for the n data, multiple-point statistics uses n-point statistics in accordance with the template τn.
Given a template τn of n data variables, the inference of the probability distribution function conditional to a data event dn requires that there are enough occurrences of the combination (τn,dn) in the data sample. Since this is well nigh impossible in practice, multiple-point statistics makes use of a training image that depicts the geometries that are expected to be present in the real situation. Training images merely reflect a prior structural concept and need not carry any locally accurate information on the studied phenomenon. Object-based algorithms freed of the constraint of data conditioning can be used to generate such images. Another possible source of training images are photographs of outcrops processed with a CAD algorithm. In this research we chose to use object-based algorithms for the construction of the training image, primarily becuase of the great difficulty implied by the construction of a 3D image that would adequately represent the deposit structure from photographs of the mining fronts (that is, reconstruct a global 3D image from 2D discrete information). Under a prior decision of stationarity, the probability of the occurrence of a data event dn associated with τn can be inferred from the training image by calculating the training proportion ck(dn) / c(dn): ProbðZ ðuÞ ¼ kjdn Þ ¼
c k ð dn Þ c ð dn Þ
ð2Þ
where c(dn) is the number of replicates of the conditioning data event in the training image, and where ck(dn) is the number of replicates associated with a central value Z(u) equal to k. Reproducing large-scale structures would require the use of a large template, at least double the size of the structure. However, because of memory constraints, it is not practical to work with a large template (S-GeMS, 2007). Under the multiple grid concept, large-scale structures of the training image can be reproduced while considering a data template with a reasonably small number of grid nodes (Gómez-Hernández, 1991; Tran, 1994). The multiple grid approach consists of simulating a number G of increasingly finer grids, where the g-th (1 b g b G) grid is constituted of each 2g − 1 -th node of the final simulation grid (g = 1). The data template associated with the search neighbourhood of the fine grid is scaled proportionally to the spacing of the nodes within the grid to be simulated. When the g-th grid simulation is completed, its simulated values are frozen as data values to be used for conditioning on of the next finer simulation grid. Proper conditioning requires assigning the original sample data to the nearest nodes in the current simulation grid, but relocating the sample data to the coarser grids may affect the local accuracy of the simulated realizations. Multiple-point statistics can also account for a secondary variable, Y(u). The integration of several sources of information under the hypothesis of conditional independence (Journel, 2002) allows Prob(Z(u)=k|dn,Y(u)) to be calculated without the joint distribution of the data event dn and Y(u). This assumption basically states that the incremental information provided by one event of Y(u) before and after knowing the other events dn of Z(u) is constant.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
53
4. Studied area and base data The study was performed in a quarry located at the southern end of the Western Asturias-Leon region, in the Hercynian chain of the Iberian Peninsula, in Cabrera (province of Leon, Spain). The deposit, located in the Agüeira Formation (MiddleUpper Ordovician), is composed of a monotonous series of finegrained slate with some sandy laminations. The principal cleavage is N110E in direction with a dip of 60 S. Overlaying the exploitable level is a dense layer of kink bands running N100E in direction with a dip of 30 N. Kink bands are microfolds arising from late deformations in the Hercynian folding. They condition the exploitability of the mass if they are frequent and measurable in decimetre terms (Taboada et al., 2006a). The mining pit covers a surface area of approximately 475 m × 275 m, with a maximum difference in depths of around 100 m. Mining begins with the stripping of surface spoil, which may consist of other lithologies or areas of altered slate with no economical value, which is excavated directly or, when necessary, by means of perforation and blasting. The cutting system applied in the slate bank is to make deep cuts using diamond wire—a horizontal cut at the foot of the bank, followed by a vertical cut perpendicular to the cleavage direction. A bank is thus obtained in the form of a parallelepiped block measuring 3.5 m in the cleavage plane direction, 15 m in the perpendicular direction, and 5 m in height. The percentage of useful slate in the bank defines the grade. Blocks are extracted from the mass in accordance with the cleavage direction using the bucket teeth of a shovel loader or backhoe. The blocks are loaded onto trucks for transportation to the processing plant, where a series of mechanical processes are applied to produce commercial-quality slate. The average block/bank yield of useful slate for the quarry, defined as the quotient between the useful slate extracted in the form of blocks and the useful slate contained in the bank, is 52.6%. Plant yield, calculated on the basis of blocks inputs to the plant and slate production, is 15.82%. Therefore, the average saleable slate/bank yield for useful slate is 8.3%. These data were obtained from statistical analyzes of material flows in a year's mining operations. In order to compare deposit reserves using the different techniques, we used data extracted from 26 continuous drill cores. The deposit was discretized in blocks oriented according to the directions of the drill holes perforated and extraction of the blocks was planned along the time. For this purpose we used the block
Table 1 Anisotropy directions for the slate in the studied deposit Azimuth/direction (degrees)
Dip (degrees)
Description Direction of the cleavage– bedding intersection Maximum slope line for the kink bands Perpendicular to the kink bands
X'
95
−9
Y'
10
30
Z'
171
50
Fig. 1. Directional variograms obtained from the drill core data.
model estimated using indicator kriging and mine planning optimization algorithms (Bastante et al., 2004). The marginal cutoff grade of the orebody, lm, calculated using the maximum plant capacity, the cost and the average selling price, is 19% (Bastante et al., 2004). This grade distinguishes between the blocks that will be sawn with diamond wire and the blocks that—given their low grade—will be extracted directly by an excavator or by small-scale blasting (spoil). 5. Implementation of the techniques 5.1. Indicator kriging (ik) A detailed explanation of the estimation of deposit slate reserves using indicator kriging is given in Bastante et al. (2005). In order to use this technique in the evaluation of slate resources the information collected from the samples is first of all transformed into the binary system, with the value 1 (Class 1) assigned to useful slate and the value 0 (Class 0) assigned to non-useful slate. We can now implement indicator kriging (the direct kriging of coded data) on blocks in such a way that the result—the mathematical expectation or expected value of the indicator for each block—is an estimate of the proportion of material that corresponds to Class 1 (exploitable slate) within the chosen support. From the information extracted from the samples we constructed the experimental semivariograms in the 3 main anisotropy directions of the deposit, that is, X', Y', Z' (Table 1). The X' axis is oriented according to the lineation produced by the cleavagebedding intersection. The Y' axis is parallel to the planes formed by the kink bands. The semivariograms in both directions is quite similar. The semivariogram in the direction parallel to the Z'-axis and perpendicular to the kink bands has the greatest variability. The theoretical semivariogram, adjusted to the respective experimental semivariograms, is depicted in Fig. 1 (the original
54
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
mining pit boundaries (fewer reserves). By working with a fictitious cutoff grade (14%) we alter the tonnage of useful slate above the real marginal grade. Note, however, that in our case this variation will be very small, as the marginal grade is low and most of the reserves are in the richest blocks. 5.2. Sequential indicator simulation (sisim) The simulation was implemented as follows: Fig. 2. Training image used for the multiple-point analysis (transversal and longitudinal representation).
system of coordinates X, Y, Z was rotated until located parallel to the main anisotropic directions X', Y', Z'). The models are spherical nested models in all 3 cases. As typically happens with geostatistics, the theoretical semivariograms are constructed from the empirical semivariograms, on the basis of experience and of knowledge of the studied phenomenon and with the help of regression techniques. The goodness-of-fit of the results is analyzed using crossed validation. There is, therefore, a certain subjectivity in the construction that is unavoidable. Kriging was implemented once the directional semivariograms have been defined. Considered as the support for the calculation was a block measuring 14 m × 15 m × 5 m. Using this model we planned an annual extraction throughout the lifetime of the mine using optimization techniques (Bastante et al., 2004). In planning we have to bear in mind that the grade distribution function for the kriging model is only a crude estimate of the real distribution function, and that introducing the marginal grade will result in a percentage of in-quarry material for wire-cutting of 90%, when in reality we know that the real value of the operation is around 60%. Therefore we cannot use the real cost in the financial model, but will have to reduce it in such a way that the product of the new cost, multiplied by the percentage of the above-marginal grade material obtained using this cost (the fictitious cutoff grade) will be equal to the product of the real cost by the real percentage of 19% + grade material. In this way we avoid introducing a generalized bias in the financial evaluation that would affect the definition of the
1. 50 equiprobable realizations were generated on a point support (2 m × 2.5 m × 1.25 m grid) representing the volume of the deposit. The 2 possible values were 0 (non-useful slate) and 1 (useful slate). 2. For each realization the mean of the point values for the support used (14 m × 15 m × 5 m) was calculated for the blocks included in mine planning. 5.3. Multiple-point statistics (snesim algorithm) This technique requires a training image (Fig. 2). In this case, a training image was generated using an object-based algorithm called ellipsim (Deutsch and Journel, 1998), a program which simulates ellipsoids of different sizes and anisotropies until the proportion of points in their interior (in this case corresponding to useful slate) reaches a specific objective value. The training image was created by combining 5 classes of ellipsoids. The anisotropy directions of the ellipsoid axes were taken as equal to those for the variograms. The eccentricity relationships between the axes of the ellipsoids were obtained from the relationship between the ranges of the semivariograms and, to a certain degree, from the transversal profiles of the deposit interpreted by the mining expert. The size and proportion of each class were determined from the histogram of lengths of continuous stretches of useful slate created from the drill core data. Finally, the variogram of the training image was calculated, and comparison with the theoretical variogram indicated good results (Fig. 3). This coherence between the model and reality is no accident but arises as a consequence of the fact that, in order to generate the pattern image, the quantitative information supplied by the variograph and
Fig. 3. Theoretical and experimental variogram of the training image.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
the continuous cores was combined with the qualitative information supplied by the expert (the interpretation of the deposit reflected in the cross-sectional profiles). In order to infer the conditioned probability distribution from the training image for an event dn and to realize the simulations, the snesim algorithm was used. This algorithm implements a sequential simulation algorithm, called the single normal equation simulation (Strebelle, 2000), which simulates categorical variables using mps. The algorithm, implemented in the S-GeMS program (2007), functions on the basis of storing all training cpdfs in advance in a dynamic structure called a search tree in order to avoid scanning the training image anew for each unsampled node. See Strebelle and Journel (2001) for a detailed explanation of the snesim algorithm and its comparison with other algorithms. Snesim requires the introduction of a series of parameters that will determine the quality of the simulations. In our particular case, the most relevant parameters are the number of grid levels and the search neighbourhood (which should be adapted to the dimensions of the structures to be reproduced) and the maximum number of conditioning data items (which improve large-scale structure reproduction, although entailing a greater RAM and CPU cost). In Liu (2006), input parameters to the snesim algorithm are analyzed and recommendations are provided on how to set these parameters appropriately. The template τn used in the simulation consists of 25 points and takes the anisotropy directions of the variograms. In order to capture the structure over greater distances, 3 grids were created in accordance with the multiple grid concept described above. The simulation procedure was the same as for the indicator simulation: 50 realizations on a point support and integration in blocks included in mine planning. 5.4. Multiple-point statistics combined with indicator kriging (snesim + ik) One problem with simulation and multiple grids is the loss of local precision on assigning original data to the nearest nodes in the coarser grid. In order to overcome this problem the results for indicator kriging were introduced in the nodes of the coarser grid as a second source of information and the simulation was performed under a hypothesis of conditional independence between the information sources. Ortiz and Deutsch (2004) propose this combination of the information provided by indicator kriging and by mps for the simulation of a porphyry copper mine. In our case we consider that the hypothesis of conditional independence is a powerful one and we only want to avoid a loss of local precision in projecting the data onto the coarser grid. The same simulation procedure was used: 50 realizations on a point support and integration in blocks included in mine planning. 6. Results analysis Fig. 4 depicts a transversal and longitudinal representation of the deposit simulated using the different techniques. The longitudinal profile corresponds approximately to the east–west axis, with the direction of the lineation produced by the cleavage–bedding intersection shown in the form of a slight pitch of the slate struc-
55
Fig. 4. Sections of the continuous model of the deposit showing useful slate (dark grey) and non-useful slate (light grey). Above: conditional simulation model. Centre: multiple-point statistic model. Below: multiple-point statistic model reinforced with ordinary kriging.
tures. The transversal profile, running north–south, clearly defines how the useful slate structures are delimited according to the direction marked by the kink bands. It can be observed that the structures that shape the useful slate follow the anisotropy direction of the variogram, which is coherent with the geological interpretation of the deposit. It is evident that the models simulated with mps compared to sisim show areas of useful slate with more regular and continuous contours (greater connectivity and less entropy, after Journel and Deutsch, 1993). This is a consequence of using the training image and an n-dimensional template to simulate values for the unsampled points; unlike the conditional simulation, this leads to realizations with short-scale variations/noise since the simulation is based on information obtained solely from the pairs of points. These noisy features often seem to be geologically unrealistic and make the realizations less credible to geologists (Deutsch, 1998). Table 2 shows the overall estimate for grade mean and standard deviation, for all the blocks in the deposit to be extracted according to the mining plan. The standard deviation for the mean for the ik model was calculated from the overall variance for the kriging. For the simulated models the results corresponded to the values obtained for the 50 simulations. Established as the reference for the mean grade was the ik model, given that this was the best unbiased linear estimator (BLUE). Note that for all the cases the standard deviation for the mean was around 4%. The snesim model, however, gives the value that is farthest from the mean with respect to ik, which we
56
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
Table 2 Statistics for the 4 techniques used to evaluate useful slate
Ik Sisim Snesim Snesim + ik
Grade mean (%)
Standard deviation (%)
49.3 51.1 46.1 48.1
4.1 3.7 3.7 4.7
The standard deviation for ik was calculated from the overall variance for the estimation.
believe to be due to the loss of local precision arising from the projection of the data in the coarsest grid (as has been commented earlier). For this reason the snesim + ik model was included in the overall estimate of the mean grade (48.1%), at the cost, however, of increasing the standard deviation (4.7%). Of interest is Fig. 5, which represents the grade distribution function for the blocks obtained using the different methods. For the simulated models the results correspond to the composition of the 50 histograms from the percentiles. Indicator kriging shows significant clustering around the mean (a small degree of dispersion), with a smoothing of the results. This is because the estimation aims at local accuracy so as to reduce the dispersion of the results in a magnitude that is approximately equal to the mean kriging variance of the support used. Basically, small values are overestimated, whereas larger values are underestimated. Sisim, which has a distribution function that is practically uniform, increases the dispersion of the grades. Although the sequential indicator simulation aims at reproducing the two-point variability of the real underlying phenomena, when large-range continuity is present, the traditional variogram-based simulation approaches do not provide good n-point reproduction. The algorithm tends to produce realizations with high entropy/low connectivity, and, in our case, this is further aggravated by the practically uniform distribution of the starting data. Combining the results for the points so as to calculate the grade for each block obtains a uniform grade distribution, as the blocks are large relative to the—noisy—simulated structures (that is, the information for the blocks contains a lot of noise). Finally, snesim and snesim + ik show more clustering at the extremes of the distribution function (a large degree of dis-
Fig. 5. Distribution functions of useful slate frequencies obtained via ik, sisim and mps.
persion). Note that the greater results variability for the snesim + ik application arises from the hypothesis of conditional independence. Unlike the sisim, snesim takes into account the moments of an order greater than 2, thereby reflecting the continuity pattern for the training image. The realizations clearly differentiate between large zones of useful slate and large areas of non-useful slate. Grouping the information in the blocks, we observe that a high proportion have very high or very low grades, with an increased dispersion in the grade distribution. With snesim + ik this effect is multiplied on originating useful slate structures of larger dimensions although fewer in number. These two models better reflect the reality of the mine, in which the miner typically encounters zones of spoil (poor in useful slate) or of mineral (rich in useful slate); in other words, it is not normal to encounter areas with grades of around 50%. Furthermore, bearing in mind that the distribution function for the point information (the drill cores) is represented by 50% of the data with values equal to 0 (waste slate) and the remaining 50% with values equal to 1 (useful slate), intuitively it can be deduced that the effect of grouping the information in a larger support would be precisely that represented by the distribution functions obtained with the mps algorithms. Fig. 5 shows how introducing the cutoff grade into each of the distribution functions affects the finances of the operation. Since the cutoff grade is low (19%), the total percentage of useful slate recovered from the blocks after diamond-wire cutting will not vary in relative terms by more than 3% from one model to another. What will be very different, however, is the total cost of perforation and diamond-wire cutting, given that the percentage of blocks that surpass the grade varies greatly between the models generated by the geostatistical techniques and using mps. This is where the reason for using a fictitious cost and cutoff grade in the mine plan for the krigged block model becomes evident—it is a way of reducing financial bias in the mine evaluation. In fact, as mentioned previously, approximately 60% of the orebody is cut, a value which is very close to that obtained with the mps technique. The difference between the models, in fact, is even smaller than can be deduced from the graph, given that a higher cutoff grade than the marginal grade of around 25% is used in the quarry. Fig. 6 depicts the omnidirectional variogram for the different estimated and simulated block models (for the latter a model
Fig. 6. Omnidirectional variograms for the kriging and simulated models.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
considered as representative was chosen). What is clear from the graph is that the fact that a simulation technique may be capable of reproducing the variogram of the points does not mean that it can reproduce the variogram on a different support when the simulated data is grouped and averaged. Nor should it be overlooked that although the simulation techniques theoretically permit the calculation of an infinite number of equiprobable realizations of the statistical model that reflect reality, the model/reality approximation will only be effective in the terms previously dictated by the model, whether approximation to the mean, the histogram, variogram function, or n-point statistics. Paradoxically, of the different techniques used in the variogram obtained, snesim was the closest to the geostatistical model, with a sill that practically coincided with the theoretical variance in dispersion calculated from the theoretical variogram. The effect of introducing the ik model information in some of the grid points of the snesim algorithm is again evident, with a sill that is substantially greater than the theoretical sill; also reflected is the lesser dispersion of sisim and, logically, of ik. These results are coherent with the arguments described above in relation to the dispersion of the distribution functions for the block grades, and, to a certain extent, confirm that the training image used is appropriate. It should be recalled that the effect of introducing the conditional independence hypothesis is to increase the size/continuity of the structures. This explains the greater range and the greater plateau for the variograms obtained from the realizations of snesim + ik. Table 3 shows deviations in the total production of saleable slate over the life of the mine—taking as reference the model ik— and also the coefficient of variation. In this case (although this cannot be established conclusively), the combination snesim + ik seems to adapt better to reality if we take into account that the ik model overestimates somewhat the total production by using an estimated model and a fictitious cutoff grade. Due to the bias produced by the projection of the data, snesim results in excessively low production. If we graph the mean annual deviation in production for each model with respect to the ik model, it will be seen that only snesim + ik locates the values at around zero (below zero in fact). Given that each year some 150 blocks of slate are extracted, this would seem to indicate that snesim + ik does not introduce local bias (at least for this time interval or number of blocks). An interesting exercise was the calculation of the distribution function for the discounted cash flows with a view to quantifying how uncertainty as to knowledge of the deposit affects uncertainty as to the profitability of the mining operation. This requires the financial model and mining sequence to be introduced into each of Table 3 Deviation in total production and coefficient of variation (CV) after applying the cutoff grade (19%)
Sisim Snesim Snesim + ik Reference model: ik.
Deviation (%)
CV
+2.7 − 7.6 − 2.9
7.6 8.3 10
57
Fig. 7. Discounted cash flows for the simulation techniques and for ordinary kriging.
the simulated deposits prior to implementing the financial analysis, this time with the original marginal cutoff grade, that is, without modifying the unit cost of perforation and diamond-wire cutting, as was done with the kriging. Fig. 7 shows the results: the X axis represents the discounted cash flows and the ordinate axis shows the percentiles deriving from the 50 simulated block models for each technique. In all cases the median is below the value obtained by kriging. This result can be interpreted as contradictory in the sisim case, as this produces a positive bias in total slate production with respect to the value obtained by kriging. However, this is not the case if we observe the grade distribution function for the blocks obtained with sisim (Fig. 5). The total cost of perforation and diamond-wire cutting assigned by sisim is practically 37% greater than the real cost (82% of the blocks will be cut in accordance with the graph, compared to the real quarry value of 60%). In contrast, the values obtained using snesim and snesim + ik are very close to reality (64% and 62% of cut blocks, respectively), and so these models will introduce no significant degree of bias in the cost. Therefore, the strong negative deviation in the value of the median presented by snesim with respect to ik is due to the bias introduced into total production (Table 3). Despite this bias that emerges with the sisim and snesim simulation techniques—displacement of the median/mean with respect to the value of ik—the variance of the distribution function for the discounted cash flows with respect to the mean in both cases is around 130%2 or, put in another way, the coefficient of variation is around 11.4%. Moreover, snesim + ik has an expected value that is quite close to that obtained with ik but with a variance value that is somewhat greater. This is due to the double effect of introducing—under the conditional independence hypothesis—certain values of ik into snesim, improving local accuracy at the cost of increasing the dispersion of the results. 7. Conclusions Traditional variogram-based geostatistics and multiple-point statistics rely on the same assumption of stationarity. Otherwise estimation or simulation would not be possible. From this
58
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
perspective, all these techniques have the same limitations, although there is a crucial difference between them: whereas in the mps there is control over the information concerning the higher-order statistics in the geostatistical simulation that's not possible (with the exception of the case in which the sequential Gaussian simulation algorithm is used under the hypothesis that the multiple-point statistics are multivariate Gaussian). In an indirect way, the algorithm is given the power to decide over these statistics. These techniques also have the same requirement, namely, the need to construct a variogram/training image that is representative of the phenomenon being studied. The quality of the results will depend fundamentally on these tools, which are the bases for calculating the simulations. In this article we analyzed the use of different techniques for modelling the reserves in a slate deposit. The results that best reflected the reality were those obtained using multiple-point statistics, which best reproduced the structure of the deposit, both from a statistical and geological perspective. The advantage of the mps algorithms compared to the sisim algorithms is that they produce simulations with a greater connectivity between points— less entropy—and that better represent the variability in the deposit, which is even more manifest when working with nonpoint supports. The main drawback is the need to create a training image, which is also a complex matter. Comparatively it is more complex to construct a training image—which has to reflect the n-point statistics for the phenomenon—than a two-point variogram. Moreover, the variogram construction task can use adjustment techniques through regression, and this is undoubtedly very useful when working in three dimensions. The crossed validation, furthermore, enables the fit of the results to be compared when different theoretical variograms are used. The mps models provide useful slate percentages—above the cutoff grade—that are much close to reality than those based on geostatistics. This is fundamental when calculating slate reserves, costs and revenue flows, and ultimately, when evaluating the mine from a financial perspective. An interesting aspect of the evaluation of the deposit is that—assuming that the remaining parameters remain constant—the financial risk profile derived exclusively from uncertainty in regard to knowledge of the deposit has, for all the techniques, a coefficient of variation that is relatively small if we consider the high variability in the variogram function with respect to the distance and number of cores used in evaluating the studied area. Returning to the evaluation of the deposit using mps, we emphasize that the loss of local precision on using multiple grids with the snesim algorithm can be avoided by introducing part of the information supplied by ik into the coarser grid used in the simulation, although at the cost of increasing the dispersion variance. Finally, it should be underlined that when modelling and evaluating a deposit, the ultimate aim of the mine planner is to design and plan mining in such a way as to maximize net present value (NPV)—in statistical terms, expected value. The importance of assessment and evaluation of the risk associated with mining (as mentioned earlier in this article) should not be taken for granted. Nonetheless, it should be pointed out that the current
approach to mining design requires a stochastic focus of the deposit; consequently, design based on optimization algorithms should start with simulated rather than estimated models of the deposit (Dimitrakopoulos et al., 2007; Bastante et al., in press). Given the availability of the software, computers, algorithms and models necessary to implement this focus, all we need are experts who combine multidisciplinary training with consolidated on-site experience. Acknowledgment This research has been partially supported by Grant BIA200766218 of the Spanish Ministry of Education and Science. References AENOR: Asociación Española de Normalización y Certificación, 2005. UNE 12326-1: 2005. Productos de pizarra y piedra natural para tejados inclinados y revestimientos discontinuos. Parte 1: Especificación de producto. Madrid. Bastante, F.G., Taboada, J., Ordóñez, C., 2004. Design and planning for slate mining using optimisation algorithms. Engineering Geology 73, 93–103. Bastante, F.G., Taboada, J., Alejano, L.R., Ordóñez, C., 2005. Evaluation of the resources of a slate deposit using indicator kriging. Engineering Geology 81, 407–418. Bastante, F. G., Taboada, J., Alejano, L., Alonso, E. in press. Optimization tools and simulation methods to design and evaluate a mining operation. Stochastic Environmental Research & Risk Assessment. DOI 10.1007/ s00477-007-0182-6. Deutsch, C.V., 1998. Cleaning categorical variable (lithofacies) realizations with maximum a-posteriori selection. Computers & Geosciences 24 (6), 551–562. Deutsch, C.V., Journel, A.G., 1998. GSLIB. Geostatistical Software Library and User's Guide, Applied Geostatistical Software, second edition. Oxford University Press, Inc., New York. Dimitrakopoulos, R., Farrelly, C., Godoy, M., 2002. Moving forward from traditional optimisation: grade uncertainty and risk effects in open pit mine design. Transactions of the IMM, Section A Mining Industry 111, A82–A89. Dimitrakopoulos, R., Martinez, L., Ramazan, S., 2007. A maximum upside / minimum downside approach to the traditional optimization of open pit mine design. Journal of Mining Science 43 (1), 73–82. Gómez-Hernández, J., 1991. A stochastic approach to the simulation of block fields conductivity upon data measured at a smaller scale. Ph.D. thesis, Stanford University, Stanford. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford Univ. Press, New York, USA. Guardiano, F., Srivastava, R.M., 1993. In: Soares, A. (Ed.), Multivariate Geostatistics: Beyond Bivariate Moments, Geostatistics-Troia, vol. 1. Kluwer Academic Publications, Dordrecht, pp. 133–144. Journel, A., 1974. Geostatistics for conditional simulation of ore-bodies. Economic Geology 69, 673–687. Journel, A.G., Huijbregts, C.J., 1978. Mining Geostatistics. Academic Press, London. Journel, A.G., Deutsch, C.V., 1993. Entropy and spatial disorder. Mathematical Geology 25 (3), 329–355. Journel, A., 2002. Combining knowledge from diverse sources: an alternative to traditional data independence hypotheses. Mathematical Geology 34 (5), 573–596. Liu, Y., 2006. Using the Snesim program for multiple-point statistical simulation. Computers & Geosciences 32, 1544–1563. Ortiz, J.M., Deutsch, C.V., 2004. Indicator simulation accounting for multiplepoint statistics. Mathematical Geology 36 (5), 545–565. S-GeMS, The Standford Geostatistical Modeling Software: Product information on the Internet. Available online at: http://sgems.sourceforge.net (accessed July 2007). Strebelle, S.B., 2000. Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology, 34, 1–22.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59 Strebelle, S.B., Journel, A.G., 2001. Reservoir modeling using multiple-point statistics. SPE International 71324, 1–11. Taboada, J., Ordóñez, C., Saavedra, A., Fiestras-Janeiro, G., 2006a. Fuzzy expert system for economic zonation of an ornamental slate deposit. Engineering Geology 84, 220–228.
59
Taboada, J., Matías, J.M., Araújo, M., Ordóñez, C., 2006b. Assessing the viability of underground slate mining by combining an expert system with a GIS. Engineering Geology 87, 75–84. Tran, T., 1994. Improving variogram reproduction on dense simulation grids. Computers & Geosciences 20, 1161–1168.
Engineering Geology 98 (2008) 50 – 59 www.elsevier.com/locate/enggeo
Comparison of indicator kriging, conditional indicator simulation and multiple-point statistics used to model slate deposits F.G. Bastante a , C. Ordóñez a,⁎, J. Taboada a , J.M. Matías b a
Department of Natural Resources and Environmental Engineering, University of Vigo, Vigo, Spain b Department of Statistics, University of Vigo, Vigo, Spain
Received 17 September 2007; received in revised form 25 December 2007; accepted 26 January 2008 Available online 8 February 2008
Abstract The resources in an ornamental slate deposit can be estimated using geostatistical estimation techniques applied to information collected from drill cores. The result, however, is a smooth approximation that fails to take account of the natural variability in mineralization, which is fundamental to proper design and evaluation of the financial viability of a mining deposit. Geostatistical simulation techniques are more useful in this respect, as they reflect different realizations of the reality and better reflect natural mineral dispersion. In this work, we evaluate the resources of a slate quarry by comparing the results obtained using two geostatistical techniques—indicator kriging (ik) and sequential indicator simulation (sisim)—with the results obtained using the single normal equation simulation (snesim) technique based on multiple-point statistics (mps), analyzing their usefulness in evaluating the financial risk derive from uncertainty in regard to knowledge of the deposit. Our results indicate that although the multiple-point statistics approach produces models that are closer to reality than the models produced by the geostatistical techniques, the simulation relies in part on information obtained via indicator kriging. © 2008 Elsevier B.V. All rights reserved. Keywords: Geostatistics; Indicator kriging; Conditional simulation; Multiple-point statistics; Mining resource evaluation
1. Introduction The mining of ornamental rocks, and especially slate for roofing, has undergone spectacular development during the last 30 years in Spain, a fact which has placed this country first among the world's slate producers. This has brought about a considerable improvement in both productive processes and the machinery used, but there are a number of technical factors that have not yet been satisfactorily resolved. One example is the exploration and evaluation of slate deposits for mining purposes. Slate is a fine-grained, homogeneous, metamorphic rock derived from an original shale-type sedimentary rock composed of clay through low-grade regional metamorphism. It is mainly
⁎ Corresponding author. E.T.S.I.MINAS, Universidad de Vigo, Campus Universitario, 36310 Vigo, Spain. Tel.: +34 986814052; fax: +34 986812201. E-mail address: [email protected] (C. Ordóñez). 0013-7952/$ - see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.enggeo.2008.01.006
composed of quartz and muscovite or illite, often along with biotite, chlorite, hematite and pyrite, and less frequently, apatite, graphite, kaolin, magnetite, tourmaline, zircon and feldspar. As methods for slate quarrying, we can distinguish clearly between the more frequently utilized surface mining method, and the underground method that is gradually becoming more popular (Taboada et al., 2006a). Slate is commonly used as a roofing and flooring material and for electrical insulation. Geological investigation of a slate deposit is mainly based on continuous drill cores, since surface alterations impede surface sampling of fresh rock. The cores are used to study the sedimentary, metamorphic and structural characteristics that define the orebody. The orebody typically contains a mass of useful slate (which is defined, as we shall see subsequently, as the equivalent of mineral grade) and areas of spoil, the latter mainly a consequence of structural discontinuities in the orebody. In slate mining as in the rest of mining, however, all knowledge is inevitably linked to uncertainty, and this is represented as risk when a mining deposit is being assessed. In order to obtain the necessary decision-making
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
elements for each of the developmental phases of the mining project, the mining company has to quantitatively calculate a value for each parameter that intervenes in the mining evaluation along with the uncertainty associated with these values. This uncertainty has a significant bearing on the mining company's final decision as to whether to assume the risk of allocating human and financial resources to the development of a mining project. The problem of grade uncertainty and risk effects in open-pit design is addressed by Dimitrakopoulos et al. (2002), who use generalized sequential Gaussian simulation and direct block simulation techniques to model and evaluate the financial risk implied by a typical, disseminated, low-grade, epithermal, quartz breccia-type gold deposit. In this article, we tackle the same problem for slate deposits, although taking a somewhat different approach. We estimate useful slate reserves using indicator kriging and compare the results with those obtained using conditional simulation techniques and multiple-point statistics (Guardiano and Srivastava, 1993). 2. Roofing slate properties: Useful slate and saleable slate Only part of the slate in a quarry can be exploited as roofing slate, as slate for this purpose must have certain physical, chemical and aesthetic properties and must not be affected by structures that interrupt the cleavage direction. This slate is referred to as useful slate. From a general perspective, the usefulness of slate for roofing is determined from samples extracted from quarry fronts which are laboratory tested for compliance with the Spanish UNE standard 12326-1: 2005 (AENOR, 2005). These tests analyze slate samples in terms of a range of physical and chemical properties such as fissility, impermeability, carbonate content, colour, alterability, inclusions, homogeneity of grain size, and roughness. (See Taboada et al., 2006b) for a more detailed description of these properties). The test results combined with data from continuous drill cores enable experts to determine which parts of the drill cores represent useful slate. Thus, the expert divides the cores into 1-m long stretches and assigns them to Class 1 if classified as useful slate and Class 0 if classified as non-useful slate. The choice of 1-m lengths is justified by the fact that mining requires a minimum block thickness of 50 cm; consequently, as thickness is given by the normal direction to the cleavage plane that dips 60°, the support can be considered to be representative of the minimum block size in that direction. As will be demonstrated in Section 5.1 below, using indicator kriging will enable us to estimate both the percentage of stretches belonging to Class 1 within the blocks into which the quarry will be discretized and the percentage of useful slate contained in this support (in other words, the mineral grade for the block). In this way gross slate reserves can be estimated and mining can be designed optimally (Bastante et al., 2005). Of the useful slate indicated as such by the drill core sampling, only part can be recovered from the quarry in the form of blocks. Only a small proportion of this recovered slate is saleable, moreover, and only after the slate has undergone various mechanical operations in the slate transformation processing plant. This slate is referred to as finished or saleable slate.
51
3. Methodologies used to evaluate slate resources In our research, we used different techniques to evaluate useful slate reserves from the information provided by cores made in a slate quarry. Below we summarize the theory underlying each of the techniques and subsequently compare and discuss the results obtained for each method. 3.1. Indicator kriging The basic approach taken by predictive statistics is to model the uncertainty associated with an unsampled value, z(u), with u as the coordinate location vector, as a random variable (RV), Z(u), the probability distribution of which characterizes the uncertainty about z(u) (Deutsch and Journel, 1998). Categorical variables such as rock types can be effectively modelled by RVs. Let Z(u) be a categorical RV that can take either of 2 outcome values k = {k1, k2}. In the particular case to be analyzed here, z(u) is the variable stretch of a metre of slate centered on u; Z(u) is the random function that models z(u); and (k1, k2) corresponds to the categories (Class 1, Class 0) or useful/non-useful slate. The conditional probability density function is: f ðu; kjðnÞÞ ¼ProbfZ ðuÞ a category kjðnÞg (n) consisting of n neighbouring data values z(uα), α = 1,…, n. The essence of the indicator approach is the binomial coding of Z(u) into either 1 or 0 depending upon its relationship to a categorical class. By considering binary indicator transforms of Z(u) defined as: 1; if Z ðuÞ a k I ðu; k Þ ¼ 0; otherwise the probability that a category k prevails at u can be expressed as: EfI ðu; kjðnÞÞg ¼ Probfðu; k Þ ¼ 1jðnÞg ¼ f ðu; kjðnÞÞ: Applying the kriging algorithm to indicator data we obtain the least-squares estimate of its conditional expectation, f ⁎(u;k| (n). Note that the bivariate (two-point) distribution of 2 RV's Z (u1), Z(u2), appears as the non-centered covariance of the indicator variables: f ðu1 ; u2 ; k1 ; k2 Þ ¼ProbfZ ðu1 Þ ¼ k1 ; Z ðu2 Þ ¼ k2 g ¼ E fI ðu1 ; k1 ÞI ðu2 ; k2 Þg and the indicator variogram function measures how often two locations that are a vector h apart belong to different categories. Indicators are useful for characterizing the spatial variability of categorical variables as the ranges and shapes of the directional indicator semivariograms that reflect geometric patterns (Goovaerts, 1997). 3.2. Sequential indicator simulation Estimating reserves by means of kriging produces a smoothed representation of the reality, and this typically underestimates real data dispersion (Journel, 1974). Since geostatistics considers the
52
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
variable that is the object of the study to be a realization of a random function, it is possible to simulate an infinite number of realizations of random functions to represent the variability and spatial correlation characteristics observed in the experimental data (Journel and Huijbregts, 1978). The conditional simulation techniques reproduce the first two experimental moments (mean and covariance, or variogram) for the real data; in other words, they reproduce a measure of real dispersion in accordance with the information available, with the variogram acting as a two-point measure of geological heterogeneity or continuity. Therefore, whereas kriging provides the best unbiased estimator for each point, the conditional simulation produces values that reproduce real fluctuations. In the sequential indicator simulation technique, the simulated value for each point (1 or 0) is obtained by applying the Monte Carlo simulation technique to the distribution function for the categorical variable, estimated by means of indicator kriging f ⁎(u;k|(n). Calculating the mean value of the simulated data for the support of interest we obtain an estimate of the proportion of class ki in this support. 3.3. Multiple-point statistics Complex geometries, such as those reflecting slate deposits, cannot be modeled accurately using only traditional two-point statistics such as the variogram. It is possible to produce images with very different structures which all share the same two-point statistics (Strebelle, 2000), but two-point statistics are not enough to fully characterize low-entropy structures. This is why higher-order statistics have to be used to model those structures. Reproduction of such random geometries requires parameterization of specific shapes or consideration of joint categorical variability at 3 or more points at a time. This is why specific geometries are poorly reproduced by traditional pixel-based algorithms—such as indicator simulation technique—which manage only to reproduce proportions and two-point variograms. In multiple-point statistics (Strebelle, 2000), the value of a categorical variable Z(u)—or a continuous variable discretized in categories—belonging to a specific domain is determined from information on the spatial structure of the values for the variable, obtained from a training image of that domain. A data dn event of size n centred at location u is constituted by: • the data geometry defined by n vectors {hα, α = 1,… n} • the n data values {z(u + hα), α = 1,… n}. If only the geometric configuration is considered for a data event dn, then we refer to a data template τn. The aim of multiple-point statistics is to calculate the probability that a variable has a specific value (or attribute) at a point conditioned by a data event (geometric configuration and value) occurring in the region of this point. ProbðZ ðuÞ ¼ kjdn Þ:
ð1Þ
Noteworthy is the fact that whereas simulation techniques based on geostatistics use two-point statistics for the n data, multiple-point statistics uses n-point statistics in accordance with the template τn.
Given a template τn of n data variables, the inference of the probability distribution function conditional to a data event dn requires that there are enough occurrences of the combination (τn,dn) in the data sample. Since this is well nigh impossible in practice, multiple-point statistics makes use of a training image that depicts the geometries that are expected to be present in the real situation. Training images merely reflect a prior structural concept and need not carry any locally accurate information on the studied phenomenon. Object-based algorithms freed of the constraint of data conditioning can be used to generate such images. Another possible source of training images are photographs of outcrops processed with a CAD algorithm. In this research we chose to use object-based algorithms for the construction of the training image, primarily becuase of the great difficulty implied by the construction of a 3D image that would adequately represent the deposit structure from photographs of the mining fronts (that is, reconstruct a global 3D image from 2D discrete information). Under a prior decision of stationarity, the probability of the occurrence of a data event dn associated with τn can be inferred from the training image by calculating the training proportion ck(dn) / c(dn): ProbðZ ðuÞ ¼ kjdn Þ ¼
c k ð dn Þ c ð dn Þ
ð2Þ
where c(dn) is the number of replicates of the conditioning data event in the training image, and where ck(dn) is the number of replicates associated with a central value Z(u) equal to k. Reproducing large-scale structures would require the use of a large template, at least double the size of the structure. However, because of memory constraints, it is not practical to work with a large template (S-GeMS, 2007). Under the multiple grid concept, large-scale structures of the training image can be reproduced while considering a data template with a reasonably small number of grid nodes (Gómez-Hernández, 1991; Tran, 1994). The multiple grid approach consists of simulating a number G of increasingly finer grids, where the g-th (1 b g b G) grid is constituted of each 2g − 1 -th node of the final simulation grid (g = 1). The data template associated with the search neighbourhood of the fine grid is scaled proportionally to the spacing of the nodes within the grid to be simulated. When the g-th grid simulation is completed, its simulated values are frozen as data values to be used for conditioning on of the next finer simulation grid. Proper conditioning requires assigning the original sample data to the nearest nodes in the current simulation grid, but relocating the sample data to the coarser grids may affect the local accuracy of the simulated realizations. Multiple-point statistics can also account for a secondary variable, Y(u). The integration of several sources of information under the hypothesis of conditional independence (Journel, 2002) allows Prob(Z(u)=k|dn,Y(u)) to be calculated without the joint distribution of the data event dn and Y(u). This assumption basically states that the incremental information provided by one event of Y(u) before and after knowing the other events dn of Z(u) is constant.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
53
4. Studied area and base data The study was performed in a quarry located at the southern end of the Western Asturias-Leon region, in the Hercynian chain of the Iberian Peninsula, in Cabrera (province of Leon, Spain). The deposit, located in the Agüeira Formation (MiddleUpper Ordovician), is composed of a monotonous series of finegrained slate with some sandy laminations. The principal cleavage is N110E in direction with a dip of 60 S. Overlaying the exploitable level is a dense layer of kink bands running N100E in direction with a dip of 30 N. Kink bands are microfolds arising from late deformations in the Hercynian folding. They condition the exploitability of the mass if they are frequent and measurable in decimetre terms (Taboada et al., 2006a). The mining pit covers a surface area of approximately 475 m × 275 m, with a maximum difference in depths of around 100 m. Mining begins with the stripping of surface spoil, which may consist of other lithologies or areas of altered slate with no economical value, which is excavated directly or, when necessary, by means of perforation and blasting. The cutting system applied in the slate bank is to make deep cuts using diamond wire—a horizontal cut at the foot of the bank, followed by a vertical cut perpendicular to the cleavage direction. A bank is thus obtained in the form of a parallelepiped block measuring 3.5 m in the cleavage plane direction, 15 m in the perpendicular direction, and 5 m in height. The percentage of useful slate in the bank defines the grade. Blocks are extracted from the mass in accordance with the cleavage direction using the bucket teeth of a shovel loader or backhoe. The blocks are loaded onto trucks for transportation to the processing plant, where a series of mechanical processes are applied to produce commercial-quality slate. The average block/bank yield of useful slate for the quarry, defined as the quotient between the useful slate extracted in the form of blocks and the useful slate contained in the bank, is 52.6%. Plant yield, calculated on the basis of blocks inputs to the plant and slate production, is 15.82%. Therefore, the average saleable slate/bank yield for useful slate is 8.3%. These data were obtained from statistical analyzes of material flows in a year's mining operations. In order to compare deposit reserves using the different techniques, we used data extracted from 26 continuous drill cores. The deposit was discretized in blocks oriented according to the directions of the drill holes perforated and extraction of the blocks was planned along the time. For this purpose we used the block
Table 1 Anisotropy directions for the slate in the studied deposit Azimuth/direction (degrees)
Dip (degrees)
Description Direction of the cleavage– bedding intersection Maximum slope line for the kink bands Perpendicular to the kink bands
X'
95
−9
Y'
10
30
Z'
171
50
Fig. 1. Directional variograms obtained from the drill core data.
model estimated using indicator kriging and mine planning optimization algorithms (Bastante et al., 2004). The marginal cutoff grade of the orebody, lm, calculated using the maximum plant capacity, the cost and the average selling price, is 19% (Bastante et al., 2004). This grade distinguishes between the blocks that will be sawn with diamond wire and the blocks that—given their low grade—will be extracted directly by an excavator or by small-scale blasting (spoil). 5. Implementation of the techniques 5.1. Indicator kriging (ik) A detailed explanation of the estimation of deposit slate reserves using indicator kriging is given in Bastante et al. (2005). In order to use this technique in the evaluation of slate resources the information collected from the samples is first of all transformed into the binary system, with the value 1 (Class 1) assigned to useful slate and the value 0 (Class 0) assigned to non-useful slate. We can now implement indicator kriging (the direct kriging of coded data) on blocks in such a way that the result—the mathematical expectation or expected value of the indicator for each block—is an estimate of the proportion of material that corresponds to Class 1 (exploitable slate) within the chosen support. From the information extracted from the samples we constructed the experimental semivariograms in the 3 main anisotropy directions of the deposit, that is, X', Y', Z' (Table 1). The X' axis is oriented according to the lineation produced by the cleavagebedding intersection. The Y' axis is parallel to the planes formed by the kink bands. The semivariograms in both directions is quite similar. The semivariogram in the direction parallel to the Z'-axis and perpendicular to the kink bands has the greatest variability. The theoretical semivariogram, adjusted to the respective experimental semivariograms, is depicted in Fig. 1 (the original
54
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
mining pit boundaries (fewer reserves). By working with a fictitious cutoff grade (14%) we alter the tonnage of useful slate above the real marginal grade. Note, however, that in our case this variation will be very small, as the marginal grade is low and most of the reserves are in the richest blocks. 5.2. Sequential indicator simulation (sisim) The simulation was implemented as follows: Fig. 2. Training image used for the multiple-point analysis (transversal and longitudinal representation).
system of coordinates X, Y, Z was rotated until located parallel to the main anisotropic directions X', Y', Z'). The models are spherical nested models in all 3 cases. As typically happens with geostatistics, the theoretical semivariograms are constructed from the empirical semivariograms, on the basis of experience and of knowledge of the studied phenomenon and with the help of regression techniques. The goodness-of-fit of the results is analyzed using crossed validation. There is, therefore, a certain subjectivity in the construction that is unavoidable. Kriging was implemented once the directional semivariograms have been defined. Considered as the support for the calculation was a block measuring 14 m × 15 m × 5 m. Using this model we planned an annual extraction throughout the lifetime of the mine using optimization techniques (Bastante et al., 2004). In planning we have to bear in mind that the grade distribution function for the kriging model is only a crude estimate of the real distribution function, and that introducing the marginal grade will result in a percentage of in-quarry material for wire-cutting of 90%, when in reality we know that the real value of the operation is around 60%. Therefore we cannot use the real cost in the financial model, but will have to reduce it in such a way that the product of the new cost, multiplied by the percentage of the above-marginal grade material obtained using this cost (the fictitious cutoff grade) will be equal to the product of the real cost by the real percentage of 19% + grade material. In this way we avoid introducing a generalized bias in the financial evaluation that would affect the definition of the
1. 50 equiprobable realizations were generated on a point support (2 m × 2.5 m × 1.25 m grid) representing the volume of the deposit. The 2 possible values were 0 (non-useful slate) and 1 (useful slate). 2. For each realization the mean of the point values for the support used (14 m × 15 m × 5 m) was calculated for the blocks included in mine planning. 5.3. Multiple-point statistics (snesim algorithm) This technique requires a training image (Fig. 2). In this case, a training image was generated using an object-based algorithm called ellipsim (Deutsch and Journel, 1998), a program which simulates ellipsoids of different sizes and anisotropies until the proportion of points in their interior (in this case corresponding to useful slate) reaches a specific objective value. The training image was created by combining 5 classes of ellipsoids. The anisotropy directions of the ellipsoid axes were taken as equal to those for the variograms. The eccentricity relationships between the axes of the ellipsoids were obtained from the relationship between the ranges of the semivariograms and, to a certain degree, from the transversal profiles of the deposit interpreted by the mining expert. The size and proportion of each class were determined from the histogram of lengths of continuous stretches of useful slate created from the drill core data. Finally, the variogram of the training image was calculated, and comparison with the theoretical variogram indicated good results (Fig. 3). This coherence between the model and reality is no accident but arises as a consequence of the fact that, in order to generate the pattern image, the quantitative information supplied by the variograph and
Fig. 3. Theoretical and experimental variogram of the training image.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
the continuous cores was combined with the qualitative information supplied by the expert (the interpretation of the deposit reflected in the cross-sectional profiles). In order to infer the conditioned probability distribution from the training image for an event dn and to realize the simulations, the snesim algorithm was used. This algorithm implements a sequential simulation algorithm, called the single normal equation simulation (Strebelle, 2000), which simulates categorical variables using mps. The algorithm, implemented in the S-GeMS program (2007), functions on the basis of storing all training cpdfs in advance in a dynamic structure called a search tree in order to avoid scanning the training image anew for each unsampled node. See Strebelle and Journel (2001) for a detailed explanation of the snesim algorithm and its comparison with other algorithms. Snesim requires the introduction of a series of parameters that will determine the quality of the simulations. In our particular case, the most relevant parameters are the number of grid levels and the search neighbourhood (which should be adapted to the dimensions of the structures to be reproduced) and the maximum number of conditioning data items (which improve large-scale structure reproduction, although entailing a greater RAM and CPU cost). In Liu (2006), input parameters to the snesim algorithm are analyzed and recommendations are provided on how to set these parameters appropriately. The template τn used in the simulation consists of 25 points and takes the anisotropy directions of the variograms. In order to capture the structure over greater distances, 3 grids were created in accordance with the multiple grid concept described above. The simulation procedure was the same as for the indicator simulation: 50 realizations on a point support and integration in blocks included in mine planning. 5.4. Multiple-point statistics combined with indicator kriging (snesim + ik) One problem with simulation and multiple grids is the loss of local precision on assigning original data to the nearest nodes in the coarser grid. In order to overcome this problem the results for indicator kriging were introduced in the nodes of the coarser grid as a second source of information and the simulation was performed under a hypothesis of conditional independence between the information sources. Ortiz and Deutsch (2004) propose this combination of the information provided by indicator kriging and by mps for the simulation of a porphyry copper mine. In our case we consider that the hypothesis of conditional independence is a powerful one and we only want to avoid a loss of local precision in projecting the data onto the coarser grid. The same simulation procedure was used: 50 realizations on a point support and integration in blocks included in mine planning. 6. Results analysis Fig. 4 depicts a transversal and longitudinal representation of the deposit simulated using the different techniques. The longitudinal profile corresponds approximately to the east–west axis, with the direction of the lineation produced by the cleavage–bedding intersection shown in the form of a slight pitch of the slate struc-
55
Fig. 4. Sections of the continuous model of the deposit showing useful slate (dark grey) and non-useful slate (light grey). Above: conditional simulation model. Centre: multiple-point statistic model. Below: multiple-point statistic model reinforced with ordinary kriging.
tures. The transversal profile, running north–south, clearly defines how the useful slate structures are delimited according to the direction marked by the kink bands. It can be observed that the structures that shape the useful slate follow the anisotropy direction of the variogram, which is coherent with the geological interpretation of the deposit. It is evident that the models simulated with mps compared to sisim show areas of useful slate with more regular and continuous contours (greater connectivity and less entropy, after Journel and Deutsch, 1993). This is a consequence of using the training image and an n-dimensional template to simulate values for the unsampled points; unlike the conditional simulation, this leads to realizations with short-scale variations/noise since the simulation is based on information obtained solely from the pairs of points. These noisy features often seem to be geologically unrealistic and make the realizations less credible to geologists (Deutsch, 1998). Table 2 shows the overall estimate for grade mean and standard deviation, for all the blocks in the deposit to be extracted according to the mining plan. The standard deviation for the mean for the ik model was calculated from the overall variance for the kriging. For the simulated models the results corresponded to the values obtained for the 50 simulations. Established as the reference for the mean grade was the ik model, given that this was the best unbiased linear estimator (BLUE). Note that for all the cases the standard deviation for the mean was around 4%. The snesim model, however, gives the value that is farthest from the mean with respect to ik, which we
56
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
Table 2 Statistics for the 4 techniques used to evaluate useful slate
Ik Sisim Snesim Snesim + ik
Grade mean (%)
Standard deviation (%)
49.3 51.1 46.1 48.1
4.1 3.7 3.7 4.7
The standard deviation for ik was calculated from the overall variance for the estimation.
believe to be due to the loss of local precision arising from the projection of the data in the coarsest grid (as has been commented earlier). For this reason the snesim + ik model was included in the overall estimate of the mean grade (48.1%), at the cost, however, of increasing the standard deviation (4.7%). Of interest is Fig. 5, which represents the grade distribution function for the blocks obtained using the different methods. For the simulated models the results correspond to the composition of the 50 histograms from the percentiles. Indicator kriging shows significant clustering around the mean (a small degree of dispersion), with a smoothing of the results. This is because the estimation aims at local accuracy so as to reduce the dispersion of the results in a magnitude that is approximately equal to the mean kriging variance of the support used. Basically, small values are overestimated, whereas larger values are underestimated. Sisim, which has a distribution function that is practically uniform, increases the dispersion of the grades. Although the sequential indicator simulation aims at reproducing the two-point variability of the real underlying phenomena, when large-range continuity is present, the traditional variogram-based simulation approaches do not provide good n-point reproduction. The algorithm tends to produce realizations with high entropy/low connectivity, and, in our case, this is further aggravated by the practically uniform distribution of the starting data. Combining the results for the points so as to calculate the grade for each block obtains a uniform grade distribution, as the blocks are large relative to the—noisy—simulated structures (that is, the information for the blocks contains a lot of noise). Finally, snesim and snesim + ik show more clustering at the extremes of the distribution function (a large degree of dis-
Fig. 5. Distribution functions of useful slate frequencies obtained via ik, sisim and mps.
persion). Note that the greater results variability for the snesim + ik application arises from the hypothesis of conditional independence. Unlike the sisim, snesim takes into account the moments of an order greater than 2, thereby reflecting the continuity pattern for the training image. The realizations clearly differentiate between large zones of useful slate and large areas of non-useful slate. Grouping the information in the blocks, we observe that a high proportion have very high or very low grades, with an increased dispersion in the grade distribution. With snesim + ik this effect is multiplied on originating useful slate structures of larger dimensions although fewer in number. These two models better reflect the reality of the mine, in which the miner typically encounters zones of spoil (poor in useful slate) or of mineral (rich in useful slate); in other words, it is not normal to encounter areas with grades of around 50%. Furthermore, bearing in mind that the distribution function for the point information (the drill cores) is represented by 50% of the data with values equal to 0 (waste slate) and the remaining 50% with values equal to 1 (useful slate), intuitively it can be deduced that the effect of grouping the information in a larger support would be precisely that represented by the distribution functions obtained with the mps algorithms. Fig. 5 shows how introducing the cutoff grade into each of the distribution functions affects the finances of the operation. Since the cutoff grade is low (19%), the total percentage of useful slate recovered from the blocks after diamond-wire cutting will not vary in relative terms by more than 3% from one model to another. What will be very different, however, is the total cost of perforation and diamond-wire cutting, given that the percentage of blocks that surpass the grade varies greatly between the models generated by the geostatistical techniques and using mps. This is where the reason for using a fictitious cost and cutoff grade in the mine plan for the krigged block model becomes evident—it is a way of reducing financial bias in the mine evaluation. In fact, as mentioned previously, approximately 60% of the orebody is cut, a value which is very close to that obtained with the mps technique. The difference between the models, in fact, is even smaller than can be deduced from the graph, given that a higher cutoff grade than the marginal grade of around 25% is used in the quarry. Fig. 6 depicts the omnidirectional variogram for the different estimated and simulated block models (for the latter a model
Fig. 6. Omnidirectional variograms for the kriging and simulated models.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
considered as representative was chosen). What is clear from the graph is that the fact that a simulation technique may be capable of reproducing the variogram of the points does not mean that it can reproduce the variogram on a different support when the simulated data is grouped and averaged. Nor should it be overlooked that although the simulation techniques theoretically permit the calculation of an infinite number of equiprobable realizations of the statistical model that reflect reality, the model/reality approximation will only be effective in the terms previously dictated by the model, whether approximation to the mean, the histogram, variogram function, or n-point statistics. Paradoxically, of the different techniques used in the variogram obtained, snesim was the closest to the geostatistical model, with a sill that practically coincided with the theoretical variance in dispersion calculated from the theoretical variogram. The effect of introducing the ik model information in some of the grid points of the snesim algorithm is again evident, with a sill that is substantially greater than the theoretical sill; also reflected is the lesser dispersion of sisim and, logically, of ik. These results are coherent with the arguments described above in relation to the dispersion of the distribution functions for the block grades, and, to a certain extent, confirm that the training image used is appropriate. It should be recalled that the effect of introducing the conditional independence hypothesis is to increase the size/continuity of the structures. This explains the greater range and the greater plateau for the variograms obtained from the realizations of snesim + ik. Table 3 shows deviations in the total production of saleable slate over the life of the mine—taking as reference the model ik— and also the coefficient of variation. In this case (although this cannot be established conclusively), the combination snesim + ik seems to adapt better to reality if we take into account that the ik model overestimates somewhat the total production by using an estimated model and a fictitious cutoff grade. Due to the bias produced by the projection of the data, snesim results in excessively low production. If we graph the mean annual deviation in production for each model with respect to the ik model, it will be seen that only snesim + ik locates the values at around zero (below zero in fact). Given that each year some 150 blocks of slate are extracted, this would seem to indicate that snesim + ik does not introduce local bias (at least for this time interval or number of blocks). An interesting exercise was the calculation of the distribution function for the discounted cash flows with a view to quantifying how uncertainty as to knowledge of the deposit affects uncertainty as to the profitability of the mining operation. This requires the financial model and mining sequence to be introduced into each of Table 3 Deviation in total production and coefficient of variation (CV) after applying the cutoff grade (19%)
Sisim Snesim Snesim + ik Reference model: ik.
Deviation (%)
CV
+2.7 − 7.6 − 2.9
7.6 8.3 10
57
Fig. 7. Discounted cash flows for the simulation techniques and for ordinary kriging.
the simulated deposits prior to implementing the financial analysis, this time with the original marginal cutoff grade, that is, without modifying the unit cost of perforation and diamond-wire cutting, as was done with the kriging. Fig. 7 shows the results: the X axis represents the discounted cash flows and the ordinate axis shows the percentiles deriving from the 50 simulated block models for each technique. In all cases the median is below the value obtained by kriging. This result can be interpreted as contradictory in the sisim case, as this produces a positive bias in total slate production with respect to the value obtained by kriging. However, this is not the case if we observe the grade distribution function for the blocks obtained with sisim (Fig. 5). The total cost of perforation and diamond-wire cutting assigned by sisim is practically 37% greater than the real cost (82% of the blocks will be cut in accordance with the graph, compared to the real quarry value of 60%). In contrast, the values obtained using snesim and snesim + ik are very close to reality (64% and 62% of cut blocks, respectively), and so these models will introduce no significant degree of bias in the cost. Therefore, the strong negative deviation in the value of the median presented by snesim with respect to ik is due to the bias introduced into total production (Table 3). Despite this bias that emerges with the sisim and snesim simulation techniques—displacement of the median/mean with respect to the value of ik—the variance of the distribution function for the discounted cash flows with respect to the mean in both cases is around 130%2 or, put in another way, the coefficient of variation is around 11.4%. Moreover, snesim + ik has an expected value that is quite close to that obtained with ik but with a variance value that is somewhat greater. This is due to the double effect of introducing—under the conditional independence hypothesis—certain values of ik into snesim, improving local accuracy at the cost of increasing the dispersion of the results. 7. Conclusions Traditional variogram-based geostatistics and multiple-point statistics rely on the same assumption of stationarity. Otherwise estimation or simulation would not be possible. From this
58
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59
perspective, all these techniques have the same limitations, although there is a crucial difference between them: whereas in the mps there is control over the information concerning the higher-order statistics in the geostatistical simulation that's not possible (with the exception of the case in which the sequential Gaussian simulation algorithm is used under the hypothesis that the multiple-point statistics are multivariate Gaussian). In an indirect way, the algorithm is given the power to decide over these statistics. These techniques also have the same requirement, namely, the need to construct a variogram/training image that is representative of the phenomenon being studied. The quality of the results will depend fundamentally on these tools, which are the bases for calculating the simulations. In this article we analyzed the use of different techniques for modelling the reserves in a slate deposit. The results that best reflected the reality were those obtained using multiple-point statistics, which best reproduced the structure of the deposit, both from a statistical and geological perspective. The advantage of the mps algorithms compared to the sisim algorithms is that they produce simulations with a greater connectivity between points— less entropy—and that better represent the variability in the deposit, which is even more manifest when working with nonpoint supports. The main drawback is the need to create a training image, which is also a complex matter. Comparatively it is more complex to construct a training image—which has to reflect the n-point statistics for the phenomenon—than a two-point variogram. Moreover, the variogram construction task can use adjustment techniques through regression, and this is undoubtedly very useful when working in three dimensions. The crossed validation, furthermore, enables the fit of the results to be compared when different theoretical variograms are used. The mps models provide useful slate percentages—above the cutoff grade—that are much close to reality than those based on geostatistics. This is fundamental when calculating slate reserves, costs and revenue flows, and ultimately, when evaluating the mine from a financial perspective. An interesting aspect of the evaluation of the deposit is that—assuming that the remaining parameters remain constant—the financial risk profile derived exclusively from uncertainty in regard to knowledge of the deposit has, for all the techniques, a coefficient of variation that is relatively small if we consider the high variability in the variogram function with respect to the distance and number of cores used in evaluating the studied area. Returning to the evaluation of the deposit using mps, we emphasize that the loss of local precision on using multiple grids with the snesim algorithm can be avoided by introducing part of the information supplied by ik into the coarser grid used in the simulation, although at the cost of increasing the dispersion variance. Finally, it should be underlined that when modelling and evaluating a deposit, the ultimate aim of the mine planner is to design and plan mining in such a way as to maximize net present value (NPV)—in statistical terms, expected value. The importance of assessment and evaluation of the risk associated with mining (as mentioned earlier in this article) should not be taken for granted. Nonetheless, it should be pointed out that the current
approach to mining design requires a stochastic focus of the deposit; consequently, design based on optimization algorithms should start with simulated rather than estimated models of the deposit (Dimitrakopoulos et al., 2007; Bastante et al., in press). Given the availability of the software, computers, algorithms and models necessary to implement this focus, all we need are experts who combine multidisciplinary training with consolidated on-site experience. Acknowledgment This research has been partially supported by Grant BIA200766218 of the Spanish Ministry of Education and Science. References AENOR: Asociación Española de Normalización y Certificación, 2005. UNE 12326-1: 2005. Productos de pizarra y piedra natural para tejados inclinados y revestimientos discontinuos. Parte 1: Especificación de producto. Madrid. Bastante, F.G., Taboada, J., Ordóñez, C., 2004. Design and planning for slate mining using optimisation algorithms. Engineering Geology 73, 93–103. Bastante, F.G., Taboada, J., Alejano, L.R., Ordóñez, C., 2005. Evaluation of the resources of a slate deposit using indicator kriging. Engineering Geology 81, 407–418. Bastante, F. G., Taboada, J., Alejano, L., Alonso, E. in press. Optimization tools and simulation methods to design and evaluate a mining operation. Stochastic Environmental Research & Risk Assessment. DOI 10.1007/ s00477-007-0182-6. Deutsch, C.V., 1998. Cleaning categorical variable (lithofacies) realizations with maximum a-posteriori selection. Computers & Geosciences 24 (6), 551–562. Deutsch, C.V., Journel, A.G., 1998. GSLIB. Geostatistical Software Library and User's Guide, Applied Geostatistical Software, second edition. Oxford University Press, Inc., New York. Dimitrakopoulos, R., Farrelly, C., Godoy, M., 2002. Moving forward from traditional optimisation: grade uncertainty and risk effects in open pit mine design. Transactions of the IMM, Section A Mining Industry 111, A82–A89. Dimitrakopoulos, R., Martinez, L., Ramazan, S., 2007. A maximum upside / minimum downside approach to the traditional optimization of open pit mine design. Journal of Mining Science 43 (1), 73–82. Gómez-Hernández, J., 1991. A stochastic approach to the simulation of block fields conductivity upon data measured at a smaller scale. Ph.D. thesis, Stanford University, Stanford. Goovaerts, P., 1997. Geostatistics for Natural Resources Evaluation. Oxford Univ. Press, New York, USA. Guardiano, F., Srivastava, R.M., 1993. In: Soares, A. (Ed.), Multivariate Geostatistics: Beyond Bivariate Moments, Geostatistics-Troia, vol. 1. Kluwer Academic Publications, Dordrecht, pp. 133–144. Journel, A., 1974. Geostatistics for conditional simulation of ore-bodies. Economic Geology 69, 673–687. Journel, A.G., Huijbregts, C.J., 1978. Mining Geostatistics. Academic Press, London. Journel, A.G., Deutsch, C.V., 1993. Entropy and spatial disorder. Mathematical Geology 25 (3), 329–355. Journel, A., 2002. Combining knowledge from diverse sources: an alternative to traditional data independence hypotheses. Mathematical Geology 34 (5), 573–596. Liu, Y., 2006. Using the Snesim program for multiple-point statistical simulation. Computers & Geosciences 32, 1544–1563. Ortiz, J.M., Deutsch, C.V., 2004. Indicator simulation accounting for multiplepoint statistics. Mathematical Geology 36 (5), 545–565. S-GeMS, The Standford Geostatistical Modeling Software: Product information on the Internet. Available online at: http://sgems.sourceforge.net (accessed July 2007). Strebelle, S.B., 2000. Conditional simulation of complex geological structures using multiple-point statistics. Mathematical Geology, 34, 1–22.
F.G. Bastante et al. / Engineering Geology 98 (2008) 50–59 Strebelle, S.B., Journel, A.G., 2001. Reservoir modeling using multiple-point statistics. SPE International 71324, 1–11. Taboada, J., Ordóñez, C., Saavedra, A., Fiestras-Janeiro, G., 2006a. Fuzzy expert system for economic zonation of an ornamental slate deposit. Engineering Geology 84, 220–228.
59
Taboada, J., Matías, J.M., Araújo, M., Ordóñez, C., 2006b. Assessing the viability of underground slate mining by combining an expert system with a GIS. Engineering Geology 87, 75–84. Tran, T., 1994. Improving variogram reproduction on dense simulation grids. Computers & Geosciences 20, 1161–1168.