Yerbilimleri, 31 (1), 45–52 Hacettepe Üniversitesi Yerbilimleri Uygulama ve Araştırma Merkezi Dergisi Journal of the Earth Sciences Application and Research Centre of Hacettepe University Geostatistical grade modeling of Choghart North anomaly iron ore deposit through disjunctive kriging Choghart kuzey demir cevheri yatağı anomalisi için ayırıcı kriging yöntemiyle jeoistatistiksel tenör modellemesi Ali Akbar DAYA, Hossein HASSANI Faculty of Mining and Metallurgy, Amirkabir University of Technology, Tehran, IRAN Geliş (received) : 01 Aralık (December) 2009 Kabul (accepted) : 02 Şubat (February) 2010 ABSTRACT This paper is devoted to the application of the disjunctive kriging method in the Choghart north anomaly iron ore deposit in Yazd province, Iran. The local distributions of the values of a regionalized attribute at unsampled locations can be assessed by disjunctive kriging. The case study consists of borehole samples measuring the iron concentration. A Gaussian isofactorial model is fitted to these data and disjunctive kriging was used to assess the local probabilities that the actual concentrations exceed a threshold value and to divide the ore into an economic and uneconomic parts. The tools and concepts are complemented by a set of computer programs and applied to the case study. The study showed that disjunctive kriging can be applied successfully for modeling the grade of an ore deposit. Keywords: Choghart north anomaly, disjunctive kriging, geostatistics, nonlinear kriging. ÖZ Bu çalışma, İran’ın Yazd bölgesinde bulunan Choghant kuzeyi demir cevheri anomalisine ayırıcı kriging yönteminin bir uygulamasını içermektedir. Bölgeselleştirilmiş bir değişkenin, örneklenmemiş noktalardaki değerlerinin yerel dağılımları ayrıca kriging ile kestirilebilmektedir. Sunulan çalışmada sondaj örneklerinde yapılan demir derişimleri kullanılmıştır. Bu verilere Gauss izofaktoryel model uydurulmuş ve gözlenen derişimlerin belirli bir eşik değerden yüksek olma olasılıklarının kestirilmesi ve cevheri ekonomik ve ekonomik olmayan kısımlara ayırabilmek amacıyla ayırıcı kriging yöntemi kullanılmıştır. Bu çalışma kapsamında kullanılan araç ve kavramlar bir dizi bilgisayar programı yardımıyla uygulanmıştır. Çalışmanın sonuçları, ayırıcı kriging yönteminin cevher tenörlerinin modellenmesinde başarıyla uygulanabildiğini göstermiştir. Anahtar Kelimeler: Choghant kuzey anomalisi, ayırıcı kriging, jeoistatistik, doğrusal olmayan kriging. H. Hassani E-mail: hhassani@aut.ac.ir 46 Yerbilimleri INTRODUCTION The aim of much drill hole sampling and mapping is to enable miners and their advisors to predict values of ore deposit properties at sites that have not been sampled. Kriging in its various forms is one way of doing this, providing estimates of minimum and known variance, and there is now a large literature on its applications in earth science. Miners must decide for each block from its sampling estimate whether to extract it for processing (if its concentration exceeds the economic threshold) or send it to waste, and for this they want to know the probability that the true value exceeds the threshold. To solve this problem Matheron (l976) developed the method known as ‘disjunctive kriging’. The main advantage of disjunctive kriging over the simpler techniques, whether of kriging or prediction from a classification, is in providing these probabilities. These probabilities should enable a miner or his advisor to assess the risks associated with imprecise estimates. The probability that a critical value is exceeded depends on the distribution function. For a normal distribution these probabilities can be read from published tables or calculated numerically using standard procedures. Disjunctive kriging involves transforming data to a normal distribution and then determining for each point of interest the probability that the true value exceeds the threshold. Yates et al. (1986a, 1986b) have provided a full derivation of Matheron’s method and illustrated it with examples in hydrology. Other geostatistical methods used to calculate the probability that the true value exceeds the threshold are conditional simulations (Chile’s and Delfiner, 1999), multigaussian kriging (Verly, 1983; Emery, 2005) or non-parametric estimators such as indicator kriging and its flavors (Deutsch and Journel, 1998). Disjunctive kriging has been applied mostly in the bivariate Gaussian model (Rendu, 1980; Yates et al., 1986a, 1986b; Yates and Yates, 1988; Webster and Oliver, 1989, 2001; Wood et al., 1990; Oliver et al., 1996; Chica-Olmo and Luque-Espinar, 2002). This paper is intended to investigate and to evaluate the potential and applicability of disjunctive kriging as a tool for modeling the grade of the Choghart North Anomaly Iron Ore Deposit. All of the statistical and geostatistical calculations and graphical output generated for this case study were made using the software system implemented by the authors. DISJUNCTIVE KRIGING Principles of the disjunctive kriging of blocks and estimating reserves were given in detail in earlier works (Chile’s and Delfiner, 1999). Basically, in disjunctive kriging, the variable z0 to be estimated is decomposed into a sum of disjoint (uncorrelated) components of sample values. When kriging of the separate components is possible, the procedure is feasible, i.e., when the joint probability density function of z0 (or the transformed Yo) and each sample za (or Ya) is of isofactorial type. In practice, a continuous variable like the iron grade of an ore deposit can always be transformed by anamorphosis into a Guassian equivalent Y, and then only a joint Guassian hypothesis for the probability density function (PDF) of samples and blocks is required. We start by transforming the measured variable, Z(x), to one Y(x) that has a standard normal distribution such that Z ( x ) = Φ[Y ( x )] (1) This is done using Hermite polynomials, which are related to the normal distribution by Rodriguez’s formula. H k ( y) = 1 d k g( y ) k k !g( y) dy (2) Where; g(y) is the normal probability density function, k is the degree of the polynomial taking values 1, 2,… and (1/(k!)0.5 is a standardizing factor. The first two Hermite polynomials are as follows. H o ( y) = 1 , H 1 ( y ) = −y ; (3) Daya and Hassani Thereafter the higher order polynomials obey the recurrence relation. H k ( y) = − k −1 1 yH k −1 ( y ) − H k −2 ( y ) k k The Hermite polynomials are orthogonal with respect to the weighting function (-y2/2) on the interval from −∞ to +∞ ; they are independent components of the normal distribution of ever increasing detail. Many functions of Y(x) can be represented as the sum of Hermite polynomials. f {Y ( x )} = fo H o {Y ( x )} + 1 H 1 {Y ( x )} + f2 H 2 {Y ( x )} +  , (5) Because the polynomials are orthogonal, the coefficients required for Eq. (1) can be calculated as given below. Z ( x ) = Φ[Y ( x )] = Φ0 H 0 {Y ( x )} + Φ1 H1 {Y ( x )}Φ2 H 2 {Y ( x )} +  = ∞ ∑ Φ H {Y (x)} k k=0 k 2 2 ∑λ ik cov[ H k {Y ( x j )}, H K {Y ( xi )}] = cov[ H k {Y ( x j )}, H k {Y ( xo )}] ; for all j (6) (7) So from n points in the neighborhood of x0 estimations of the Hermite polynomials are: (9) Or alternatively the following expression can be written. n ∑λ ρ i =1 ik k ( xi − x j ) = ρ k ( x j − x0 ); for all j (10) In particular, the procedure enables us to estimate Z(x0) by Zˆ ( x0 ) = Φ{Yˆ ( x0 )} = Φ0 + Φ [ Hˆ k { y( x )}] + 1 1 0 Φ [ Hˆ { y( x )}] + k 0 2 2 (11) The kriging variance of Hˆ k {Y ( x )} is given by the following expression. n The transform is invertible, and so the results can be expressed in the original units of measurement. To krige the variable of interest, Z(x), we simply krige the Hermite polynomials separately and sum their estimates to give the disjunctive kriging estimator. Zˆ DK ( x ) = Φ0 + Φ1 Hˆ 1k {Y ( x )} + Φ Hˆ k {Y ( x )} +  n i =1 (4) 47 σk2 = 1 − ∑ λik ρ k ( xi − z0 ) (12) i =1 The disjunctive kriging variance of fˆ[Y ( x0 )] is given below. ∞ σ ( x0 ) = ∑ fk2σk2 ( x0 ) 2 DK (13) k =1 Once the Hermite polynomials have been estimated at x0, the conditional probability, that the true value there exceeds the critical value, Zc can be estimated. The transformation Z(x) =F[Y(x)] means that zc has an equivalent yc on the standard normal scale. Since the two scales are monotonically related their indicators are the same. n Hˆ kk {Y ( xo )} = ∑ λik H k {Y ( xi )} i =1 (8) The next step is inserting them into Equation 7. λik are the kriging weights, which are found by solving the simple kriging equations. [ Z ( x ) ≤ zc ] = [Y ( x ) ≤ yc ] (14) [Y ( x )  yc ] , which is the complement of [Y ( x ) ≤ yc ] , the kth Hermite coefficient is as below. For 48 Yerbilimleri fk = ∫ +∞ −∞ yc H k ( y )g( y ) dy ∫ CASE STUDY [ y ≤ yc ]H k ( y )g( y ) dy  −∞ (15) The coefficient for k=0 is the cumulative distribution to yc. f0 = G ( yc ) (16) And for larger k fk = 1 H k −1 ( yc )g( yc ) k (17) The indicator can be expressed in terms of the cumulative distribution and the Hermite polynomials: Ω [Y ( x ) ≤ yc ] = G ( yc )  ∞ 1 H k −1 ( yc )g( yc ) H k {Y ( x )} ∑ k k =1 (18) Its disjunctive kriging estimate is obtained by as follow. ˆ DK Ω [ y( x0 ) ≤ yc ] = G ( yc )  L 1 H k −1 ( yc )g( yc ) Hˆ kk { y( x0 )} ∑ k k =1 (19) Where, L is some small numbers. The kriged k estimates Hˆ k { y( x0 } approach to zero rapidly with increasing k, and so summation need extend over only few terms. This is the same as ˆ DK Ω [ z( x0 ) ≤ zc ] . In this instance, we are interested in the probability of excess, and so we compute. Ω̂ DK [ z( x0 )  zc ]  Ω̂DK [ y( x0 )  yc ]  1 Hk −1( yc )g( yc ) Hˆ kk { y( x 0)}. (20) 1 − G ( yc ) −∑ k k 1 L The deposit under scrutiny in the case study is an iron deposit, which is in the northwest of the famous Choghart mine in the Central Iran Zone (Figure 1). The deposit is composed of magnetite-hematite with varying amounts of alkali amphiboles. It occurred in large low-grade orebodies. There is no general consensus regarding the origin of this iron oxide deposit. Some authors (Moor and Modabberi, 2003) believe that it has been formed directly from magmas filling volcanic diatrems or flowing as lavas, while others (Samani, 1988) suggest metasomatic replacement of preexisting rocks by hydrothermal (deuteric) solutions charged with iron that leached from cooling felsic plutons. The orebody and the metamorphosed country rock are cut by several diabasic dikes. Hematite is the second most ubiquitous mineral after magnetite. Although some primary hematite is also found in the drill cores, most of the hematite is secondary in origin. Some goethite and hydrous iron oxide occur on the surface, but disappear rapidly with increasing depth. Calcite, dolomite, secondary hematite and talc occur throughout the orebody as veinlets and cementing material of oxidized ore. Rutile and goethite are probably the results of a total transformation of the earlier formed martite. According to field observations, the host rocks are metamorphic rocks, hornblendite–pyroxenite, volcanic rocks (rhyolite, andesite) and the pyroclastic fill of the vent. Principally, this deposit was explored by 26 boreholes (Figure 2). The data file gives the name of each drill, the coordinate of drills, the grade of each element, the measure depth and azimuth of each drill, the inclination of boreholes and level, and lithology coding etc. In general, the drilling grid is irregular (see Figure 2). As the bench height for mining is fixed at 12 m, borehole samples were regularized at 12.5 m intervals. Daya and Hassani Figure1. 49 Figure 1. The location of the North Anomaly of Choghart deposit in central Iran. Şekil 1. Choghart cevheri kuzey anomalisinin Orta İran’daki konumu. Figure2. APPLICATION Data Transformation Figure 2. Borehole sample location map of North Anomaly of Choghart deposit. Şekil 2. Choghart cevheri kuzey anomalisinden alınan örnek yerleri. Usually, the random field that represents the regionalized attributes under study does not have any univariate distributions, and therefore must be transformed into a field with a distribution suited to a known isofactorial model. In the Gaussian case, this procedure is known as normal scores transformation. The problem Figure3. 50 Yerbilimleri in this case is that the variable is not normally (Gaussian) distributed, and therefore, it has to be transformed from the original distribution into a standard Gaussian distribution, the process being known as a Gaussian anamorphosis transformation. The histogram of the iron concentrations (Figure 3) is not normal. Under these conditions, as suggested by Chile`s and Delfiner (1999), a Guassian model is preferred. In a first attempt, the data are transformed to values with a standard normal distribution. As, to give an example from the data of the authors, the raw distribution of Fe is shown in Figure 3, and Figure 4 shows the transformed Gaussian distributions. Variography Figure2. Taking Figure2. spatial variety and randomicity into account, the variogram function can reflect the spatial variable structure of a regional variable. Regarding the variogram analysis of the normal scores data, Figures 5 and 6 of the variogram in different directions do not allow one to detect an anisotropy, so only an omnidirectional model is considered. The model consists of a pure nugget effect with 0.32 plus a spherical scheme with sill 0.90 and range 150 m. This model is required since disjunctive kriging estimation will be based on it as well as on the block anamorphosis calculated before. Figure3. Figure3. Figure5. Figure 3. Raw histogram of borehole sample iron concentrations. Şekil 3. Kuyu örneklerine ait demir derişim yatay sütun diyagramı. Figure 5. Variogram in different azimuths forFigure Chaghart 6. North anomaly. Şekil 5. Choghart cevheri kuzey anomalisinde farklı azimutlar için varyogramlar. Figure4. Figure4. Figure 7 Figure 4. Transformed histogram of borehole sample iron concentrations. Şekil 4. Kuyu örneklerinde demir derişimlerine ait dönüştürülmüş histogram. Figure 6. Variogram in different dips for Chaghart North anomaly. Şekil 6. Choghart cevheri kuzey anomalisinde farklı eğimler için varyogramlar. Figure 8 Daya and Hassani Estimation by Disjunctive Kriging To verify the practical usefulness of disjunctive kriging to estimate block values, the theory has been applied to the Choghart north anomaly iron ore deposit. The estimate of the grade of the survey samples only is not enough. The fact that the grade of a block is estimated at 10%, e.g., does not mean that this grade, as it will be known at the time of extraction, cannot exceed a threshold of 20%. So forecasting the probability of the Y(x) past the threshold of 20% is significant. Consider that in this case study 20% is cut off. It means that more than 20% is ore and lower is waste. Figures 7 and 8 show the grade and estimated error variance maps computed by disjunctive kriging. According to these two figures, a decision can be made. Figures 9 and 10 give the probability and probability error variance maps for a threshold of 20 %. CONCLUSIONS 51 distributions in the values of a spatial attribute (Fe in this case), using disjunctive kriging instead of simple or ordinary kriging. Implementation of disjunctive kriging is fairly simple when done within the scope of a bivariate isofactorial model. In practice, given the wide range of possibilities in the modeling, trial-and-error procedures are often necessary to choose a type of bivariate distribution and fit its parameters. This choice may be based on the experience of the practitioners and their understanding of the data. From the case study, we concluded that disjunctive kriging can be used to model the uncertainty of mapping iron ore concentrations in an ore deposit. It is hoped that this example, taken from very different application fields will encourage practitioners to apply disjunctive kriging with models that extend the capabilities of the bigaussian model. We end by saying that if the data are relatively marginally skewed and the goal is to predict nonlinear functional, then disjunctive kriging would be suitable. This work addressed the problem of assessing the uncertainty and estimating the conditional Figure 7 Figure 7 Figure 7. Estimated block averages of iron by DK. Şekil 7. Ayırıcı kriging yöntemiyle tahmin edilen demir ortalamaları. Figure 8. Estimated block averages error variance by DK. Şekil 8. Ayırıcı kriging yöntemiyle tahmin edilen ortalamaların hata varyansları. Figure 9. Conditional probabilities that block averages of iron are more than 20%. Figure 8 Şekil 9. Figure Demir ortalamalarının %20’den büyük oldu8 ğu koşullu olasılıklar. Figure 10. Conditional probabilities error variance that block averages of iron are more than 20%. Şekil 10. Demir ortalamalarının %20’den büyük olduğu hata varyanslarının koşullu olasılıkları. Figure 9 Figure 9 F 52 Yerbilimleri REFERENCES Chica-Olmo, M., and Luque-Espinar, J.A., 2002. Applications of the local estimation of the probability distribution function in environmental sciences by kriging methods. Inverse Problems, 18 (1), 25–36. Chile’s, J.-P., and Delfiner, P., 1999. Geostatistics: Modeling Spatial Uncertainty. Wiley, New York, NY. Emery, X., 2005. Simple and ordinary multigaussian kriging for estimating recoverable reserves. Mathematical Geology, 37 (3), 295–319. Matheron, G., 1976. A simple substitute for conditional expectation: the disjunctive kriging. In: Applied Geostatistics in the Mining industry, M. Guarascio, M. David and C. Huijbregts (eds.), Reidel, Dordrecht., pp. 221-236. Moor, F., and Modabberi, S., 2003. Origin of Choghart iron oxide deposit, Bafq mining district, Central Iran: new isotopic and geochemical evidence, Journal of Science, Islamic Republic of Iran, 14 (3), 259-269. Oliver, M.A., Webster, R., and McGrath, S.P., 1996. Disjunctive kriging for environmental management. Environmetrics, 7 (3), 333–357. Rendu, J.M., 1980. Disjunctive kriging: comparison of theory with actual results. Mathematical Geology, 12 (4), 305–320. Samani, B.A., 1988. Metallogeny of the Precambrian in Iran. Precambrian Research, 39, 85-105. Verly, G., 1983. The multigaussian approach and its applications to the estimation of local reserves. Mathematical Geology, 15 (2), 259–286. Webster, R., and Oliver, M.A., 1989. Optimal interpolation and isarithmic mapping of soil properties: VI. Disjunctive kriging and mapping the conditional probability. Journal of Soil Science, 40, 497–512. Webster, R., and Oliver, M.A., 2001. Geostatistics for Environmental Scientists. Wiley, Chichester. Wood, G., Oliver, M.A., and Webster, R., 1990. Estimating soil salinity by disjunctive kriging. Soil Use and Management, 6 (3), 97–104. Yates, R.W., Arricka,. W., and Myersd, .E., 1986a. Disjunctive kriging: I. Overview of estimation and conditional probability. Water Resources Research, 22, 615-621. Yates, R.W., Arricka,. W., and Myersd, .E., 1986b. Disjunctive kriging. II. Examples. Water Resources Research, 22, 623-630. Yates, S.R., and Yates, M.V., 1988. Disjunctive kriging as an approach to management decision making. Soil Science Society of America Journal, 52 (6), 1554–1558.