GEOSTATISTICAL INTERPOLATION AND SIMULATION OF RQD MEASUREMENTS by Yang Yu A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF APPLIED SCIENCE in THE FACULTY OF GRADUATE STUDIES (Mining Engineering) THE UNIVERSITY OF BRITISH COLUMBIA (Vancouver) April 2010 © Yang Yu, 2010 ii Abstract The application of geostatistics to strongly skewed data has always been problematic. The ordinary geostatistical methods cannot deal with highly skewed data very well. Multi-Indicator methods are potential candidates for the interpolation of this type of datasets, but the workload associated with them is sometimes intimidating. In this study, six geostatistical estimators, namely, ordinary kriging, simple kriging, universal kriging, trend-only kriging, lognormal kriging and indicator kriging, as well as two deterministic estimating techniques – inverse distance weighted and nearest neighbor – were applied to a highly negatively skewed RQD dataset to determine which one is more appropriate for interpolating a geotechnical model, based on their summary statistics. Universal kriging was identified to be the best method, with the trend being modeled by a quadratic drift item and the residual being estimated with simple kriging. Two simulators – sequential Gaussian and sequential indicator – were applied in this thesis for the purpose of identifying the probabilistic distribution at each location and joint-distribution for multiple locations. A transfer function was defined to transform each of the realizations to a single conceptual cost value for the sake of risk analysis, based on the relation between RQD and tunnel support specifications summarized by Merritt (1972). The distribution of the cost values corresponding to all of the realizations are quasi-normal and no estimators except KT produced a cost value that was less than two standard deviations away from the mean, once again proving the smoothing effect of all linear weighted estimators. If high values account for a dominant percentage in the sample dataset, the estimated RQD models are likely to be over- pessimistic compared with their simulated counterparts. GSLIB was used in this thesis, which in its original form does not impose a limit on the maximum number of samples coming from the same borehole when performing either estimation or simulation. It is recommended to customize the source code of GSLIB to implement this constraint to check what effect it would have on the results, if some commercial FORTRAN compiler can be made handy. iii Table of Contents Abstract ...................................................................................................................................... ii Table of Contents...................................................................................................................... iii List of Tables .............................................................................................................................. v List of Figures ........................................................................................................................... vi List of Nomenclature .............................................................................................................. viii Acknowledgements ................................................................................................................... ix 1 Introduction ....................................................................................................................... 1 1.1 Background............................................................................................................. 1 1.2 Objectives ............................................................................................................... 4 1.3 Literature Review ................................................................................................... 4 1.4 Outline of Thesis..................................................................................................... 4 2 Exploratory Data Analysis................................................................................................ 6 2.1 Geological Setting .................................................................................................. 6 2.2 Drill Core Data ....................................................................................................... 7 2.2.1 Data Verification .................................................................................................... 9 2.2.2 Geological Database ............................................................................................... 9 2.3 Compositing.......................................................................................................... 10 2.4 Trend Analysis...................................................................................................... 15 3 Geostatistical Interpolation ............................................................................................ 16 3.1 Introduction........................................................................................................... 16 3.2 Sample Variogram Calculation............................................................................. 17 3.3 Model Fitting ........................................................................................................ 20 3.4 Cross Validation ................................................................................................... 28 3.5 Kriging.................................................................................................................. 33 3.5.1 Simple Kriging – SK ............................................................................................ 35 3.5.2 Ordinary Kriging – OK......................................................................................... 38 3.5.3 Trend Only Kriging – TOK .................................................................................. 40 3.5.4 Universal Kriging – KT ........................................................................................ 42 3.5.5 Lognormal Kriging – LK...................................................................................... 45 3.5.6 Indicator Kriging – IK .......................................................................................... 49 iv 3.5.7 Comparisons and Discussion................................................................................ 52 4 Deterministic Interpolation ............................................................................................ 56 4.1 Inverse Distance Weighted – IDW....................................................................... 56 4.2 Nearest Neighbor – NN ........................................................................................ 58 4.2.1 Comparisons and Discussion................................................................................ 59 5 Simulation......................................................................................................................... 62 5.1 Sequential Gaussian Simulation ........................................................................... 63 5.1.1 Introduction........................................................................................................... 63 5.1.2 Response Value Distribution ................................................................................ 65 5.2 Sequential Indicator Simulation ........................................................................... 69 5.2.1 Introduction........................................................................................................... 69 5.2.2 Response Value Distribution ................................................................................ 70 5.2.3 Post Processing ..................................................................................................... 71 6 Conclusions and Recommentations................................................................................ 74 6.1 Simple Kriging vs. Ordinary Kriging ................................................................... 74 6.2 Ordinary Kriging vs. Universal Kriging............................................................... 74 6.3 Summary............................................................................................................... 75 6.4 Estimation or Simulation ...................................................................................... 76 6.5 Recommendations................................................................................................. 77 7 Bibliography..................................................................................................................... 78 v List of Tables Table 2-1 Composite Summary Statistics.................................................................................. 14 Table 3-1 Experimental Variogram Directions ......................................................................... 19 Table 3-2 Model Sample Variograms - Parameters .................................................................. 23 Table 3-3 Cross-Validation Parameter File ............................................................................... 29 Table 3-4 Summary Statistics for the Error Distributions of the Three Cases .......................... 30 Table 3-5 Comparison of Statistics of True and Estimated Values........................................... 30 Table 3-6 Grid Setup Parameters............................................................................................... 34 Table 3-7 Anisotropic Variogram Models in Reduced Form.................................................... 51 Table 3-8 Summary Statistics for the Estimates and the True Values ...................................... 53 Table 3-9 Summary Statistics for the Error Distributions ......................................................... 54 Table 4-1 Summary Statistics for the Estimates and the True Values ...................................... 59 Table 4-2 Summary Statistics for the Error Distributions ......................................................... 60 Table 4-3 Summary Statistics for All of the Grid Cell Values.................................................. 60 Table 5-1 Comparison of RQD and Support Requirements for a 6m-wide Tunnel .................. 65 Table 5-2 Hypothetical U/G Support Unit Costs Used in the Transfer Function...................... 65 Table 5-3 Support Costs Calculated by Estimators and Conditional Simulators ...................... 67 Table 5-4 Probability Distribution @Cutoff = 25 ..................................................................... 72 Table 5-5 Probability Distribution @Cutoff = 50 ..................................................................... 72 Table 5-6 Probability Distribution @Cutoff = 75 ..................................................................... 73 Table 5-7 Probability Distribution @Cutoff = 90 ..................................................................... 73 vi List of Figures Figure 1-1 Various Views of an OK Estimated Geotechnical Model ......................................... 2 Figure 1-2 Positively Skewed Distribution of WO3 grades......................................................... 3 Figure 1-3 Procedure for Measurement and Calculation of RQD (After Deere, 1989.) ............. 3 Figure 2-1 Property Location Map (Wardrop Engineering Inc., 2009)....................................... 6 Figure 2-2 N-S Cross Section Looking West (Wardrop Engineering Inc., 2009)....................... 7 Figure 2-3 Scatter Plot of RQD vs. Down-hole Depth................................................................ 8 Figure 2-4 Histogram and Statistics of Interval Lengths of RQD Field Measurement ............. 10 Figure 2-5 Ordinary Compositing Procedure ............................................................................ 11 Figure 2-6 Histogram and Statistics of Customized RQD Composites..................................... 13 Figure 2-7 Histogram and Statistics of Naïve RQD Values ...................................................... 13 Figure 2-8 Lognormal Probability Plot and Statistics of Customized RQD Composites.......... 14 Figure 2-9 Trend Curve Superimposed on Scatterplot of RQD vs. Elevation .......................... 15 Figure 3-1 Sample Points Pairing Controlling Bandwidths Plot ............................................... 18 Figure 3-2 Variograms in 20 Different Directions .................................................................... 21 Figure 3-3 Down-the-hole Variogram Plot................................................................................ 24 Figure 3-4 Fitted Models of the 3 Directional Variograms of the 1 st Structure ........................ 26 Figure 3-5 Fitted Models of the 3 Directional Variograms of the 2 nd Structure........................ 26 Figure 3-6 Scatter Plot of Cross-validated Values vs. True Values .......................................... 31 Figure 3-7 Distribution of the Cross Validation Errors ............................................................. 32 Figure 3-8 Scatter Plot of Cross Validation Errors vs. Estimated Values ................................. 32 Figure 3-9 Boreholes Relative to the Grid................................................................................. 34 Figure 3-10 Frequency Distribution of the RQD Values Estimated by SK .............................. 37 Figure 3-11 Maps of Three Orthogonal Slices – SK ................................................................. 38 Figure 3-12 Frequency Distribution of OK Estimated RQD Values......................................... 39 Figure 3-13 Maps of Three Orthogonal Slices – OK ................................................................ 40 Figure 3-14 Frequency Distribution of TOK Estimated Trend ................................................. 41 Figure 3-15 Maps of Three Orthogonal Slices – TOK.............................................................. 42 Figure 3-16 Models of the 1 st Structure for KT Residuals ........................................................ 43 Figure 3-17 Models of the 2 nd Structure for KT Residuals ....................................................... 44 Figure 3-18 Frequency Distribution of KT Estimated RQD Values ......................................... 44 vii Figure 3-19 Maps of Three Orthogonal Slices – KT................................................................. 45 Figure 3-20 Histogram and QQ-Plot of the Log-transformed Composites ............................... 46 Figure 3-21 Fit Ranking of the Log-transformed Dataset ......................................................... 46 Figure 3-22 Fitted Models of the 3-D Variograms of the 1 st & 2 nd Structures of LK ............... 47 Figure 3-23 Frequency Distribution of LK Estimated RQD Values ......................................... 48 Figure 3-24 Maps of Three Orthogonal Slices – LK................................................................. 49 Figure 3-25 Frequency Distribution of IK Estimated RQD Values .......................................... 51 Figure 3-26 Maps of Three Orthogonal Slices – IK .................................................................. 52 Figure 3-27 Illustration of Calculation of True Values of Blocks Containing Samples ........... 53 Figure 4-1 Frequency Distribution of IDW Estimated RQD Values ........................................ 57 Figure 4-2 Maps of Three Orthogonal Slices– IDW ................................................................. 57 Figure 4-3 Frequency Distribution of NN Estimated RQD Values........................................... 58 Figure 4-4 Maps of Three Orthogonal Slices – NN .................................................................. 59 Figure 5-1 Fitted Models of the Variograms of the 1 st & 2 nd Structures of Simulation ............ 64 Figure 5-2 Distribution of 100 Sequential Gaussian Simulation Realizations .......................... 66 Figure 5-3 QQPlots of Realizations 27 and 56 of SGSIM vs. KT Estimated Values ............... 67 Figure 5-4 YZ-Slices of 2 SGSIM Realizations @441650m and KT Estimated Values .......... 68 Figure 5-5 Distribution of 100 Sequential Indicator Simulation Realizations .......................... 70 Figure 5-6 QQPlots of Realizations 51 and 11 of SISIM vs. KT Estimated Values ................. 71 Figure 5-7 YZ-Slices of 2 SISIM Realizations @441650m...................................................... 71 viii List of Nomenclature CDF: Cumulative Distribution Function CCDF: Conditional Cumulative Distribution Function GSLIB: Geostatistical Software Library IDW: Inverse Distance Weighted IK: Indicator Kriging IQR: Interquartile Range KT: Kriging with a Trend LK: Lognormal Kriging MAE: Mean Absolute Error MSE: Mean Squared Error NN: Nearest Neighbor OK: Ordinary Kriging PDF: Probability Distribution Function QQ: Quantile – Quantile RF: Random Function RQD: Rock Quality Designation RV: Random Variable SGSIM: Sequential Gaussian Simulation SISIM: Sequential Indicator Simulation SK: Simple Kriging SSE: Sum of Squared Error TOK: Trend only Kriging ix Acknowledgements First of all, I would like to express my thankfulness to my advisor, Dr. Scott Dunbar, for giving me the opportunity to carry out this research and guiding me throughout the process. It is through studying with him that I got exposed to geostatistics. I am very grateful to MR. Britt Reid and Mr. Dave Tenney, Chief Operating Officer and Exploration Manager of North American Tungsten Corp. Ltd. respectively, who generously provided me with practical engineering data based on which this research was carried out. I appreciate the guidance of Dr. Bern Klein who offered me the access to an ongoing research project, and the software support of Dr. Michael Hitch. Finally, I would like to thank my parents for their selfless support and encouraging confidence in me. 1 1 Introduction 1.1 Background Geostatistical analysis is critical in the development of a good understanding of many mining projects. It can provide important information about the orebody that can affect the value and economic risks associated with the project and be decisive in determining the most profitable mining methods. Geostatistics is intended to be used to estimate the values of spatially distributed properties on the assumption that these properties, subject to the controlling random processes, will vary with directions and distances. These properties are modeled as regionalized variables. A 3-D geotechnical model has application for many major mining ventures that require a detailed comprehension of the variability in rock mass conditions. The development of 3-D geotechnical models can provide certain guiding information for the mining phase. Using the model, one can evaluate mining areas to predict the rock mass quality. Blast and processing plant design can be derived from the rock mass quality predictions. This information can also be used for mine-life scheduling and evaluation, cost analysis, production optimization, as well as pit slope design. The application of such a model could lower costs and improve efficiencies by utilizing drilling, blasting and mining equipment according to the rock mass conditions. A precise budget expenditure estimate cannot be realized without a comprehensive knowledge of the rock mass variability. A 3-D geotechnical model can be created by using the same set of geological model exploration boreholes. In addition to equipment selection and plant design mentioned above, the information derived from the model can be used for underground layout design and prediction of support requirements that depend on rock mass conditions. The geotechnical model estimated using ordinary kriging and the dataset provided by North American Tungsten Corporation Ltd. (NATCL) is shown in Figure 1-1. The transparent view of the model was created by setting those blocks with a RQD value falling somewhere between 80% and 100% transparent. Such 3-D views were created by transporting output from GSLIB to another application equipped with a graphical user interface. 2 Figure 1-1 Various Views of an OK Estimated Geotechnical Model The measured RQD values provided by NATCL are from 28 boreholes that are located on the west side of the deposit area. Most of the samples are about 3m long and the distribution of RQD measurements is highly negatively skewed. This is in contrast to most mining datasets which usually have a small fraction of high values and a large amount of small values, as demonstrated in Figure 1-2 which shows the distribution of tungsten grades of samples from a subset of the boreholes along which the RQD measurements were made. Only a selected group of boreholes were sent to the assay laboratory for grade analysis. 3 Figure 1-2 Positively Skewed Distribution of WO3 grades RQD (Rock Quality Designation, Deere et al. 1967) is defined as the percentage of the core length recovered in pieces larger than 10cm with respect to the length of the core run. The procedure for RQD measurement is illustrated in Figure 1-3. Figure 1-3 Procedure for Measurement and Calculation of RQD (After Deere, 1989.) 4 1.2 Objectives The goal of this study is to evaluate various geostatistical and deterministic estimation algorithms applied to a geotechnical dataset having an extremely negatively skewed distribution in order to identify the estimator which returns the best summary statistics and which, in comparison with the distribution of response values produced by the realization transfer function, returns the lowest but realistic tunnel support cost. Linear-weighted algorithms applied to positively skewed datasets, typical of most mining datasets, are likely to overestimate the data while on the contrary, with negatively skewed datasets, distributions tend to be underestimated. Attempts will be made to check whether non-linear algorithms, such as lognormal kriging, are a better candidate than their linear counterparts. 1.3 Literature Review Earth sciences events have spatial, temporal, or spatiotemporal variabilities. In addition to proposing the methodology of analyzing and modeling the spatial variability possibly existing in an earth sciences sample dataset, Sen (2009) described in detail how to obtain the PDF of RQD by way of Monte Carlo simulation, based on the given PDF of intact length. Both the independent intact length formulation and the dependent intact length formulation are examined. It was discovered that the simulated RQD values deteriorated by 15% with each increment in the number of fractures. In the article titled “three-dimensional rock mass characterization for the design of excavations and estimation of ground support requirements” by Cepuritis, collected in (Villaescusa and Potvin, 2004), it was suggested that the geotechnical data needs to be in the order of at least 35% - 40% of total geological resource data to develop a robust 3-dimensional block model that accounts for local variations in rock mass conditions and that can be used to develop site specific excavation dimensions and estimate local rock reinforcement and ground support requirements. 1.4 Outline of Thesis Chapter One provides background information on the application of geostatistics in the mining industry and the value of a 3-dimensional geotechnical model. 5 Chapter Two describes the geological settings of the study area, raw sample data verification and validation, the compositing method customized to account for the non-additivity of RQD values, and the spatial characteristics of the dataset. Chapter Three describes the application of each of the geostatistical estimators to the dataset. Six methods are proposed – simple kriging, ordinary kriging, universal kriging, trend only kriging, lognormal kriging, and indicator kriging. Summary statistics of the interpolation results are provided at the end of this chapter for the sake of comparison. Chapter Four describes the performance of two deterministic estimators applied to the same dataset – inverse distance weighted and nearest neighbor. Similarly, Summary statistics of the interpolation results are provided for the sake of comparison of all eight methods. Chapter Five introduces the application of two simulators to the dataset – sequential Gaussian and sequential indicator, and how they are used for risk analysis and evaluation of estimators. Chapter Six concludes the thesis by highlighting the key results. 6 2 Exploratory Data Analysis 2.1 Geological Setting A good understanding of the geological setting is vital to the interpretation of spatial variability. For example, qualitative information such as the azimuth and dip angles of the ore body should be respected when it is inconsistent with the direction of the continuity obtained by way of spatial analysis that is largely affected by the quality of the sampling program. The Mactung deposit owned by North American Tungsten Corporation Limited (NATCL) is the world’s largest undeveloped tungsten skarn-type deposit (Wardrop Engineering Inc., 2009). The Mactung property is located in the Selwyn Mountain Range and covers the area around Mt. Allan on the Yukon/NWT border, approximately eight kilometers northwest of MacMillan Pass. Figure 2-1 is the location map of the deposit. Figure 2-1 Property Location Map (Wardrop Engineering Inc., 2009) The study area is located in the eastern Selwyn basin that refers to a region of deep-water, off- shelf sedimentation. The rocks in the study area are part of a west-trending fold belt. Folding is tight and a narrow imbricate fault zone of southerly directed east-west trending thrust faults 7 repeats Lower Cambrian to Devonian stratigraphy. South of the imbricate belt, open to closed folds and steep faults are the dominant structures. Stratigraphy in the general deposit area trends generally E-W, and dips from 10° to 40° to the south. The axes of large folds also trend E-W and may have a shallow westerly plunge. Several ages of high-angle normal faulting, of various orientations, are known in the area. A total of 9 units have been identified in the study area. The deposit is cut and offset by numerous steeply dipping northerly trending faults. RQD itself does not account for joint orientation, continuity, etc.; it must be used with other geological input. RQD is sensitive to the orientation of joints with respect to the core axis which is in this case almost parallel to the fault system. This is probably one of the reasons why the sample RQD values are relatively high. Figure 2-2 shows the N-S cross-section of the unit sequence. In the figure, the intrusive rocks are the youngest while unit 1 is the oldest. Figure 2-2 N-S Cross Section Looking West (Wardrop Engineering Inc., 2009) 2.2 Drill Core Data In the database provided by North American Tungsten Corporation Ltd (NATCL), there are a total of 33 drill holes that have RQD measurements. 28 out of these 33 holes, that is, all of the 8 MS-series drill holes except MS-153 and MS-154 which are too short, as well as 5 MT80000- series holes that are in the vicinity of the MS-series holes, were used to make the composites and interpolate the geotechnical block model. Most of the drill holes were drilled parallel to the faults which dip to the north by 70°, almost perpendicular to the orebody slightly dipping to the south by 20°. The average spacing between these holes along the strike is about 60m and 40m to 60m along the dip direction which has more variability. A relatively dense grid like this was intended to delineate the continuity of the mineralization and will be used as a reference when determining the lag distance for the semi-variogram later on. For the most part, the NQ size core recovery was close to 100%. Figure 2-3 is a scatter plot showcasing the distribution of RQD values along the depth of various drill holes. The distribution of low RQD values is sparse whatever the depth is while high RQD values are evenly distributed along the depth. This characteristic can be used as a strong indication that the whole borehole-covered area can be treated as one population. Thus, taking features such as lithology and constraining structures into consideration will cause this thesis to lose its generality since each project has its own geological settings. Figure 2-3 Scatter Plot of RQD vs. Down-hole Depth 9 2.2.1 Data Verification The provided RQD data and assay data are in Excel format. Generally, the holes seem to have been properly located and surveyed. Almost 100% of the core has been recovered. A few discrepancies were identified: • Drill holes MS-146 and MS-165 each has a positive dip angle for one section of their drill-hole paths, perhaps due to data entry errors, leading to the borehole trace abruptly turning around while going downwards. • Some drill core intervals have a higher-than-100% recovery caused by the move of footage tag in the field. This results in RQD percentages that are more than 100%. To cater for this problem, all of these outliers were capped to 100%. • All of the drill core intervals are measured for RQD except drill hole MS156. Only half of this hole has RQD measurements. It is either because the bottom half of the core run might have been lost or because the data for this particular hole are incomplete. 2.2.2 Geological Database The data supplied is in plain text format. As stated in the previous section, RQD values were capped to 100%. Data were transferred to an MS Access database which could check and maintain the integrity of the data. Two tables are mandatory: • Collar table: contains the collar easting, northing and elevation coordinates, as well as the initial azimuths and dips. • Survey table: contains down-the-hole measurements of depths, azimuths, and dips for each drill hole, and is an interval-type table. Each drill hole can correspond to multiple records each of which represents the measurement made at certain depth. A drill hole with only two records in the table will have a straight hole path. The difference between 2 adjacent measurements will be equally spread along the corresponding interval to calculate the coordinates of the composites, as discussed in section 2.1.3. 10 Other datasets, such as assay data, geotechnical data (RQD measurements), can be saved in optional tables. All the other tables are linked to the collar table by the “hole identity” field. Identities not included in the collar table but appearing in any other table will be regarded as illegal. 2.3 Compositing Compositing refers to the process of choosing an appropriate sample interval or support and splitting the original measure length intervals into these regularized intervals. Compositing can be thought of as re-cutting the drill core into pieces of equal length. This is required so that samples will have the same weight or influence when used to calculate block values (Hustrulid and Kuchta, 2006). Sample lengths vary from 0.20m to 13.73m. About 38% of the samples are 3 meters in length; 34% 3.10m in length; 14% 3.05m in length. A total of 91.6% of drill core intervals have a length that is equal to or more than 2.25m. Figure 2-4 shows a histogram and its summary statistics depicting the distribution of sample lengths. Figure 2-4 Histogram and Statistics of Interval Lengths of RQD Field Measurement For geotechnical modeling, the sample RQD measurements must have the same length, or rather, be on the same support, to make physical or statistical distances the only variable 11 affecting the weights assigned to each sample. For additive attributes such as grade, the compositing principle, a weighted linear combination process, is illustrated by Figure 2-5. Figure 2-5 Ordinary Compositing Procedure However, according to the definition of RQD by Deere (1967), it is not appropriate to calculate the average RQD of a composite spanning two intervals this way, although it has been interpolated like grade in a few engineering practices. The property that allows the mean of some variables to be calculated by a simple linear average is known as additivity (Coward, Vann, Dunham, and Stewart, 2009). In view of the length distribution depicted by the histogram in Figure 1-1, a decision was made to customize the compositing procedure by ignoring intervals less than 2.25m and only calculating the coordinates of the centroid of longer pieces. This way, the averaging of multiple independent intervals can be avoided. Intervals longer than 3.75m will be split into as many 3m-long subsections as possible, with each subsection having the same RQD as the parent interval. Recovery rate was used as an optional weighting field to filter out those intervals which satisfy the length requirement but have a lower-than-85% recovery rate. The pseudo code of the compositing procedure is given below: 1 st 3-m composite 2 nd 3-m composite L1 L2 Suppose the grade of piece #1 is grd 1 and the grade of piece #2 is grd 2 , then the average grade of the 2 nd composite is: 3 2211 1 1 LgrdLgrd L Lgrd g n i i n i ii mean ∗+∗ == ∑ ∑ = = Core Piece #1 Core Piece #2 12 // calculate the coordinates of each 1cm-long section of each drill hole For each collar in the Collar table { Select all the survey measures corresponding to the current collar from the Survey table; For each survey measure / section { Split the down-hole length into 1cm-long pieces; Determine the quadrant of the current point to calculate increment in azimuth; Spread the differences in azimuths and dips of the 2 section ends equally to the pieces; Add the azimuth and dip increments to those of the previous point 1cm above; Calculate the coordinates of the current point using the new azimuth and dip angles; Insert the coordinates and the down-hole accumulated length into a database table; } } // calculate the coordinates of the centroid of each RQD interval satisfying the criteria For each interval satisfying the compositing criteria (length and recovery rate) { Calculate the down-hole length of the centroid of the current interval; Select the record having the same down-hole length from the table containing xyz coordinates; Insert the centroid coordinates, down-hole length between the centroid and the collar, and RQD value into another table; } The purpose of dividing the borehole path into equal length pieces and spreading the changes in azimuth and dip is to produce a smoothly curved hole trace. Shown in Figure 2-6 and Figure 2-7 are the histogram of composites generated using the customized compositing method and the histogram of naïve RQD data, respectively. As the summary statistics shown in the legends indicate, after compositing, or rather, after eliminating those short samples and poorly recovered samples, the composite dataset became more concentrated as the Interquartile range and the standard deviation both decrease, but the distribution became more negatively skewed. The mean of the distribution generated using the customized compositing method is about 4% higher than that of the original samples. The same effect would have been obtained, as shown in Table 2-1, if the conventional compositing method had been used which is just a linear averaging process. 13 Figure 2-6 Histogram and Statistics of Customized RQD Composites Figure 2-7 Histogram and Statistics of Naïve RQD Values Table 2-1summarizes the statistics of customized composites and ordinary composites, with conventional composites referring to those sections generated as shown in Figure 2-5. There are not many differences between these two statistics because both the sample lengths and the RQD sample values contain an element that accounts for a large percentage of the set. 14 Table 2-1 Composite Summary Statistics Customized Compositing Conventional Compositing Samples 2337 2544 Minimum 0% 0% Maximum 100% 100% IQR 23.3% 22.9% Mean 82.8% 80.5% Median 90.3% 89.4% Std Dev 22.4 24.2 Skewness -1.91 -1.77 Figure 2-8 shows how the customized composites will look like if log transformation has been applied. It is obvious that the transformed data do not follow a normal distribution since the lognormal probability plot even has a higher skewness and kurtosis than the original data. The original data contain 36 records that have a RQD value of zero. To check the lognormality, these zeros were replaced with 0.001. What was used here is natural logarithm. The statistics shown in the figure enclosed in a frame are about the frequency distribution of the log- transformation, but with all zero values deleted instead of being replaced by a very small number. The skewness and kurtosis have both improved obviously but are still much higher than before the log-transform. Certain kriging methods works best with an approximately normally distributed data input. Figure 2-8 Lognormal Probability Plot and Statistics of Customized RQD Composites 15 2.4 Trend Analysis Trend analysis is about collecting information and spotting a trend or pattern in the information and is done using ordinary least squares. If a trend is shown to exist in the data, then it must be removed before kriging and added back to the residual after interpolation to get the complete estimates. Figure 2-9 is similar to Figure 2-3 but with RQD on the Y axis. The fitted curve has a parabola-like shape, indicating the existence of a quadratic drift. It is recommended to model the trend as a low-order )2(≤ polynomial of the coordinates in the absence of a physical interpolation (Deutsch and Journel, 1998). The low correlation indicates that the residual component of each measured location is not trivial in scale. The trend formula identified in the figure 2-9 was used in KT kriging. Although the quadratic drift in elevation has a small coefficient, it is better to keep it since the elevations are high, rendering the quadratic drift almost as significant as the linear component. More sampling is recommended to refine the trend. Figure 2-9 Trend Curve Superimposed on Scatterplot of RQD vs. Elevation 16 3 Geostatistical Interpolation 3.1 Introduction The next step after data analysis is the quantitative modeling of the statistics of the population. Applications, such as estimation and risk analysis, rely on the model to be used. A model is nothing but a representation of reality (Goovaerts, 1997). But depending on the amount of data available, a few of models can be built to represent the unique reality. Models can be classified into two groups: • Deterministic models This kind of models assume that the estimated value for an unsampled location is the true value without recording the error - )()( uZuZ −∗ . In decision making, it is essential to know the distribution of the estimation error. • Probabilistic models Probabilistic models provide a series of possible values as well as the corresponding probabilities of occurrence by treating each location as a random variable. In the case of continuous variables, a cumulative distribution function would be set up while in the case of category variables, a probability would be given to all of the instances of a category variable. Geostatistics is primarily based on the concept of random function which can be regarded as a collection of dependent random variables. A random function is described by the distribution of the random variable at each place in the region along with the dependencies. This so-called spatial distribution law is quite demanding and is seldom needed in mining. Most geostatistical applications are based on the stationarity of the first two moments of the law, i.e., the moments do not vary spatially so that the joint distribution of any vector of random variables is location- independent. Otherwise, many realizations of the joint distribution are needed to predict the population distribution. In most cases, there is only one realization – the orebody under study. Therefore, stationarity is an essential assumption. Three types of stationarity are used in geostatistics. 17 • Strong stationarity With this stationarity, any two k-component vectors, each component separated by the same distance h, will have the same distribution. If F is the cumulative distribution function of the joint distribution of random variables kxx ,...,1 , the mathematical statement of strong stationarity is: ),...,(),...,( 11 hxhxFxxF kk ++= (3-1) • Second-order stationarity Second-order stationarity arises when strict stationarity is only applied to pairs of random variables. A random function is assumed to be second-order stationary if all the random variables have the same expected value m and the covariance 2)]()([)( mhxZxZEhC −+∗= (3-2) exists and is a function of the distance h. • Intrinsic hypothesis A random function is assumed to be intrinsic if the expected value exists and remains constant everywhere, and the variance of the difference between two random variables exists. )(2]))()([()]()([ 2 hxZhxZExZhxZVar γ=−+=−+ (3-3) It must be emphasized that stationarity is not a property of the physical dataset, but a property of random functions. 3.2 Sample Variogram Calculation The sample variogram can be calculated using the following equation: ∑ = −+= )( 1 2)]()([ )(2 1 )( hN i ii xZhxZ hN hγ (3-4) where ix are the locations of the samples, )( ixZ are the values at the locations, )(hN is the number of pairs separated by the distance h. 18 To attain a reasonable model of 3D spatial variability, 30 or so directional sample variogram have to be examined in most cases. These variograms should be in indifferent directions. In this thesis, a total of 37 variograms were calculated by incrementing the dip angle from 0° to - 90° and the azimuth angle from 0° to 330°, with each increment being 30°. The total number of directions is equal to the number of azimuth directions multiplied by the number of dip directions. For the dip angle of -90°, there is only one direction as the corresponding azimuth is 0°. Figure 3-1 illustrates the bandwidths constraining the pairing of sample points. They are 25m horizontally and 5m vertically. Note that the searching strategy shown in the figure uses a rectangle as the bottom of the cone, with the bandwidths being equal to the side lengths respectively, while in many other implementations the cone has a circular bottom with the bandwidth being equal to the radius of the circle. Therefore, the number of pairs found by the latter case is higher. Figure 3-1 Sample Points Pairing Controlling Bandwidths Plot Listed in Table 3-1 are the directions in which the experimental variograms were calculated. Vertical Bandwidth Tail Point Horizontal Bandwidth Separation Vector Angle Tolerance 19 Table 3-1 Experimental Variogram Directions Azimuth Angle Tolerance Dip Lag (m) Lag Tolerance (m) 0° 22.5° 0° 50 22.5 30° 22.5° 0° 50 22.5 60° 22.5° 0° 50 22.5 90° 22.5° 0° 50 22.5 120° 22.5° 0° 50 22.5 150° 22.5° 0° 50 22.5 180° 22.5° 0° 50 22.5 210° 22.5° 0° 50 22.5 240° 22.5° 0° 50 22.5 270° 22.5° 0° 50 22.5 300° 22.5° 0° 50 22.5 330° 22.5° 0° 50 22.5 0° 22.5° -30° 50 22.5 30° 22.5° -30° 50 22.5 60° 22.5° -30° 50 22.5 90° 22.5° -30° 50 22.5 120° 22.5° -30° 50 22.5 150° 22.5° -30° 50 22.5 180° 22.5° -30° 50 22.5 210° 22.5° -30° 50 22.5 240° 22.5° -30° 50 22.5 270° 22.5° -30° 50 22.5 300° 22.5° -30° 50 22.5 330° 22.5° -30° 50 22.5 0° 22.5° -60° 50 22.5 30° 22.5° -60° 50 22.5 60° 22.5° -60° 50 22.5 90° 22.5° -60° 50 22.5 120° 22.5° -60° 50 22.5 150° 22.5° -60° 50 22.5 180° 22.5° -60° 50 22.5 210° 22.5° -60° 50 22.5 240° 22.5° -60° 50 22.5 270° 22.5° -60° 50 22.5 300° 22.5° -60° 50 22.5 330° 22.5° -60° 50 22.5 0° 22.5° -90° 3 1.35 20 3.3 Model Fitting Given below are the definitions of some concepts relevant to variogram modeling: • Nugget effect The vertical jump caused by such factors as sampling error from the value of zero at the origin to the value of the variogram at extremely small separation distances is called the nugget effect. • Sill The plateau the variogram reaches at the range is called the sill. It is usually equal to the variance of the population. • Range The distance at which the variogram reaches the plateau is called the range. The directional variograms mentioned above just describe the variability along the anisotropy axes, to perform kriging, a 3D variogram model is needed. One method for combining the various directional models into an equivalent isotropic model is to reduce all directional variograms to a common isotropic model with a standardized range of 1 (Isaaks and Srivastava, 1987). For geometric anisotropy, the reduced distance corresponding to the actual vector h is 2221 )()()( z z y y x x a h a h a h h ++= (3-5) where zyx hhh ,, are the components of the vector h along the three anisotropy axes, zyx aaa ,, are the ranges of the structure along the three axes. If each directional model consists of more than one structure, then each structure should have a reduced distance calculated with Equation 3-5. The generalized equation for 3D anisotropic model composed of n structures plus one nugget is ∑ = += n i ii hwhwh 1 100 )()()( γγγ ni ,...,1= (3-6) where iw is the contribution of the thi structure, ih is the thi reduced distance 21 There are two types of anisotropies – geometric and zonal. In geometric anisotropy, all directional models are assumed to have the same sill values. All the directional models must consist of the same amount and the same types of structures. A zonal anisotropy where the sill value changes with direction while the range remains constant seldom appears alone. Shown in Figure 3-2 are the experimental variograms along twenty different directions. It seems that all of the variograms level off at a variogram value of 500, namely the variance of the samples, but varying with different ranges. Some variogram nodes, especially those calculated with larger lag distances and varying significantly, might not be reliable if the number of data pairs for each of them is too small. Figure 3-2 Variograms in 20 Different Directions 22 In this study, one nugget and two spherical structures were used to model the three anisotropic directional sample variograms. One spherical model was used to describe the short-range behavior while the other model to depict the long-distance feature of the random function model. Among the legitimate mathematical models, the spherical model, together with the exponential model, is the most common variogram model. • The spherical model is a transition model which means that it flattens out and reaches the sill at the range, but before reaching the sill, it has a quasi-linear behavior. It has a standardized equation given by:     − = 1 )(5.05.1 )( 3 a h a h hγ ah ah ≥ < (3-7) • The exponential model is a transition model as well. It reaches its sill asymptotically at the effective range. Its standardized equation is given by: ) 3 exp(1)( a h h −−=γ (3-8) • The equation for the nugget effect is given by    = 1 0 )(0 hγ otherwise h 0= (3-9) This model is indeed a step function which has a value of zero at the origin but can be significantly larger than zero at small separation deviation. The common practice is to express it as a constant. The sum of it and the contributions of other structures is equal to the sill which in turn is equal to the variance )0(C . Table 3-2 contains the details of the two structures which have independent anisotropies, rendering the modeling more accurate and flexible. 23 Table 3-2 Model Sample Variograms - Parameters Parameters Values Estimator Correlogram Sill 1 Nugget 0.282 Spherical #1 Spherical #2 Contribution 0.392 0.326 Rotation around X-axis -43.6° 8.3° Rotation around Y-axis 39.1° -17.4° Rotation around Z-axis -5.1° -3.5° Range along rotated X 48m 77m Range along rotated Y 138m 241m Range along rotated Z 21m 401m Anisotropy Axes Azimuth Dip Azimuth Dip Rotated X 114° -27° 89° 17° Rotated Y 355° -44° 356° 8° Rotated Z 45° 34° 242° 71° SSE 0.0107 The sill value is one because the correlogram was used as the estimator which has a value of one according to its definition: ∑ = +−+−+ −== )( 1 )()()()( /)( )( 1 )0( )( )( hN i hhhhhii mmZZ hNc hc h σσρ (3-10) where )0(c refers to the variance of the dataset and is also equal to the sill of the model and the tail mean )( hm − is given by: ∑ = − = )( 1 )( )( 1 hN i ih Z hN m (3-11) and the head mean is similarly given by: ∑ = ++ = )( 1 )( )( 1 hN i hih Z hN m (3-12) and )()( , hh +− δδ are the standard deviations of the tail and the head respectively. 24 The nugget was determined using the down-the-hole variogram shown in Figure 3-3. Figure 3-3 Down-the-hole Variogram Plot The nugget value is equal to the intercept of the fitted model with the Y axis. Since most of the drill holes dip to the north, the down-the-hole variogram is in fact a directional variogram. No tolerance angles and bandwidths were used. Only samples from the same drill hole were paired. A final variogram value is the result of averaging the values from various holes supposing their lag distances are identical. Normally, the anisotropy axes do not coincide with the dataset axes due to rotation. The determination of the dataset axes is arbitrary while the orientations of the anisotropy axes are determined by the spatial variability. Coordinates transformations by rotation are not unique. Each estimation or simulation algorithm uses a different rotation convention. The variogram model must be using the same convention in order for the kriging algorithms to calculate the variogram correctly. In this thesis, the LRR convention for Z, X, Y axes was utilized, with L referring to the left-hand rule and R the right-hand rule. The order of the axes specified here cannot be altered, which means that the Z axis must be rotated first, followed by X and by Y. The three ranges of the anisotropy axes correspond to the maximum, medium, and minimum continuities of the phenomenon. 25 The rotation angles about various axes for each structure can be used in coordinate transformation. hRh ∗=' (3-13) The transformation matrix R is given by (Isaaks and Srivastava, 1987):           ∗−∗− − ∗∗ = )cos()sin()sin()sin()cos( 0)cos()sin( )sin()cos()sin()cos()cos( φφαφα αα φφαφα R (3-14) whereα is the clockwise rotation angle about the Z axis, φ is the clockwise rotation angle about the new Y axis caused by the first rotation. A left-hand-determined azimuth angle should change its sign since the clockwise direction is defined by using the right hand rule. The azimuth and dip angles shown in the “Anisotropy Axes” frame are the final locations of the rotated anisotropy axes for each structure and are measured in the original data coordinate system. SSE is short for total sum of squared error. This is a parameter used to evaluate how close the regression line is to the lag points. Figure 3-4 and Figure 3-5 show the fitted models of the two structures. As is obvious, the semi-variograms along some directions are quite erratic. In the first structure, it is more continuous along the Y axis of the anisotropy coordinate system since its range is the longest. Corresponding to this, the range of the first structure of the black curve line in Figure 3-4 is the highest. It should be noted that the azimuth and dip angles shown in the plots are the ones that are most identical to those listed in Table 3-2. As we only created directional models for the multiples of 30°, there is not a variogram plot, for example, in the direction Azimuth=242°, Dip=71°, the closest direction is Azimuth=60°, Dip=-60°. The calculation process is demonstrated below. The dip angle must be equal to or less than 0°.    °−≈ °≈ ⇒    °−= °= ⇒    °−= °+°= ⇒    °= °= 60 60 71 62 71 180242 71 242 Dip Azimuth Dip Azimuth Dip Azimuth Dip Azimuth 26 Displayed in the legend of Figure 3-4 are the closest azimuths and dips angles for the orientation of the anisotropy axes of the first structure found in Table 3-2, while what are displayed in the legend of Figure 3-5 are for the second structure. Figure 3-4 Fitted Models of the 3 Directional Variograms of the 1 st Structure Figure 3-5 Fitted Models of the 3 Directional Variograms of the 2 nd Structure 27 Their corresponding model equations are: )(326.0)(392.0282.0)30,240( )(326.0)(392.0282.0)30,30( )(326.0)(392.0282.0)30,120( 3.941.19 3.2185.41 9.1372.43 zz yy xx hSphhSphDipAz hSphhSphDipAz hSphhSphDipAz ++=°−=°= ++=°−=°= ++=°−=°= γ γ γ (3-15) )(326.0)(392.0282.0)60,60( )(326.0)(392.0282.0)0,180( )(326.0)(392.0282.0)30,270( 4.32960 4.2422.30 6.8523 zz yy xx hSphhSphDipAz hSphhSphDipAz hSphhSphDipAz ++=°−=°= ++=°=°= ++=°−=°= γ γ γ (3-16) The reduced distances 21 ,hh are calculated as follows:             ••                 =           = ' 1, ' 1, ' 1, 1 1, 1, 1, 1 21 1 00 0 138 1 0 00 48 1 z y x z y x h h h R h h h h (3-17)             ••                 =           = ' 2, ' 2, ' 2, 2 2, 2, 2, 2 401 1 00 0 241 1 0 00 77 1 z y x z y x h h h R h h h h (3-18) where ' , ' , ' , ,, iziyix hhh )2,1( =i are the components of vectors defined in the data coordinate system, 21 ,hh are the two reduced distance vectors, 21 ,RR are the two transformation matrices defined with equation 3-19 and 3-20.           °°∗°−°∗°− °°− °°∗°°∗° = )1.39cos()1.39sin()1.5sin()1.39sin()1.5cos( 0)1.5cos()1.5sin( )1.39sin()1.39cos()1.5sin()1.39cos()1.5cos( 1R (3-19)           °−°−∗°−°−∗°− °°− °−°−∗°°−∗° = )4.17cos()4.17sin()5.3sin()4.17sin()5.3cos( 0)5.3cos()5.3sin( )4.17sin()4.17cos()5.3sin()4.17cos()5.3cos( 2R (3-20) The values used in Equation 3-15 to Equation 3-20 are from Table 3-2. 28 The length of 1h is given by: 2 1 2 1, 2 1, 2 1,1 ])()()[( zyx hhhh ++= (3-21) and that of 2h is given by: 2 1 2 2, 2 2, 2 2,2 ])()()[( zyx hhhh ++= (3-22) The anisotropic auto variogram model in reduced form is given by: )(326.0)(392.0282.0)( 2111 hSphhSphh RangeRange == ++=γ (3-23) 3.4 Cross Validation Cross-validation is used to estimate how good the model is. By comparing the estimated values with the true values, the best approach can be determined. Considering that in most cases the exhaustive data sets are not available, the cross-validation technique was devised to make use of the sample data set to determine the better search strategy and the better variogram models. The estimation errors from different models or strategies are good indicators of whether improvements are necessary, as shown below. This technique functions by temporarily removing the sample value at a location each time and predicting a value for the same location using the remaining samples. This procedure is repeated for all samples. The resulting predicted values are then compared to the measured values to evaluate the previous decisions on anisotropy and search strategy. The criteria used to evaluate the re-estimation scores are: • The distribution of the cross validation errors },,1),()({ * Niuzuz ii L=− should be symmetric, with a zero mean and a minimum spread (Deutsch and Journel, 1998). • The scatter plot of the errors versus the estimated values should be centered on the zero error line. This is a property called “conditional unbiasedness”. The plot should have an equal spread, that is to say, the error variance should not correlate to the magnitude of the true values. This is another property called “homoscedasticity of the error variance” (Deutsch and Journel, 1998). 29 Table 3-3 shows the parameters required by the kt3d program of GSLIB to perform cross- validation. The rotation convention of GSLIB is “ZXY – LRR”. The kt3d program is capable of performing not only cross validation but also jackknife and various kriging interpolations, be it linear or non-linear. Table 3-3 will be re-used later on when kriging the block model, with modifications to certain parameters. Table 3-3 Cross-Validation Parameter File START of PARAMETERS: ./rqd.dat 0 1 2 3 4 0 -1.0 1.0e21 1 jackknife.dat 1 2 3 4 0 1 RQD_Cross-validation.dbg RQD_Cross_validation.out 60 441440 6 70 7017370 6 70 1700 5 1 1 1 2 8 0 315 126 126 0(9.6) 0(47.2) 0(-7.6) 1 0 0 0 0 0 0 0 0 0 0 0 external_drift.dat 4 2 0.282 1 0.392 -5.1 -43.6 39.1 138 48 21 1 0.326 -3.5 8.3 -17.4 241 77 401 -Data file in Geo-EAS format -X, Y, Z, variable columns -Trimming limits -1=cross validation -File with jackknife data -Debug level -nx, xmn, xsize -ny, ymn, ysize -nz, zmn, zsize -Discretization -min, max data for kriging -max per octant -Search radii -Search angles -ikrige 0=SK, 1=OK -Drift -File with drift -# of structures and nugget -Spherical structure #1 -Spherical structure #2 Several variations of the parameter file were used as inputs to the kt3d program by varying the interpolation methods and the search angles. The purpose is to identify the better estimator and search strategy. The one that has the settings given in Table 3-3 returns the best statistics – lowest error mean and mean squared error. Shown in the Table 3-4 are the summary statistics 30 of errors of the three cases. The first data column corresponds to the case when the search ellipsoid is aligned with the shared-by-all-structures anisotropic axes; statistics given in the second column are for the case when all the angles are reset to zero; the last one summarizes the results when simple kriging was used. As can be seen, the statistics of both ordinary kriging cases are better than those of simple kriging estimated values. In all of the cases, the medians are different from the means, indicating the existence of skewnesses. Shown in Table 3-5 is a comparison of statistics of the true values and the estimated values. Although the dispersions of the two OK cases are somewhat higher with higher variances and higher IQRs, their statistics are more similar to those of the true values. Moreover, the first two cases are more correlated with the true if more effective numbers had been used. For the sake of comparisons of various estimators, the same isotropic search strategy was utilized, namely, all of the search angles were reset to zeros. Table 3-4 Summary Statistics for the Error Distributions of the Three Cases Statistics OK OK – Zero Search Angles SK (m=78.9%) Number of Samples 2337 2337 2337 Mean Error 0.19 0.21 0.10 Standard Deviation 14.23 14.23 14.27 min -55.35 -55.35 -53.14 Median -1.61 -1.67 -2.16 IQR 11.85 11.77 11.91 max 91.71 91.71 91.08 Skewness 0.99 0.97 1.10 MAE 9.57 9.58 9.69 MSE 202.54 202.52 203.46 Table 3-5 Comparison of Statistics of True and Estimated Values Statistics True OK OK – Zero Search Angles SK (m=78.9%) Number of Samples 2337 2337 2337 2337 Mean 82.84 83.03 83.05 82.93 Standard Deviation 22.40 17.54 17.52 16.77 min 0 4.97 4.97 8.97 IQR 23.3 17.53 17.35 16.76 Median 90.3 89.34 89.36 88.98 max 100 100 100 99.29 Correlation 0.77 0.77 0.77 31 MAE and MSE are two statistics that incorporate both the bias and the spread of the distribution of error. Their equations are given by: Mean Absolute Error = MAE = ∑ = n i r n 1 1 (3-24) Mean Squared Error = MSE = ∑ = n i r n 1 21 (3-25) Figure 3-6 is a scatter plot of the measured values versus the estimated values. The slope of the regression line is less than 1 because kriging tends to underestimate the large values and overestimate the small values. The dispersion on the lower true value side is larger than that of the right side, visually proving the aforementioned property of kriging. Figure 3-6 Scatter Plot of Cross-validated Values vs. True Values Figure 3-7 shows the error distribution of the OK case where all of the search angles are zero. It has a quasi-normal shape. The expected value is almost zero and the standard deviation is relatively small with more than 95% of the errors falling in the range (-28, 28). The distribution perfectly abides by the first criterion mentioned above. Figure 3-8 shows the scatter plot of the errors versus the estimated values. As shown in the figure, the plot is almost centered around the zero line, indicating that there is no global bias but that the conditional bias persists since the lowest estimated values tend to be too low while the highest ones tend to be too high. 32 Conditional unbiasedness is difficult to achieve in practice. There are always overestimations or underestimations for some ranges of data. Figure 3-7 Distribution of the Cross Validation Errors Figure 3-8 Scatter Plot of Cross Validation Errors vs. Estimated Values 33 3.5 Kriging Kriging refers to a group of generalized least-squares regression algorithms (Goovaerts, 1997). All types of kriging methods associate a probabilistic value to an unsampled location since they are based on the assumption of a statistical model. Kriging methods rely on the concept of autocorrelation, a special case of correlation, which summarizes the tendency of one variable’s residual values at different locations. All kriging methods are variants of the linear regression estimator given by: ∑ = ∗ −+= )( 1 )]()()[()()( un umuZuumuZ α αααλ (3-26) where )(uαλ is the weight assigned to datum )( αuz . )(um and )( αum are the expected values of the random variables defined for each location. Both the number of samples and the weights are location-dependent. The common purpose of all kriging methods is to minimize the estimation variance )(2 uEσ given by: )}()({)(2 uZuZVaruE −= ∗σ (3-27) But depending on how the trend of the random function model is modeled, a couple of kriging variants can be proposed: • Simple kriging (SK) assumes that the expected value )(um is constant everywhere. • Ordinary kriging (OK) reduces the area or volume in which the mean is assumed to be static within the neighborhood of the unknown location. The mean value has to be dynamically calculated for each neighborhood. • Kriging with a trend (KT) models the drift with a linear combination of functions )(uf k of the coordinates. The most commonly seen functions )(uf k are the 2 nd -order polynomial functions. ∑ = ∗= K k kk ufum 0 )()( α (3-28) where kα is constant but unknown. 34 In this chapter, six kriging variants were examined: • Simple kriging • Ordinary kriging • Trend only kriging • Universal kriging • Lognormal kriging • Indicator kriging Table 3-6 gives the location and dimensions of the grid. The cell size is determined based on the dimensions of underground tunnels. Figure 3-9 shows all the boreholes relative to the grid defined in Table 3-6 Table 3-6 Grid Setup Parameters Original Cell Size Number of Cells Easting 441440 6 60 Northing 7017370 6 70 Elevation 1700 5 70 Figure 3-9 Boreholes Relative to the Grid 35 3.5.1 Simple Kriging – SK With the mean being constant throughout the study area, Equation 3-26 can be rewritten as: muuZuuZ un un ∗      −+= ∑ ∑ = = ∗ )( 1 )( 1 )(1)()()( α α ααα λλ (3-29) The )(un weights )(uαλ are chosen to minimize the estimation variance given in Equation 3-26, and can be calculated using Equation 3-31 given below in the matrix format: kuK =)(λ (3-30) where K is the covariance matrix of the samples contained within the neighborhood:           −− −− = )()( )()( )()(1)( )(111 ununun un uuCuuC uuCuuC K L MMM L (3-31) )(uλ is the vector of )(un weights:           = )( )( )( )( 1 u u u unλ λ λ M (3-32) k is the vector data-to-unknown covariances:           − − = )( )( )( 1 uuC uuC k un M (3-33) If the covariance matrix of samples is stable and invertible, the weight vector can be written as: kKu 1)( −=λ (3-34) Figure 3-10 shows the histogram of the estimated values classified into 30 groups. As is obvious, the spread of the distribution is very small. Figure 3-11 shows the color-scale maps of three slices in the XY, XZ, and YZ orientations respectively. The box plot is shown under the X axis (outside the whiskers at the 95% probability interval, inside the box at the 50% probability interval, vertical solid line at the median, and the black dot at the reference value 36 specified in the parameter file of the histogram generating program). The input parameters of the kt3d program are identical to those given in Table 3-3 except for the kriging type indicator and its argument - skmean: • ikrige = 0 • skmean = 78.9% The arithmetic means of the composited samples and of the original samples are 82.8% and 78.9% respectively. 78.9% is an unbiased estimator of the mean of the population since the actual number of samples is much more than theoretically necessary, in spite of the distribution followed by the samples. However, as the dataset, if reflected, follows a quasi-lognormal distribution which does not conform to the central limit theorem, this statistic, although unbiased, might not be the best estimator. Sichel (1952) derived the formula for calculating the mean of the original values with the mean and standard deviation of the logarithms. ) 2 exp( 2σ µ + (3-35) The mean calculated with equation 3-35 happened to be equal to that of the original samples, that is, 78.9%. The determination of the mean of the whole study area is critical and decision- making demanding. An average of 82.8% was simple kriged as well at the end of this chapter for the sake of comparison. The maximum search radii were set based on the known geological information and the anisotropy. The anisotropy has the maximum range along the Y axis and the ratios of major/median and major/minor are both 2.5. The azimuth, dip, and rake of the search ellipsoid are all 0°, in consistence with the LRR rotation rule. They were chosen based on the following criteria: • Setting the minimum number of sample data points can effectively prevents the algorithm from performing estimation in areas where data are sparse. Too low a setting of this parameter will render the estimation unreasonably variable. The minimum number of samples utilized in this thesis is two. • In kriging, all data points within the neighborhood will be utilized to estimate the unknown point, even if they are not correlated with the unknown point. Setting 37 the maximum number of data too high will not only worsen this problem but also smooth the estimate. On the other hand, if this parameter is set too low, the estimate would be unreliable. The maximum number of samples utilized in this thesis is eight. The purpose of utilizing a small number like this is to avoid making unrealistic assumption of stationarity over large search neighborhood. • The search radius needs to be set large enough to guarantee a stable estimate, but too large a search radius will make the estimate too smooth. Generally speaking, the search radius is a bit larger than the variogram range to allow the remote data to contribute to the calculation of local mean. As mentioned in the first chapter, all of the samples belong to the same population; therefore, the samples included in the neighborhood are relevant. A comparison of all the kriging variants will be performed at the end of this chapter by examining their error distributions and the statistics of the estimates. Figure 3-10 Frequency Distribution of the RQD Values Estimated by SK 38 Figure 3-11 Maps of Three Orthogonal Slices – SK 3.5.2 Ordinary Kriging – OK The difference between ordinary kriging and simple kriging is that ordinary kriging allow one to account for the variation of the mean by confining the stationary domain to the local neighborhood which is usually an ellipsoid centered on the location to be estimated. To put it in another way, the expected values of the random variables are only constant within the neighborhood. The unknown mean is eliminated from Equation 3-28 by forcing the kriging weights to sum to one. This is also the non-bias condition of ordinary kriging which calls for the introduction of Lagrange parameter to convert a constrained minimization problem into an unconstrained one. The set of weights that minimize the kriging variance under the condition that they sum to one should satisfy the following equation:        = =+ ∑ ∑ = = n i i n j iijj w CCw 1 1 0 1 µ ni L,1=∀ (3-36) 39 Equation 3-36 can be written in matrix format as:             =             ∗             1011 1 1 0 101 1 111 nnnnn n C C w w CC CC MM L L MMOM L µ (3-37) The covariances of each pair of samples and between each sample and the location being estimated can be read from the variogram model given in section 3.3. Figure 3-12 shows the histogram of the estimated values classified into 30 groups. Figure 3-13 shows the color-scale maps of three orthogonal slices. The input parameters of the kt3d program are identical to those given in Table 3-3 except for the kriging type indicator which was set to one. As both the histogram and the scale map show, OK estimates are more variable than SK estimates with a much higher variance and minimum value, due to the varying means. Figure 3-12 Frequency Distribution of OK Estimated RQD Values 40 Figure 3-13 Maps of Three Orthogonal Slices – OK 3.5.3 Trend Only Kriging – TOK The estimated value can be regarded as the sum of two components: the trend component and the residual. If the residual is assumed to have no correlation, namely 0)( =hCR , one can estimate only the trend component. The modified universal kriging system for trend is given by:        == ==+− = ∑ ∑ ∑ ∑ = = = = n kk m n K k k m kR m n m KT Kkufufu nufuuuCu uZuum 1 )( 1 0 )()( 1 )(* ,,0),()()( ,,1,0)()()()( )()()( β ββ β ααββ α αα λ αµλ λ L L (3-38) where the sm ')(αλ are the KT weights and the s m k ' )(µ are Lagrange parameters. Note that the covariance of the residual has been set to zero in the first summation equation. 41 This algorithm computes the least-squares fit of the trend model and removes the random component. The direct KT estimation of the trend component can be regarded as a low-pass filter that eliminates the random, high-frequency residual. Ordinary kriging is a special case of universal kriging, whose mean is constant locally but variant throughout the area. Figure 3-14 shows the histogram of the estimated values and Figure 3-15 shows the color-scale maps of three orthogonal slices. The trend was estimated by setting all of the linear, quadratic, and cross-quadratic drifts indicators to zero and trend indicator to one in Table 3-3. This setting will inform kt3d to estimate with ordinary kriging. The results seem to be more smoothed with a lower variance and a higher minimum value. Figure 3-14 Frequency Distribution of TOK Estimated Trend 42 Figure 3-15 Maps of Three Orthogonal Slices – TOK 3.5.4 Universal Kriging – KT Universal kriging is one of the non-stationary geostatistical estimators. With universal kriging, the mean is no longer assumed to be stationary either globally or within the search neighborhood and the trend component is modeled as a smoothly varying deterministic function of either the coordinates vector u or some auxiliary variables, whose unknown parameters are computed from the data: ∑ = = K k kk ufum 0 )()( α (3-39) Ordinary kriging corresponds to the case when 0=k . The mean is equal to 0a which is constant but unknown. The residual is usually modeled as a stationary random function with zero mean. The universal kriging system for trend is given by: 43        == =−=+− = ∑ ∑ ∑ ∑ = = = = n kk KT n K k RkkR KT n KT KT Kkufufu nuuCufuuuCu uZuuZ 1 )( 1 0 )( 1 )(* ,,0),()()( ,,1),()()()()( )()()( β ββ β αααββ α αα λ αµλ λ L L (3-40) Considering that the sample data were obtained from a relatively small area, a comprehensive and representative data set is necessary to extract a precise model for the deterministic component of the RF. In most cases, OK and SK work well without considering the trend. In this case, universal kriging was carried out by first simple kriging the residual computed by subtracting the trend defined in Figure 2-10 from the original sample values, and then adding the trend back to the interpolated values. In contrast to TOK, it is the residual that should be modeled with a random function. Figure 3-16 and 3-17 show the closest azimuths and dips angles for the orientation of the anisotropy axes of the two structures respectively. Figure 3-18 shows the distribution of the estimated values while Figure 3-19 shows the color-scale maps of three orthogonal slices. As shown in Figure 3-18, the variance is satisfying compared with that of simple kriging, although SK was indeed used to estimate the residual. Compared with the results of OK, KT’s estimates have a higher smoothing effect. Figure 3-16 Models of the 1 st Structure for KT Residuals 44 Figure 3-17 Models of the 2 nd Structure for KT Residuals Figure 3-18 Frequency Distribution of KT Estimated RQD Values 45 Figure 3-19 Maps of Three Orthogonal Slices – KT 3.5.5 Lognormal Kriging – LK If the data are very much skewed, a linear estimator may not be the best option, since the conditional expectation may be non-linear. Lognormal kriging is a non-linear kriging algorithm. It is applied to the logarithms of data if such a transform produce a normal distribution. All non-linear kriging algorithms are linear kriging applied non-linear transforms of the original data. The back-transform is sensible to the error since it amplifies any error by exponentiating it. This sensitivity has made direct non-linear kriging less and less popular (Deutsch and Journel, 1998). Equation 3-41 was used to back-transform the estimated values which relies on correction by kriging variance. ) 2 )( )(exp()( 2 µ σ −+= ∗∗ u uyuz OK (3-41) where )(2 uOKσ is the ordinary lognormal kriging variance andµ is the Lagrange multiplier. 46 Figure 3-20 shows the histogram of the transformed data and its normal QQ plot. The data had to be reflected before transformation by subtracting each datum from 101. This is used to convert a negatively skewed distribution to its counterpart whose tail is on the right side. Figure 3-20 Histogram and QQ-Plot of the Log-transformed Composites Shown in Figure 3-21 is a @RISK-generated ranking list of probability distributions that can be fitted to the reflected histogram. As can be seen, the Normal distribution is the best fit. Figure 3-21 Fit Ranking of the Log-transformed Dataset A three-parameter logarithm transformation with the constant being set to five seemed to be able to reduce the spike to a larger extent and return relatively better statistics, but the best fit 47 in this case is the Uniform distribution instead of Gaussian. Lognormal kriging is based on the strict assumption that the dataset is log normal, an assumption that virtually cannot be verified unless a comprehensive knowledge of the data set has been achieved. However, despite the fact that the histogram does not have a standard bell shape due to the large proportion of zero values, a grid was interpolated using the transformed data with ordinary kriging. One nugget and two spherical models were used to model the variogram of the transformed data, as shown in Figure 3-22. The two-parameter logarithm transformation was applied here. Figure 3-22 Fitted Models of the 3-D Variograms of the 1 st & 2 nd Structures of LK 48 The anisotropic auto variogram model in reduced form is given by: )(457.0)(157.0386.0)( 2111 hSphhSphh RangeRange == ++=γ (3-42) Figure 3-23 shows the distribution of the LK estimated values while Figure 3-24 shows the color-scale maps of three orthogonal slices. Although theoretically correct, nine times out of ten, LK does not perform as expected most likely due to the sensitivity of the back transform to the exponentiated kriging variance and Lagrange parameter that both depend on the variogram modeling. Some other factors seem to have to be considered than the kriging variance and the Lagrange multiplier, such as panel variance, but these ingredients are very project-specific. As shown in Figure 21, there is a high spike on the left side of the reflected histogram, indicating the presence of a large amount of censored data. In this case, the transform cannot push the spike to the center of the distribution, rendering LK inappropriate for this kind of dataset. In this case, the LK estimated average value is obviously higher than those of other estimators. As mentioned earlier, direct non-linear kriging techniques such as LK have fallen into disuse. Figure 3-23 Frequency Distribution of LK Estimated RQD Values 49 Figure 3-24 Maps of Three Orthogonal Slices – LK 3.5.6 Indicator Kriging – IK All of the estimators examined above are designed for the estimation of a mean value. IK, on the contrary, is used to model the uncertainty about any particular unknown value with a probability distribution. Such a model is dependent on the relations between the known information and the unknown. There are two approaches to estimating distributions. • The nonparametric approach that calculates cumulative probabilities at several thresholds. • The parametric approach that determines a priori a model that describes the distribution for any value. The nonparametric approach does not rely on a random function model to determine the statistics of up to the second order, but it does require appropriate assumptions of models for interpolation and extrapolation of values other than the specified thresholds. To calculate the 50 proportion of values below the threshold cv , values are transformed to a set of corresponding indicator variables )(,),(1 cnc vivi L , as shown below.     > ≤ = cj cj cj vvif vvif vi 0 1 )( (3-43) The cumulative proportion of values below any threshold can be expressed as: ∑ = ∗= n j cjjc viw n vF 1 )( 1 )( (3-44) where jw are weights assigned to each nearby sample and can be calculated using either SK or OK depending on the configuration and number of samples. The thresholds chosen are vital for the determination of the posteriori distribution. Normally speaking, five to twelve thresholds are enough, as suggested in (Isaaks and Srivastava, 1987). For this study, a total of six thresholds were applied – 25%, 50%, 75%, 90%, 95%, and 97%, each being depicted by a separate variogram model of distinct anisotropy. The linear model was utilized for both interpolation and extrapolation. The order relation correction and the bound checking are automatically overlooked by ik3d, the indicator kriging program. As mentioned above, IK is used to model the local uncertainty by calculating and maintaining a probability matrix for each grid cell, each element of the matrix corresponding to probability of the cell value being equal to or less than the threshold. To get an estimate for an unknown location, the E-type estimate, which is the mean of the ccdf, can be defined by the Stieltjes integral (Pollard, 1920) and calculated as a discrete sum: ))](|;()(|;([))(|;()]([ 1 1 1 '* nzuFnzuFznzuzdFuz kk K k kE − + = +∞ ∞− −≈= ∑∫ (3-45) where Kkzk ,,1, L= are the K thresholds specified, and max1min0 , zzzz K == + are the lower and upper bounds of the attribute. The E-type estimate is a least squares type estimate and is usually different from the OK estimate in that OK only uses the original data and the corresponding covariance model and is not able to qualify the uncertainty of estimate with just a error variance. Given in Table 3-7 are 51 the variogram models of each of the thresholds. Figure 3-25 shows the distribution of the IK estimates while Figure 3-26 shows the color-scale maps of three orthogonal slices. Table 3-7 Anisotropic Variogram Models in Reduced Form Thresholds Variogram Models 25% )(497.0)(033.0470.0)( 2111 hSphhSphh RangeRange == ++=γ 50% )(427.0)(237.0336.0)( 2111 hSphhSphh RangeRange == ++=γ 75% )(249.0)(255.0496.0)( 2111 hSphhSphh RangeRange == ++=γ 90% )(277.0)(154.0569.0)( 2111 hSphhSphh RangeRange == ++=γ 95% )(303.0)(181.0516.0)( 2111 hSphhSphh RangeRange == ++=γ 97% )(115.0)(243.0642.0)( 2111 hSphhSphh RangeRange == ++=γ Figure 3-25 Frequency Distribution of IK Estimated RQD Values 52 Figure 3-26 Maps of Three Orthogonal Slices – IK 3.5.7 Comparisons and Discussion To draw a conclusion as to which technique is better, the true cell values must be computed. In this case, the true value of each cell was calculated as the arithmetic mean of the composite points that fell within that cell. The 2,337 composite points were classified into 1,290 cells. Obviously, the sparser the sample points are, the less accurate the true values are and the more favorable it is to the deterministic methods which calculate the weight assigned to each sample based on its Euclidian distance to the location under estimation. For example, if a block only contains one sample, its true value will be equal to its estimated value in the case of the polygonal method, as illustrated in Figure 3-27. If a block contains more than one sample, its true value will be equal to the average of the samples enclosed within it. The combination of the sample length that is 3m or so and the block size that is 6m x 6m x 5m determines that there are at most two samples falling in any of the blocks. 53 Figure 3-27 Illustration of Calculation of True Values of Blocks Containing Samples The summary statistics for the estimates are given in Table 3-8 and the summary statistics of the estimation error are given in Table 3-9. Table 3-8 Summary Statistics for the Estimates and the True Values True Values SK OK TOK KT LK IK n 1,290 1,290 1,290 1,290 1,290 1,290 1,290 Mean 87.01 86.83 87.20 87.27 87.08 91.69 85.93 )(σStdev 17.75 12.06 13.19 12.55 13.19 10.53 13.05 CV 0.20 0.14 0.15 0.14 0.15 0.11 0.15 Min 0.001 19.51 14.55 19.43 14.54 25.45 21.05 1Q 83.36 83.86 83.71 83.77 83.51 89.81 82.11 Median 93.30 90.68 91.40 91.23 91.01 95.66 90.00 3Q 98.40 94.72 95.92 95.67 95.77 98.24 94.76 Max 100 99.01 100 100 99.93 100 98.60 ρ 0.83 0.82 0.74 0.87 0.80 0.84 54 Table 3-9 Summary Statistics for the Error Distributions SK OK TOK KT LK IK n 1,290 1,290 1,290 1,290 1,290 1,290 Mean -0.19 0.19 0.26 0.06 4.68 -1.08 )(σStdev 10.32 10.16 11.94 9.00 11.34 9.85 Min -42.47 -45.12 -53.06 -39.22 -31.84 -46.58 1Q -5.04 -4.47 -5.24 -4.13 -0.97 -5.72 Median -1.96 -0.96 -1.15 -0.95 1.39 -2.20 3Q 2.56 3.20 3.63 2.81 6.86 1.92 Max 68.06 67.59 76.11 55.51 86.31 66.69 MAE 6.68 6.46 7.63 5.80 6.81 6.67 MSE 106.42 103.26 142.56 80.91 150.38 98.13 It seems that universal kriging returned the best estimation followed by ordinary kriging and trend only kriging. The KT estimated values are more correlated to the true values and less skewed with a smaller IQR and CV. What is worth mentioning is that the variance of the KT estimates is almost identical to that of the OK estimates, in spite of the fact that SK was used to estimate the residuals, that is, a dynamically changing mean can counteract the smoothing effect brought about by the SK estimator to the residuals. Simple kriging, as expected, gave the worst prediction, since it assumes a much more optimistic stationarity than the other variants while the other methods either assume a constant mean within the local search neighborhood or make the mean a function of the coordinates. The color-scale map of the SK estimates looks less varied and its histogram is more skewed. GSLIB does not provide a drift item to represent the shift. Simply setting the quadratic drift item in Z produced certain amount of estimates that were out of the bounds of RQD limits. The results must be post-processed to cap the larger- than-100% estimates to 100% and to reset those negative estimates to zero. More drilling that covers a larger area on a reasonably dense grid is necessary to better define the trend. The current data do indicate certain trend existing down-the-hole. The statistics of trend only kriging are almost identical to those of ordinary kriging with the exceptions being in the minimum value and the correlation of the estimates. It seems that the trend plays a major role in the estimation while the residual, though trivial, does not have a zero covariance. Moreover, the MAE and MSE of ordinary kriging are lower than their 55 counterparts of trend only kriging, indicating that OK’s error distribution is less biased and more converged. The statistics of lognormal kriging are obviously different than those of any other estimator. It has much higher means for both the estimated values and the errors, as well as the smallest standard deviation for the estimates. It is not as ideal as expected for skewed dataset probably because of the high spike on the left side of the reflected, log-transformed histogram. Lognormal kriging does not seem to work well with censored distribution in which either the lowest or the highest legal value dominates the dataset. Indicator kriging is primarily designed for the analysis of local uncertainty, but its statistics for the E-type estimates are no inferior to those of all of the other estimators than KT. Its error distribution has the second lowest variance and MSE. The shortcoming of IK is that if median indicator kriging is not applicable, sample variograms must be calculated and modeled for each threshold. 56 4 Deterministic Interpolation Deterministic techniques perform interpolation using the measured points, based on the extent of similarity or the degree of smoothing, whereas the geostatistical techniques make use of the statistical properties of the measured points. There are two types of deterministic techniques: global and local. Global techniques make use of all of the measured data points in the study area while local techniques only utilizes those points contained within the specified neighborhood. Inverse Distance Weighted is an example of local interpolators. 4.1 Inverse Distance Weighted – IDW IDW entirely relies on the assumption that entities that are closer to one another are more likely to be similar than those that are separated further apart. The sampled points closer to the location whose value is to be estimated will be assigned more weights than those further away, that is to say, the influence a measured point can apply to the unknown point is inversely proportional to the distance raised to the power p. The inverse distance estimator is given in Equation 4-1: ∑ ∑ = == n i p i n i ip i d v d v 1 1 1 1 (4-1) The same search strategy was used for IDW as for all the kriging variants. Each block was discretized into 3 x 3 x 5 points. Figure 4-1 and 4-2 show the distribution of the IDW estimated values and the color-scale maps of three orthogonal slices. The distribution in Figure 4-1 looks more negatively skewed than that of ordinary kriging estimated values. A comprehensive comparison of all of the interpolation techniques investigated in the thesis will be made at the end of this chapter by looking at their error distributions as well as the distribution of their estimates. 57 Figure 4-1 Frequency Distribution of IDW Estimated RQD Values Figure 4-2 Maps of Three Orthogonal Slices– IDW 58 4.2 Nearest Neighbor – NN With this method, the sample value that is closest to the point to be estimated is chosen. This polygonal method is equivalent to a weighted linear combination that assigns all of the weight to the closest sample value. The shortcoming of this method is that much of the nearby information is ignored and therefore in most cases leads to the worst results compared with other techniques. Surpac developed by Gemcom was used to carry out this estimate. Figure 4-3 and 4-4 show the distribution of the polygonal method estimated values and the color-scale maps of three orthogonal slices. As is obvious, the cross-section views are very irregular and the histogram is very similar to that of the original sample distribution as the estimates are nothing but a subset of the samples. Figure 4-3 Frequency Distribution of NN Estimated RQD Values 59 Figure 4-4 Maps of Three Orthogonal Slices – NN 4.2.1 Comparisons and Discussion The summary statistics for the estimates are given in Table 4-1 and the summary statistics for the estimation errors are given in Table 4-2. For the convenience of comparison, the statistics of universal kriging estimated values are duplicated here. Table 4-1 Summary Statistics for the Estimates and the True Values True Values IDW NN KT n 1,290 1,290 1,290 1,290 Mean 87.01 87.23 87.31 87.08 )(σStdev 17.75 14.63 18.29 13.19 CV 0.20 0.17 0.21 0.15 Min 0.001 3.59 0.00 14.54 1Q 83.36 83.84 83.90 83.51 Median 93.30 91.99 93.50 91.01 3Q 98.40 96.66 100 95.77 Max 100 100 100 99.93 ρ 0.93 0.92 0.87 60 Table 4-2 Summary Statistics for the Error Distributions IDW NN KT n 1,290 1,290 1,290 Mean 0.21 0.30 0.06 )(σStdev 6.66 7.30 9.00 Min -26.96 -61.30 -39.22 1Q -2.47 0.00 -4.13 Median -0.35 0.00 -0.95 3Q 1.93 0.88 2.81 Max 43.07 41.60 55.51 MAE 4.10 3.53 5.80 MSE 44.41 53.39 80.91 The estimations by inverse distance weighted and even by nearest neighbor seem better than that of universal kriging. This is due to the methodology how the true average values were calculated. This method apparently favors IDW and NN since these two methods calculate weights based on the reciprocal of the physical distance raised to any specified power instead of the stochastic distance. The composite points within a cell always receive the highest weights, making the estimate similar to the average of the sample values enclosed in the cell. The higher the power, the more the sample points farther away from the cell will be screened. In the case of nearest neighbor, all of the weights are assigned to the sample that is closest to the cell. To prove the above analysis, the more holistic statistics of all the 294,000 cell values are given in Table 4-3. Table 4-3 Summary Statistics for All of the Grid Cell Values IDW NN OK KT SK TOK LK IK n 294,000 294,000 294,000 294,000 294,000 294,000 294,000 294,000 Mean 81.03 79.45 80.81 82.30 82.56 80.93 85.75 79.56 )(σStdev 19.18 27.36 19.70 15.73 6.63 19.44 17.98 19.07 CV 0.24 0.34 0.24 0.19 0.08 0.24 0.21 0.24 Min 2.48 0.00 0.94 6.94 15.61 1.01 0.00 11.66 1Q 76.97 74.00 76.54 76.03 81.86 76.46 83.13 74.39 Median 88.24 90.30 88.23 88.15 83.52 88.22 92.94 86.92 3Q 93.67 100 93.87 93.65 85.72 93.90 96.97 92.51 Max 100 100 100 99.93 99.01 100 100 98.60 Skew -1.60 -1.70 -1.60 -1.38 -2.40 -1.60 -2.07 -1.50 Kurtosis 1.84 1.90 1.74 1.25 9.77 1.77 3.97 1.38 61 The statistics of IDW are similar to those of OK except that the standard deviation and IQR of IDW are both smaller. With a slightly higher kurtosis, the IDW distribution is more normal than that of OK. It seems that as long as the data are not clustered, the inverse distance estimator can always produce an estimate that is as good as that provided by kriging, whether the distribution of samples is negatively skewed or not. Deterministic interpolation techniques do not provide the probabilities of all of the potential values for the unknown. Instead, they just calculate a number for each cell without providing any measure of the reliability of the estimates. In contrast, the favorable advantage of probabilistic interpolation techniques is that an error variance is associated with each estimate. However, unless the distribution of the spatial errors is Gaussian, the error variance can not be used to determine the confidence interval. In many cases, the Gaussian distribution model is assumed for the errors, which can be fully characterized by the mean and variance, but such an assumption is questionable in nature-controlled environment. For the purpose of risk analysis and decision making, an uncertainty-modeling probability distribution instead of a particular estimate should be constructed by way of non-parametric geostatistics – indicator kriging and sequential indicator simulation. The sequential indicator simulator will be discussed in the last chapter. 62 5 Simulation Interpolation methods produce one value for every location. In reality, there are many, equally probable values that could occur at each unsampled location (the true value for each unsampled location is unknown). Geostatistical simulation produces multiple maps that mimic the real phenomenon and provide these possible values, providing a basis for risk analyses, economic decision making and other estimations involving uncertainty. Just as the distribution of a random variable can be inferred from a set of simulated realizations, a random function is also depicted by its realizations denoted by{ } LlSuuz l ,...,1,),()( =∈ . A realization can be treated as a numerical model of the possible distribution in space of the variable under study. As described in Section 5.1, generating realizations is not the goal of simulation. The goal is to define some sort of single-input single-output transfer function to transform these realizations to produce a distribution of response values, such as tunnel support costs, that can be used in risk analysis and decision making (Birkhoff and Rota, 1978). As all of the realizations are equi-probable and conform to the same covariance function, the distribution of response values can reflect the probability of reproducing the spatial structures. The purpose of doing simulation in this thesis is justified by what can be achieved from a series of realizations generated: • The probability distribution of the variable at any location can be derived from the histogram of the L realizations while a single estimation can only provide a determined value, even enhanced by the calculation of error variance which does not imply any error distribution information. • The joint probability of any two locations at certain cut-off value can be evaluated by the ratio of the number of realizations having higher-than-the-cut-off values at these two locations to the total number of realizations. Correlations among the locations have been taken into account when generating the realizations, which means that this ratio is more precise than what can be attained from the independent marginal cumulative distribution functions for each location. 63 5.1 Sequential Gaussian Simulation 5.1.1 Introduction The sequential Gaussian simulation algorithm requires the input data be normal. The original data must be transformed if they are not normal. The normal score iy corresponding to the thi largest data iz is calculated as: ) 2 ( 11 −− + = iii cc Gy (5-1) where )(xG is the standard normal cdf (cumulative distribution function), )(1 cG − is the corresponding standard normal quantile, ic being the cumulative probability associated with iz . Back transformation is usually done by the simulation program and is calculated as: ))((1 ii yGFz −= where )(zF is the cdf of the original data. A realization can be generated following the steps given below, supposing that there are N nodes to be simulated: • Determine a random path along which each of the grid cells will be visited only once. • Model the ccdf of the first location conditional to the n original data, where ccdf stands for conditional cumulative distribution function. A SK or OK system has to be computed for each node to determine the mean and variance of the normal ccdf. • Draw from the ccdf a value which will becomes a conditioning datum for all subsequent drawing. Drawing a value from a distribution is realized by reading the quantile value )(pqz = that corresponds to a random number p uniformly distributed in [0, 1]: ]1,0[))(()()()( 1 ∈==→== − ppqFzFpFpqz (5-2) • Repeat the steps two and three until all of the N nodes have been assigned a simulated value. 64 Multiple realizations can be achieved by repeating the above 3 steps the specified times. The cumulative distribution of the series of realizations Llpqz ll ,....,1),( )()( == reproduces the target cumulative distribution function (CDF) if all of the sp l ')( are independent. Figure 5-1 shows the semi-variograms to be used in the simulation. They were created using the normal score values of the sample values. Figure 5-1 Fitted Models of the Variograms of the 1 st & 2 nd Structures of Simulation 65 5.1.2 Response Value Distribution It has long been found that RQD is closely related to the estimate of support requirements for rock tunnels. Summarized in Table 5-1 are the support specifications, as a function of RQD and tunnel width, proposed by Merritt and others (Bieniawski, 1989). To be in accordance with the tunnel width, the dimension of the grid cell is defined as 6m x 6m x 5m, as tunnels are 5m high or less. Table 5-1 Comparison of RQD and Support Requirements for a 6m-wide Tunnel No Support or Local Bolts Pattern Bolts Steel Ribs Deere et al. (1970) RQD 75-100 RQD 50-75 (1.5-1.8m spacing) RQD 25-50 (0.9-1.5m spacing) RQD 50-75 (light ribs on 1.5-1.8m spacing s alternative to bolts) RQD 25-50 (light to medium ribs on 0.9- 1.5m spacing as alternative to bolts) RQD 0-25 (medium to heavy circular ribs on 0.6-0.9m spacing) Cecil (1970) RQD 82-100 RQD 52-82 (alternatively 40- 60mm shotcrete) RQD 0-52 (ribs or reinforced shotcrete) Merritt (1972) RQD 72-100 RQD 23-72 (1.2-1.8m spacing) RQD 0-23 The unit costs given in Table 5-2 are defined as the support costs necessary for various RQD classifications. The transfer function, which is a simplified version, adds up all of the relative costs assigned to each of the grid cells based on their estimated or simulated RQD value. Table 5-2 Hypothetical U/G Support Unit Costs Used in the Transfer Function RQD (%) Rock Quality Relative Cost <25 Very poor 5 25-50 Poor 4 50-75 Fine 3 75-90 Good 2 90-100 Excellent 1 The distribution of the sums of all realizations is shown in Figure 5-2 which has a quasi-normal shape. The summary statistics displayed in the legend further proves its Gaussian nature, with the mean being similar to the median value, the skewness being approximately zero, and the 66 kurtosis being close to three. The best scenario corresponds to realization #27 which returns the minimum cost and therefore indicates the highest RQD values. The worst scenario corresponds to realization #56 which returns a cost most similar to the ordinary-kriging- estimated cost denoted by a vertical arrow on the histogram of Figure 5-2. As the KT estimated cost is more than two standard deviations away from the mean of the distribution, with 95% of confidence, it can be deemed as conservative to declare the study area unsafe if only the estimated populations by whatever estimators are accounted for. The economic impact on the project might be considerable since the all the estimators return too pessimistic a result. Figure 5-2 Distribution of 100 Sequential Gaussian Simulation Realizations Table 5-3 shows the costs computed by various estimators. The same search strategy was used in both estimation and simulation, although different variograms were used for some kriging variants and the conditional simulation. The lognormal kriging method returns the lowest cost most likely because the distribution input to kriging is not normal thanks to the high spike. Simple kriging seem to have the most conditional bias as indicated by its seriously smoothed scale map. Although the mean of the SK estimated values is similar to those by the other kriging variants, its cost is the highest, indicating that most of the estimates should be fluctuating around the mean. This can be verified by the SK histogram which has a very small standard deviation. The result also depends on what mean value was specified. The cost calculated by the inverse distance method is similar to those by ordinary kriging and trend only 67 kriging. This proves that so long as the data is not clustered, IDW’s performance is not inferior to OK’s performance. The lowest cost returned by universal kriging, enhanced by its better statistics, suggests that KT, or rather, the non-stationary random function model is more favorable for this kind of dataset, supposing the trend identified is reliable and conforms to the geological interpretation. Table 5-3 Support Costs Calculated by Estimators and Conditional Simulators Items Ranks Costs Lognormal kriging 1 490,905 Universal Kriging 2 548,864 Inverse Distance Weighted 3 564,672 Trend only ordinary kriging 4 565,539 Ordinary kriging 5 566,558 Nearest neighbor 6 570,423 Indicator kriging 7 590,748 Simple kriging – 78.9% 8 619,114 Realization #27 446,403 Realization #56 529,645 It is the smoothing effect of all linear weighted estimation methods that renders their RQD estimations lower than it does their simulated counterpart. The QQPlots of the two bounding realizations shown in Figure 5-3 indicate that the solid lines always fall below the diagonal line after some point, that is, the simulated values are higher than the estimated values because all of the estimators underestimate the high values. Figure 5-3 QQPlots of Realizations 27 and 56 of SGSIM vs. KT Estimated Values 68 The color-scale maps shown in Figure 5-4 indicate that the result of KT is a smooth set that has less variance than the actual data, and that the simulated images, in contrast, capture the variability of the actual data. There is certain amount of grid cells in blue color in the simulated images but not in the KT estimated image. This further demonstrates the nature of kriging of overestimating the low values. The continuity of the high value zones in the simulated images to a large extent coincides with those found in the KT estimated image. Figure 5-4 YZ-Slices of 2 SGSIM Realizations @441650m and KT Estimated Values Sequential Gaussian simulation requires simple kriging be used since the mean and variance of the conditional cumulative distribution function at each node are only mathematically equal to the mean and the kriging variance of simple kriging (Goovaerts, 1997). If ordinary kriging had been used, the realizations costs would have been higher, but they were still lower than the cost estimated by universal kriging by at least 2 standard deviations. 69 5.2 Sequential Indicator Simulation 5.2.1 Introduction Sequential Gaussian simulation requires the check of the normality of the distribution to a level of at least two points (Goovaerts, 1997). If the biGaussian assumption can not be validated, other procedures should be considered to determine the local conditional cumulative distribution function, such as the sequential indicator simulation algorithm. Very few mining- related datasets can satisfactorily meet the aforementioned condition. Sequential indicator simulation is implemented as follows: • Discretized the continuous attribute RQD into several classes using the specified K cut-off values. Then transform each datum into an indicator vector:    = 0 1 );( kzui α otherwise zuzif k≤)( α Kk ,,1L= (5-3) • Define a random path which runs through each grid cell only once. • At each node, determine the K ccdf values using whichever kriging method is appropriate - in this thesis, simple kriging was used to be in line with the Gaussian simulation; adjust any order relation deviation; draw a simulated value from the ccdf; add the simulated value to the existing dataset. • Proceed to the next node along the random path. If more than one realization is requested, the above procedure will be repeated the specified number of times. The random path can be different from iteration to iteration. Indicator kriging is a non-parametric approach where several values of the cumulative distribution function are used to estimate the complete distribution. Interpolation and extrapolation are required for values other than the defined cut-offs. For this study, the linear model was used for the extrapolations in both the lower tail and the upper tail, and for the interpolation of the middle sections of the distribution. The use of the power model for the lower tail and the hyperbolic model for the upper tail did not make any difference than the linear model in this case. 70 5.2.2 Response Value Distribution The distribution of the sums of all realizations is shown in Figure 5-5. Full kriging was performed instead of median indicator kriging which requires more system resources since a kriging system had to be set up for each threshold. The best scenario corresponds to realization #51 which returns the minimum cost and therefore indicates the highest RQD values. The worst scenario corresponds to realization #11 which returns a cost most similar to the ordinary kriging estimated cost denoted by the right-most vertical arrow on the histogram of Figure 5-5. Again the OK estimated cost is more than two standard deviations away from the mean of the distribution. The statistics indicate that the average cost computed by the sequential indicator algorithm is about 13% higher than that by the sequential Gaussian algorithm. The indicator approach deals with spatial correlation of extreme values better than the Gaussian model which does not allow for significant spatial correlation of extremely large or small values (Goovaerts, 1997). The histogram of the original sample dataset and the cross-section views of the realizations indicate that connected large clusters of RQD values exist in the deposit. Figure 5-5 Distribution of 100 Sequential Indicator Simulation Realizations Figure 5-6 shows the QQPlots of the two bounding realizations. It seems that the distribution of sisim-simulated values is more similar to that of KT-estimated values, compared with sgsim. 71 Figure 5-6 QQPlots of Realizations 51 and 11 of SISIM vs. KT Estimated Values Figure 5-7 shows the color-scale maps of the two bounding realizations. It seems that the sisim-simulated values are less variable than those simulated by sgsim. The clustered high values found in the sgsim images do not appear in the indicator simulated images. The sisim- simulated values are more concentrated with a large amount of middle-ranked values showing up in the images. Figure 5-7 YZ-Slices of 2 SISIM Realizations @441650m 5.2.3 Post Processing Post processing allows for a number of summaries to be extracted from the simulated realizations, such as the E-type estimate which is the point-by-point average of the realizations, the variance of the conditional distribution, the probability of exceeding a specified cutoff, the 72 average values above and below that cutoff, and the value where a fixed conditional cumulative function value p is reached. For this study, the 3 rd type of summaries was computed with the RQD cutoffs specified as 25, 50, 75, and 90. Listed in Tables 5-4 to 5-7 are the distribution of the probabilities of grid cell RQD values being greater than the corresponding threshold. “-999” denotes the non-existence of the values for the corresponding items. Table 5-4 Probability Distribution @Cutoff = 25 Probability Proportion Average Probability Mean > cutoff Mean < cutoff 0.0 -> 0.1 0.0187% 0 -999 8.6 0.1 -> 0.2 0.0000% 0 0 0 0.2 -> 0.3 0.0000% 0 0 0 0.3 -> 0.4 0.0014% 0.36 64.24 11.9 0.4 -> 0.5 0.0048% 0.47 68.87 12.68 0.5 -> 0.6 0.0259% 0.56 70.41 12.47 0.6 -> 0.7 0.0789% 0.66 71.98 12.42 0.7 -> 0.8 0.1963% 0.76 74.03 12.55 0.8 -> 0.9 1.6490% 0.87 78.45 12.57 0.9 -> 1.0 84.5803% 0.97 83.35 12.51 1.0 13.4449% 1 84.78 -999 Table 5-5 Probability Distribution @Cutoff = 50 Probability Proportion Average Probability Mean > cutoff Mean < cutoff 0.0 -> 0.1 0.0429% 0 -999 25.05 0.1 -> 0.2 0.0014% 0.16 87.57 20.63 0.2 -> 0.3 0.0095% 0.25 84.86 24.15 0.3 -> 0.4 0.0207% 0.35 83.48 24.49 0.4 -> 0.5 0.0602% 0.45 81.89 23.21 0.5 -> 0.6 0.1313% 0.55 82.32 24.28 0.6 -> 0.7 0.4327% 0.66 82.16 25.84 0.7 -> 0.8 1.3548% 0.76 83.33 27.25 0.8 -> 0.9 39.0337% 0.88 86.06 30.51 0.9 -> 1.0 58.1442% 0.94 87.31 31.19 1.0 0.7687% 1 89.74 -999 73 Table 5-6 Probability Distribution @Cutoff = 75 Probability Proportion Average Probability Mean > cutoff Mean < cutoff 0.0 -> 0.1 0.1102% 0 -985.53 48.89 0.1 -> 0.2 0.0412% 0.16 93 41.68 0.2 -> 0.3 0.1520% 0.25 92.86 43.51 0.3 -> 0.4 0.3901% 0.35 92.78 47.25 0.4 -> 0.5 1.2146% 0.45 93.05 49.38 0.5 -> 0.6 4.3527% 0.56 92.9 51.25 0.6 -> 0.7 21.8973% 0.66 92.36 51.1 0.7 -> 0.8 42.6136% 0.75 92.02 51.16 0.8 -> 0.9 25.1160% 0.84 91.81 50.36 0.9 -> 1.0 3.7031% 0.93 92.38 48.93 1.0 0.4092% 1 93.31 -999 Table 5-7 Probability Distribution @Cutoff = 90 Probability Proportion Average Probability Mean > cutoff Mean < cutoff 0.0 -> 0.1 0.2279% 0.01 -905.81 64.19 0.1 -> 0.2 0.1626% 0.15 96.14 51.73 0.2 -> 0.3 0.5721% 0.26 96.1 54.18 0.3 -> 0.4 3.0667% 0.36 96.14 60.19 0.4 -> 0.5 30.1463% 0.46 96.19 64.66 0.5 -> 0.6 52.9282% 0.54 96.21 66.53 0.6 -> 0.7 11.1310% 0.63 96.29 69.37 0.7 -> 0.8 1.1738% 0.74 96.58 72.72 0.8 -> 0.9 0.2827% 0.83 96.79 73.48 0.9 -> 1.0 0.0099% 0.92 96.87 73.71 1.0 0.2990% 1 96.69 -999 74 6 Conclusions and Recommentations 6.1 Simple Kriging vs. Ordinary Kriging Ordinary kriging is usually preferred to simple kriging because it requires neither knowledge nor stationarity of the mean over the entire area. It amounts to applying SK to each neighborhood, with the local mean instead of the global mean being used. The difference between the SK and OK estimates of z at u is caused by a departure of the local mean from the global mean. The relation between the Ok and SK estimators is: ])()[()()( *** mumuuZuZ OK SK mSKOK −+= λ (6-1) More precisely, since )(uSKmλ is usually positive, the OK estimate is smaller than the SK estimate in low-valued areas where the local mean is smaller than the global mean. Conversely, the OK estimate is larger than the SK estimate in high-valued areas where the local mean is larger than the global mean. In this sense, the commonly used OK algorithm with local search neighborhoods already accounts for trends (varying mean) in the z-data values. The discrepancy between the two estimates increases as the location u being estimated gets farther away from the data locations. 6.2 Ordinary Kriging vs. Universal Kriging Ordinary kriging amounts to estimating the local mean within each search neighborhood, then performing SK on the corresponding residuals. Similarly, KT amounts to estimating within each search neighborhood the trend component, then performing SK on the corresponding residuals. Both OK and KT estimates can be expressed as the sum of two terms: the same linear combination of sample data and a function of the trend estimates. ∑ ∑ = = −+= )( 1 ))(( )( 1 *** * )()()()()()( un umf un OK SK OK SK OK OK umuumuzuuz α α ααα λλ 44444 34444 21 (6-2) ∑ ∑ = = −+= )( 1 ))(( )( 1 *** * )()()()()()( un umf un KT SK KT SK KT KT umuumuzuuz α α ααα λλ 4444 34444 21 (6-3) 75 The difference between the two is thus: ∑ = −−−=− )( 1 ****** )]()()[()]()([)()( un KTOK SK KTOKKTOK umumuumumuzuz α αλ (6-4) Any difference between the OK and KT estimates originates from a difference between the two trend estimates, hence from the usually arbitrary decision of modeling the local trend within a neighborhood as a constant or a particular polynomial of order K. 6.3 Summary The use of a trend consisting of a quadratic and a linear item in Z makes KT an outstanding method among all of the estimators examined in this thesis. If a population can not be divided into sub zones due to the lack of samples or of geological confining structures, KT is recommended however weak the correlation between the fitted regression curve and the data points is because the assumption of the existence of a trend in Z is supported by the qualitative interpretation of the lithology. The E-type estimate is usually different from the OK estimate which does not qualify the estimate for uncertainty. To generate a meaningful posterior distribution, however, more decisions have to be made, e.g., what thresholds to use and what extrapolation models to apply. The workload would be heavy if the median IK method is not appropriate. The statistics of IK are a bit better than those of OK with smaller MAE and MSE of the error distribution and smaller skewness and excess kurtosis of the population distribution, although IK tends to have smaller variances in both distributions, but the difference is not paramount considering the workload of IK aforementioned. LK seems to be the worst geostatistical method as it expects the log-transformed data to be normal. Direct non-linear algorithms are falling out of use nowadays owing to their sensitivity to back-transform. The kriging methods involved in the two simulation techniques do not make use of back-transform. The statistics of IDW are similar to those of OK since the sample dataset is not clustered as the samples are generally evenly distributed along the hole paths. The drawback of deterministic methods like IDW is that no probability is associated with an estimated value. 76 6.4 Estimation or Simulation Interpolation algorithms used in estimation tend to smooth out local details of the spatial variation of the attribute. Typically, small values are overestimated whereas large values are underestimated. Such conditional bias is a serious shortcoming when trying to detect patterns of extreme attribute values, such as zones of high RQD values. Another drawback of estimation is that the smoothing is not uniform as it depends on the local data configuration: smoothing is minimal close to the data locations and increases as the location being estimated gets farther away from data locations. A map of kriging estimates appears more variable in densely sampled areas than in sparsely sampled areas. Thus the kriging map may display artifact structures. Smooth interpolated maps should not be used for applications sensitive to the presence of extreme values and their patterns of continuity. Stochastic simulation generates a map or a realization of z-values which reproduces the statistics deemed most consequential for the problem in hand. Typical requirements for simulated maps are: • Data values are honored at their locations. The realization is then said to be conditional. • The histogram of the simulated values reproduces closely the declustered sample histogram. • The covariance model or the set of indicator covariance models for various thresholds are reproduced. The variance of kriged estimates is normally smaller than the sample variance. The semivariogram of kriged estimates has a much smaller relative nugget effect than the semivariogram model, which reflects the underestimation of the short-range variability of z- values. Since kriging methods focus more on local accuracy while simulation techniques attach more importance to the reproduction of spatial structures, the estimation method which generates a realization that has a transferred value closer to the mean of the realization distribution generated by a simulator ought to be preferred. 77 6.5 Recommendations The GSLIB developers encourage the customization of the source code to accommodate any special needs. In this study, the samples used in the estimation of each cell should have been prohibited from coming from the same borehole, that is, a parameter should have been incorporated into the implementation of the algorithm to cast a limit on the maximum number of samples that could come from the same borehole. However, due to the fact that all of the FORTRAN compilers are commercial, the free source code available for customization could not be compiled into executables. It would be worthwhile and recommended to examine what the effects would be if such a constraint is introduced into the system. The dimension of the block model is relatively large for the boreholes available. There are cells, especially those located in the corners of the block model, which are far away from the nearest borehole. Such a large distance, as mentioned above, will increase the smoothing effect. To avoid this, more drilling is recommended not only for this concern but also to define the trend more precisely. There are two reasons to use the sisim algorithm as the primary simulation tool. Firstly, checking the multiGaussian assumption is difficult and dependent on whether the sampling is sufficient; secondly, the Gaussian model does not allow for any significant spatial correlation of extremely large or small values, a property known as destructuration effect and associated to the maximum entropy property of the Gaussian random function model. Low entropy patterns such as connected patches of extreme values often correspond to features that are worth specific attention. In this case, both the samples and the interpolated images indicate that there exist connected clusters of large RQD values which render sgsim less appropriate than sisim. If safety is the major concern, an estimator returning the minimum value within the neighborhood or a lower quantile value might be the best. 78 7 Bibliography Anderson, T.W. 2003. An Introduction to Multivariate Statistical Analysis, 3 rd edn, John Wiley & Sons, New York. Barnes, R.J. 1991. The variogram sill and the sample variance. Mathematical Geology, 23, 673-678. Bieniawski, Z.T. 1989. Engineering rock mass classifications: a complete manual for engineers and geologists in mining, civil, and petroleum engineering. John Wiley & Sons, New York Bilgin, N. 1989. Applied Rock Cutting Mechanics for Civil and Mining Engineers. Birsen Publication, Istanbul. Birkhoff, G. and Rota, G. 1978. Ordinary differential equations. John Wiley & Sons, New York Bourgault, G. and Journel, A.G. 1994. Gaussian or indicator-based simulation? Which variograms is more relevant? In International Association for Mathematical Geology Annual Conference, p. 32-37. Bratley, O.B.L. Fox and Schrage, L.E. 1987. A Guide to Simulation. Springer-Verlag, New York. Chatterjee, S. and Price, B. 1977. Regression analysis by Example. Wiley, New York. Christensen, R. 1990. The equivalence of predictions from universal kriging and intrinsic random function kriging. Mathematical Geology, 22, 655-664. Clark, I. 1979. Practical Geostatistics. Applied Science Publishers, London. Coward, S., Vann, J., Dunham, S., and Stewart, M. 2009. The Primary-Response Framework for Geometallurgical Variables. Seventh International Mining Geology Conference. Cressie, N. and Zimmerman, D.L. 1992. On the stability of the geostatistical method. Mathematical Geology, 24, 45-59. David, M. 1977. Geostatistical Ore Reserve Estimation. Amsterdam, Elsevier. David, M. 1988. Handbook of Applied Advanced Geostatistical Ore Reserve Estimation. Amsterdam, Elsevier. Davis, J.C. 1986. Statistics and Data Analysis in Geology. 2 nd edn, John Wiley & Sons, New York. 79 Deere, D.U. and Deere, D.W. 1988. The rock quality designation (RQD) index in practice. In Rock classification systems for engineering purposes, (ed. L. Kirkaldie), ASTM Special Publication 984, 91-101. Philadelphia: Am. Soc. Test. Mat. Deere, D.U. 1989. Rock quality designation (RQD) after 20 years. U.S. Army Corps Engrs Contract Report GL-89-1. Vicksburg, MS: Waterways Experimental Station. Deere, D.U., Hendron, A.J., Patton, F.D. and Cording E.J. 1967. Design of surface and near surface construction in rock. In Failure and breakage of rock, proc. 8 th U.S. symp. rock mech., (ed. C. Fairhurst), 237-302. New York: Soc. Min. Engrs, Am. Inst. Min. Metall. Petrolm Engrs. Delfiner, P. 1976. Linear estimation of nonstationary spatial phenomena, in Guarascio, M., et al. (Eds), Advanced Geostatistics in the Mining Industry, Reidel, Dordrecht, 49-68. Deutsch, C.V. 1994a. Algorithmically-defined Random function Models. In R. Dimitrakopoulos, editor, Geostatistics for the Next Century, Kluwer, Dordrecht. Deutsch, C.V. and Journel, A.G. 1998. GSLIB: Geostatistical Software Library and User’s Guide. 2 nd edn, Oxford University Press, New York. Hustrulid, W.A. and Kuchta, M. 2006. Open Pit Mine Planning and Design. 2 nd edn, Taylor Francis, New York. Gardner, W.A. 1986. Introduction to Random Processes. Macmillan, New York. Gomez-Hernandez, J. and Cassiraga, E. 1994. Theory and practice of sequential simulation. In M. Armstrong and P.A. Dowd, editors, Geostatistical Simulation, p. 111-124. Kluwer Academic Publishers, Dordrecht. Goovaerts, P. 1997. Geostatistics for Natural Resources Evaluation, Oxford University Press, New York. Guibal, D. and Remacre, A. 1984. Local estimation of the recoverable reserves: comparing various methods with the reality on a porphyry copper deposit, Geostatistics for Natural Resource Characterization, D. Reidel Publishing Company. Isaaks, E.H. and Srivastava, R.M. 1989. An Introduction to Applied Geostatistics. Oxford University Press, New York. Isaaks, E.H. and Srivastava, R.M. 1987. Spatial continuity measures for probabilistic and deterministic geostatistics, MGUS conference. Journel, A.G. 1980. The lognormal approach to predicting local distributions of selective mining unit grades. Mathematical Geology, 12, 285-303. 80 Journel, A.G. 1983. Nonparametric estimation of spatial distribution. Mathematical Geology, 15, 445-468. Journel, A.G. and Huijbregts, C.J. 1978. Mining Geostatistics. Academic Press, London. Journel A.G. and Issaks, E.H. 1984. Conditional indicator simulation: application to a Saskatchewan uranium deposit. Mathematical Geology, 22(8), 1011-1025. Journel, A.G. and Rossi, M.E. 1989. When do we need a trend model in kriging? Mathematical Geology, 21(7), 715-739. Koch, G. and Link, R. 1986. Statistical Analysis of Geological Data. 2 nd edn, Wiley, New York. Mardia, K.V. and Marshall, R.L. 1984. Maximum likelihood estimation of models for residual covariance in spatial regression. Biometrika, 71, 135-146. Matheron, G. 1976. Transfer functions and their estimations. Advances Geostatistics in the Mining Industry, D. Reidel Publishing Company. Merritt, A.H. 1972. Geologic prediction for underground excavations. Proc. North American. rapid excav. tunneling conf., Chicago, (eds K.S. Lane and L.A. Garfield) 1, 115-132. New York: Soc. Min. Engrs, Am. Inst. Min. Metall. Petrolm Engrs. Olkin, I., Gleser, L. and Derman, C. 1980. Probability Models and Applications. Macmillan, New York. Pollard, H. 1920. The Stieltjes Integral and its Generalizations. Quarterly Journal of Pure and Applied Mathematics 19. Rendu, J.M. 1979. Normal and lognormal estimation. Mathematical Geology, 11, 407-422. Rendu, J.M. 1984. Geostatistical modeling and geological controls; Trans. Inst. Min. Metall, Sect. B., v. 93, p. B166-B172. Rendu, J.M. 1986. How the geologist can prevent a geostatistical study from running out of control: some suggestions; in Ranta, D.E., ed., Applied mining geology: ore reserve estimation; society of Mining Engineers of AIME, p. 11-17. Ripley, B.D. 1981. Spatial Statistics. Wiley, New York. Rivoirard, J. 1994. Introduction to disjunctive kriging and non linear geostatistics. Oxford University Press. Sen, Z. 2009. Spatial Modeling Principles in Earth Sciences. Springer-Verlag, New York. 81 Sichel, H.S. 1952. New Methods in the Statistical Evaluation of Mine Sampling Data. Trans of the London Inst of Mining and Metallurgy, 61, 261-288. Sinclair, A.J. and Giroux, G.H. 1984. Geological controls of semi-variograms in precious metal deposits; in Verley, G., et al., Geostatistics for Natural Resources Characterization, Part 2, p. 965-978, D. Reidel Pub. Co., Dordrecht, Netherlands. Starks, T. and Fang, J. 1982a. The effect of drift on the experimental semi-variogram: J. Int. Assoc. Mathematical Geology, 14, 57-64. Stein, M.L. 1989. The loss of efficiency in kriging prediction caused by misspecification of the covariance structure, in Armstrong, M. et al. (Eds.), Geostatistics vol. 1: Kluwer, 273-282. Sullivan, J. 1984. Conditional recovery estimation through Probability Kriging, in Verly, G. et al. (Eds.), Geostatistics for Natural Resources Characterization: Reidel, Dordrecht, 365-384. Verly, G. 1983. The multiGaussian approach and its application to the estimation of local reserves. Mathematical Geology, 15, 259-286. Villaescusa, E. and Potvin, Y. 2004. Ground Support in Mining & Underground Construction. Taylor and Francis, London. Wackernagel, H. 1995. Multivariate Geostatistics: An introduction with applications. Springer-Verlag, Berlin, 256 p. Wardrop Engineering Inc. 2009. Technical Report on the Mactung Deposit – Yukon, Canada Weisberg, S. 1985. Applied Linear Regression. Wiley, New York. Zimmerman, D.L. and Cressie, N. 1991. Mean squared prediction error in the spatial linear model with estimated covariance parameters: Ann. Ins. Stat. Math.