**PAPERS OF GENERAL INTEREST**

**Accounting for a spatial trend in fine-scale ground-penetrating radar data: a comparative case study**

**Y. Dagasan ^{I}; O. Erten^{I}; E. Topalf^{II}**

^{I}Mining and Metallurgical Engineering Department, WA School of Mines, Curtin University, Kalgoorlie, WA, Australia

^{II}School of Mining and Geosciences, Nazarbayev University, Astana, Kazakhstan

]]>

**SYNOPSIS**

In geostatistics, one of the challenges is to account for the spatial trend that is evident in a data-set. Two well-known kriging algorithms, namely universal kriging (UK) and intrinsic random function of order *k *(IRF-k), are mainly used to deal with the trend apparent in the data-set. These two algorithms differ in the way they account for the trend and they both have different advantages and drawbacks. In this study, the performances of UK, IRF-k, and ordinary kriging (OK) methods are compared on densely sampled ground-penetrating radar (GPR) data acquired to assist in delineation of the ore and waste contact within a laterite-type bauxite deposit. The original GPR data was first pre-processed to generate prediction and validation datasets in order to compare the estimation performance of each kriging algorithm. The structural analysis required for each algorithm was carried out and the resulting variograms and generalized covariance models were verified through cross-validation. The variable representing the elevation of the ore unit base was then estimated at the unknown locations using the prediction data-set. The estimated values were compared against the validation data using mean absolute error (MAE) and mean squared error (MSE) criteria. The results show although IRF*-k *slightly outperformed OK and UK, all the algorithms produced satisfactory and similar results. MSE values obtained from the comparison with the validation data were 0.1267, 0.1322, and 0.1349 for IRF-k, OK, and UK algorithms respectively. The similarity in the results generated by these algorithms is explained by the existence of a large data-set and the chosen neighbourhood parameters for the kriging technique.

**Keywords**: ground-penetrating radar, geostatistics, nonstationarity, universal kriging, ordinary kriging, intrinsic random function of order *k.*

**Introduction**

In the geological modelling of the laterite-type deposits, the exploration boreholes may be sparsely spaced for two major reasons: (1) the attributed grades do not tend to vary significantly across the deposit; (2) the overall exploration costs must be minimized. However, since the geological contact between the ore and underlying waste unit fluctuates in a rather complex manner, one should not rely solely on the interpolation of the sparsely-spaced borehole data, as the estimates would not reflect the actual variability in the surface. Therefore, it is imperative that the sparsely-spaced borehole data be supplemented with appropriate geophysical information for more accurate resource estimations (Erten *et al., *2015, 2013). Ground-penetrating radar (GP R) has widely been used to acquire a complementary dense data-set to better delineate the interface between two geological units (bauxite/ferricrete) (Francke and Nobes 2000; Francke and Parkinson 2000; Watts 1997). Due to the tropical weathering and leaching mechanisms that generate laterite-type deposits, the geological interface measured through GPR appears to be well correlated with the easting (X) and northing (Y) coordinates of the surface, which indicates the presence of a spatial trend (McLennan, Ortiz, and Deutch, 2006; Leuangthong, Lyall, and Deutsch, 2002). Geostatistical techniques are the main tools which are used to estimate the elevation of the interface at non-sampled locations, and there are different approaches used to account for the trend in a data-set. In this paper, we estimate the elevation of the geological interface through ordinary kriging (OK), universal kriging (UK), intrinsic random function of order *k *(RF-k) methods and compare the estimation performances.

Geostatistical techniques are based on the theory of regionalized variables (RVs) and are used to estimate an attribute of interest at non-sampled locations (Goovaerts, 1997; Journel and Huijbregts, 1978). The idea behind the theory is that a RV z**(u**_{α}),{α+1.....n} is considered to be the realization of an order-two random function (RF) *Z*** (u) **and is assumed to have been generated according to a probability density function (Matheron, 1971; Olea, 1974). Due to this characteristic of the RV, there is a spatial correlation between the samples, which allows the prediction of

*Z**

**(u)**at each non-sampled location

**u.**The prediction of Z

***(u)**requires the covariance function C

**(h)**of Z

**(u)**to be known, and this statistical inference can practically be made from an available realization if the realization exhibits stationary characteristics. If there is a trend in the realization/data-set and the realization exhibits nonstationary characteristics, the change in the value of RV is no longer dependent on the lag distance

**h,**but is also dependent on the location of the RV. In the case of a nonstationarity, there are basically two geostatistical methods to account for the trend in the data-set: universal kriging and intrinsic random function of order

*k*(IRF-k). These two methods differ in the way they detect the trend and the type of structural function used to describe the spatial relationship.

In order to estimate the value of an attribute at a non-sampled location, the kriging algorithm requires the computation of a system of equations with a known variogram *γ***(h). **The problem with UK arises when the underlying variogram or covariance function is not known. In UK, the RF Z**(u) **is comprised of a trend m**(u) **and a residual R**(u) **component. The variogram of residuals can be used to calculate the underlying variogram, but this requires the determination of the trend component *m**** (u). **Calculation of m

***(u),**on the other hand, requires the variogram to be known. One solution to this circular problem is refining both the variogram and the trend estimates iteratively, as mentioned in Neuman and Jacobson (1984). However, it has some associated drawbacks, such as underestimation of the underlying variogram and extreme difficulties in determining the degree of trend or the underlying variogram from the residuals (Armstrong, 1984; Cressie, 1993). Another solution to these problems exists if there is a particular direction or subzone where there is a sub-stationary zone in a nonstationary data-set. Journel and Rossi (1989) proposed that experimental variograms inferred in these directions or subzones can be used to estimate the trend optimally. Chilès (1976) gives an example of such practice assuming an isotropic variogram model computed from the sample variogram in stationary directions or subzones.

IRF-k was introduced by Matheron (1973) due to the practical difficulties in the application of the UK approach. It basically decomposes the trend and covariance structure through increments of a sufficient order to filter out the trend and to achieve stationarity (Chiles and Delfiner, 2012). The kriging system of the IRF*-k *method is identical to the UK method except that the variogram employed in UK is replaced by the generalized covariance (GC) in the IRF-k method. Contrary to UK, where the trend is required to be estimated beforehand as a linear combination of known, linearly independent functions to obtain the stationary residuals, nonstationarity in IRF*-k *is accounted for through the calculation of the GC. The advantage of using GC is that it has wider class of admissible functions compared to the variogram and the automatic detection of the parameters of GC makes the application easier (Delhomme, 1978). On the other hand, the use of GC creates some hurdles in practical modelling since the method requires identification of the order *k *and it is difficult to interpret the GC (Cressie, 1993). Both UK and IRF*-k *methods have advantages and disadvantages and as a result, the choice between the two methods is based heavily on the practical difficulties of fitting functions (Buttafuoco and Castrignano, 2005). Christensen (1990) states that the IRF*-k *and UK are identical provided that the GC is identified correctly.

*-k,*and OK have been compared in several publications. Journel and Rossi (1989) compared OK and UK in a case study in which they regionalized seam thickness and coal quality variables. The results indicated that the OK and UK methods gave similar results. It was also concluded that any kriging algorithm with moving data windows is equal to considering a nonstationary random function model with a mean re-estimated at each new location. Similarly, Zimmerman

*et al.*(1999) compared the performances of four interpolation algorithms, two of which were OK and UK. They stated that although UK was expected to outperform OK in situations where trends exist, OK performed slightly better than UK. Odeh, McBratney, and Chittleborough (1994) modelled soil variables by using different interpolation methods, and stated that OK was the most inferior of all the methods implemented, including UK. Odeh, McBratney, and Slater (1997) compared the performance of several prediction models, including OK and IRF-k, in prediction of soil parameters. The IRF

*-k*method performed slightly better than OK. It is therefore apparent that the performance ranking of UK, IRF

*-k,*and OK algorithms varies between different investigations.

The objective of this study is to predict the variability in elevation of the base of the ore unit and compare the performance of different kriging estimators using densely sampled GPR data. The contribution of this paper is the implementation of, and comparison of, the performance of, different kriging algorithms in the case of a bauxite laterite deposit. This was mainly done in order to ascertain which kriging algorithm is more suitable for bauxite base elevation data with a spatial trend. Since GPR data represents the elevation of an interface surface, nonstationarity was intuitively expected. This was also confirmed by the omnidirectional variograms computed in the initial data analysis. The omnidirectional variograms revealed the existence of nonstationarity in the prediction data-set, requiring handling of the trend by UK and IRF-k methods. Therefore, these methods were employed to account for the spatial trend in the data-set along with the OK algorithm, which considers the spatial trend to be constant. In order to evaluate the performances of each kriging algorithm, GPR data representing the elevation was resampled to form validation and prediction data-sets. The validation data-set was used to assess the performance of the UK, IRF-k, and OK estimators.

**General review of methods**

*Geostatistical theory*

The RF can be represented by the following model:

where m**(u) **represents the trend component and *R*** (u) **represents the random part having a covariance function in two-dimensional space

*u**=(x,y).*The covariance C

**(h)**is defined as the following:

where **h **is the lag distance. Considering the second-orderstationarity assumption, *R***(u) **is assumed to have a zero mean value and the trend *m*** (u) **is assumed to be constant (Oliver and Webster, 2015). Hence, C

**(h)**is equal to:

In situations where the mean is not constant, the covariance cannot exist. In these cases, the assumption of stationarity is weakened to the one called *intrinsic stationarity *(Matheron, 1963), where the expected differences E[Z**(u+h**)-Z**(u**)]are equal to zero and the covariance of the random part, which is used to measure the spatial relations of the residuals, is replaced by the variance of the differences:

where 7**(h) **represents the variogram at lag **h. **The following relationship between the variogram 7**(h) **and covariance C**(h) **function exists for a RF which is second-order stationary:

where C**(0) **represents the variance σ^{2}of the RF.

*Ordinary kriging (OK)*

In OK, the RV is estimated at an unsampled location **u **as a linear combination of available neighbouring data {z**(u**„),α=1.....n} (Goovaerts, 2000).

As with all kriging algorithms, the objective of OK is to minimize the estimation variance as well as ensure the unbiasedness of the estimator. OK weights *Xß*** (u) **are determined to satisfy this unbiasedness and minimized variance goals. The following OK system is used to obtain the kriging weights:

where μ**(u) **represents the Lagrange parameter used to minimize the variance. The only information needed for the OK system is the variogram value corresponding to every lag distance **h. **The kriging variance of OK is calculated by:

*Universal kriging (UK)*

In UK, the trend component is modelled as a smoothly varying deterministic function of **u **and is expressed as:

where m**(u) **represents the local mean, *at, *the unknown coefficients of the trend function, and *f***(u) **represent the known functions of the coordinates **u **and are called *trial *or *basejunctions *(Kitanidis, 1997; Rossi and Deutsch, 2014). The residual component R**(u) **is modelled as a second-order stationary RF with a zero-mean, *E*[Z**(u**)]=0. Combining the trend and the residual, the RV is represented by the following equation:

The trend can be modelled as a low-order polynomial function of the spatial coordinates **u. **Previous studies have indicated that increasing the number of order of trend functions does not lead to better modelling of the trend (Journel and Rossi, 1989).

To satisfy the unbiasedness condition E[Z***(u**)-Z**(u**)]=0, Equation [9] can be rewritten as:

where the set of *L *constraints is termed the *universality *or *unbiasedness conditions. *The kriging system satisfying these requirements is then defined by:

where *λ*_{β} are the UK weights and *μι*are the Lagrange parameters used to determine the coefficients *λ*_{β}minimizing the error variance.

The UK estimate variance is defined as:

The inference of the residual variogram y** _{R}(h) **is, however, not straightforward, as the only available data is in fact

*Z*values, not

*R*values. The experimental variogram of

*Z*values is defined by:

For UK, the underlying (trend-free) variogram is assumed to be known (Armstrong, 1984). To predict this variogram, the form of the trend should be known. However, to estimate the form of the trend, the variogram must be known. Various approaches have been developed to resolve this circular problem. The most common approach is to select a direction or a subzone in which the trend m**(u) **can be negligible and calculating the experimental covariance along these selected directions or within this subzone (Atkinson and Lloyd. 2007; Chihi *et al., *2000; Journel and Rossi, 1989; Myers, 1989).

*Intrinsic random function of order k (IRF-k)*

*k*was developed by Matheron (1973) and later refined by Delfiner (1976). It is an alternative tool used in nonstationary data to remove the trend by filtering out the low-order polynomials and is commonly used when it is difficult to infer the underlying

variogram of the variable of interest (Wackernagel, 2002). In the IRF-k method, the deterministic functions *f _{1}*

**(u)**used to represent the trend model

*m*

**in Equation [9] are restricted to those which are only translation-invariant and pairwise orthogonal (Wackernagel, 2002).**

*(u)*Considering a set of weights *λ*_{α}applied to particular points **u****„****, **a discrete measure **u****„** is defined as:

where *δ** _{ua}*represents the Dirac measure at point

**u**

**„**

**.**Any linear combinations of the weights with RV at locations

**u**are defined as:

_{a }The expression shown in Equation [16] is called the allowable linear combination of order *k(ALC-k) *in Equation [17]

holds true for all monomials of order <k. Considering a nonstationary Z**(u), **if the expression given in Equation [16] is second-order stationary regardless of any translation **h **and whatever the ALC-k λ, the RV Z**(u) **is called as IRF-k.

The variogram in IRF-k is replaced by a new function called the GC function (Chiles and Delfiner, 2012). GC is denoted by *K*** (h) **and is used to describe the correlation structure of the random part R

**(u).**The GC function of an IRF-

*k*is a symmetric function, K

**(h**)=K

**(-h),**and satisfies the condition:

An example to GC function is known the as *polynomial generalized covarianceJunction. *The equation of the polynomial generalized covariance function is given as:

The conditions on *b** _{p}*are satisfied if bp<0 Vp. The intrinsic kriging system minimizing the variance is expressed in terms of the GC and shown in the following:

The estimation variance is defined by:

**Case study**

*Geological setting*

The bauxite deposit in the mine area is thought to have formed from *in-situ *chemical weathering of kaolinite, quartz, and iron oxide minerals (Loughnan and Bayliss, 1961). The occurrence of the alumina-rich was horizon controlled by climate, vegetation cover, chemical conditions, bedrock composition and texture, groundwater circulation, relief, time, and tectonic conditions (Gow and Lozej, 1993). There is a regolith zone in the mine area comprising, from top to bottom, post-weathering sediments (red soil), bauxitic cement, pisolitic bauxite, nodular ferricrete, the kaolinite zone, and the saprolitic zone (Bardossy and Aleva 1990).

*The GPR survey*

An electromagnetic (EM) wave that travels through shallow ground shows different responses to subsurface structures with varying electromagnetic properties such as dielectric permittivity, conductivity, and electromagnetic permeability. GPR utilizes the dielectric permittivity contrast that exists between the geological structures (ASTM D 6432-99 2005). The EM wave emitted into the ground from the transmitting antenna at the surface is reflected back when there is a difference in the electrical properties of the subsurface structures. This reflected wave is then received by a receiving antenna and recorded as a function of time (Davis and Annan, 1989). Being a function of depth, antenna spacing, and average radar-wave velocity, the time taken for the wave to travel to the interface and back up to the surface is called the *two-way travel time. *Knowing the radar wave velocities, this two-way travel time is then converted into depth.

*Selection of GPR method to map the lateral variability at the mine site*

The main aim of the survey was to map the lateral variability at the bauxite/ferricrete interface. As the area to be surveyed was large (360 χ 800 m), the method chosen needed to be easy to implement as well as provide fast data acquisition so that the surveying results could be checked immediately. Other considerations in selecting a suitable geophysical method were spatial resolution capability, cost-effectiveness, and data processing requirements (Erten, 2012).

In order to be sure about the applicability of the method selected, the petrophysical properties established in previous work on samples collected from another location within the mine site were considered. The laboratory results revealed that the conductivities of the bauxite and ferricrete are rather low and there is a dielectric permittivity contrast between the bauxite and the ferricrete. This indicated that the emitted waves would reflect from the bauxite interface, favouring the use of the GPR method.

*The data*

The mine area chosen in this case study is approximately 360 χ 800 m in size. The data-set comprises the GPR pick-points acquired from the surface of the mine area using a radar device. The specifications of the radar device are given by Francke and Utsi (2009). The total areal coverage of the GPR survey at the mine area was 142 300 m2 with GPR profiles 17 940 m in length, which provided 64 670 GPR pick-points distributed along the GPR profiles at 0.25 m interval (Figure 1). The GPR profiles were arranged in a square grid with a line spacing of about 15 χ 15 m (Erten, 2012).

]]>

Due to the nature of GPR data acquisition, it is expected that there might be multiple points that have the same coordinates. Therefore, the data-set was first pre-processed to mask the potential duplicates in order to avoid any kriging matrix instability in the estimation process. This was carried out by masking all the sample points that were within 0.25 m of the other data-points. The number of raw GPR pick-points, after processing, reduced to 30 630. In order to compare the performance of each kriging method, the data was split into two parts: prediction and validation data-sets. The prediction data-set was generated by re-sampling the GPR pick-points randomly on a regular 15 χ 15 m grid. This process yielded 735 pick-points as a prediction/training data-set to be used in the estimation. The remaining 29 895 pick-points were kept for the validation. The flow chart of the data processing and the methodology are presented in Figure 2.

*Exploratory data analysis*

The target variable to be regionalized in this study is the elevation of the bauxite/ferricrete interface. This variable was obtained through the GPR survey, the depth to the interface being calculated from the two-way travel time. This depth was then subtracted from the ground surface elevation to obtain the elevation of the bauxite/ferricrete interface. The unit of the variable is metres (m).

The histograms of both prediction and validation data corresponding to the elevations of the base of the ore unit are shown in Figure 3.

As it can be seen in Figure 3, the histograms of the prediction and validation data are closely similar. This resemblance is also revealed in the descriptive statistics of these two data-sets given in Table I.

]]>

The calculated coefficients of skewness of the prediction/training and validation data-sets are -0.70 and -0.64, respectively.

Scattergrams plotted against X and Y coordinates in Figure 4 reveals that the elevation variable decreases towards the east and increases towards the north, which suggests a possible trend dipping in southeasterly direction.

*Structural analysis*

The raw omnidirectional variograms of the ore unit base elevation were computed from the prediction data with 50 lags, having a lag distance of 15 m and a lag tolerance of 7.5 m. This lag distance was the average distance between the samples in the prediction data-set. The experimental omnidirectional variogram is shown in Figure 5.

]]> The apparent increase in the variogram parallel to the increasing lag distance also confirms the existence of a trend or nonstationarity of the elevation variable. This conclusion is supported by the scattergrams shown in Figure 4.

A variogram map was computed to detect the maximum and minimum spatial continuity directions as well as any possible anisotropies. The variogram map and experimental variograms in these directions are shown in Figure 6.

The maximum spatial continuity in the variogram map is seen in the N20° direction as the variogram values do not change significantly along this direction. On the other hand, perpendicular to this direction, which is N110°, the elevation varies more rapidly, indicating the minimum spatial continuity direction. In the N70° direction, there are relatively high variogram values. As seen in Figure 6b, the variogram values in the N70° direction are higher only for the lag values greater than 180. Since the GPR data is densely sampled, only samples that are around 80-100 m apart from the estimation grid node are considered for the estimation and are of greatest importance. This fact confirms the selection of the N110° direction as the minimum spatial continuity direction, since the experimental variogram values in the N110° direction are higher in the 0-100 lag distance range.

*Structural modelling for UK*

UK was the first approach used in this study to account for the nonstationarity seen in the data. UK requires a prior determination of the *L *trend functions *F**i*** (u) **and the covariance Cr

**(h)**of the residual component R

**(u)**inferred from the residual variogram yR

**(h)**(Goovaerts, 1997). The the

*L*trend functions were identified by testing low-order polynomial functions, k<2, by locally fitting the polynomial functions through the ordinary least-squares method. The trend component identification of the elevation variable is summarized in Table II.

In order to identify the order *k *of the trend component, the degree of trend, a number of varying orders of polynomials were selected and fitted to the data-points. In this process, errors obtained from the fitting of different orders of *k *polynomials are calculated at each point and ranked by order of absolute magnitude. The ranks obtained from each point are then averaged and the order *k *having the lowest average rank selected as true *k. *In this case study, polynomials of order 0, 1, and 2 were tried. As is seen from Table II, the smallest mean rank value was 1.850, and it was yielded by the first-order polynomial, k=1. Therefore, the first-order trend function, which is comprised of 1, *x,y *monomials, was selected as the best-fit trend. It would, therefore, be expected that the OK, UK, and IRF*-k *estimates would be similar due to the linear model fitted to the trend.

*n*times,

*n*being a user-defined maximum iteration number. Variogram fitting for UK was carried out by using ISATIS software. The mathematical expression of the model fitted to the underlying variogram model is as follows:

The mathematical model fitted to the underlying variogram is given in Figure 7.

*Structural modelling for IRF-k*

The second approach used to account for the nonstationarity was the IRF-k method. Nonstationary modelling with the IRF-k method involves three steps: (1) determination of order k, (2) inference of the optimal GC model, and (3) kriging the elevation variable based on the inferred GC. The GC model is readjusted (if necessary) based on the comparison of the obtained error with the theoretical standard deviations (Delfiner, 1976). Trend analysis has already been carried out for the UK case and the results were summarized in Table II. It was found that the best fit trend for the variable was the order 1 trend (k=1).

In order to infer the optimal GC, several arbitrary models of GCs were proposed and tested. The optimality of a GC model was determined on the basis of the ratio called *'JackniJe' *which is basically the number indicating the ratio of the theoretical variance to the experimental one (Farkhutdinov *et al., *2016). The calculated Jacknife numbers were ranked in ascending order and the GC model corresponding to the Jacknife number closest to unity was selected as the optimal GC. The tested GC models and their scores can be seen in Table III.

It can be seen in Table III that the order 1 GC function yielded the Jacknife score closest to unity. Hence, order 1 type of the GC function was selected with a sill value of 3.148 and a range value of 170 m. Any possible anisotropy that may exist in the data-set is accounted for by the polynomial function which filters out the trend (Delfiner, 1976).

*Structural modelling for OK*

In the OK case, the residual variogram yR **(h) **was inferred by calculating the experimental variogram of Z**(u) **along the direction in which the trend m**(u) **was deemed negligible. The deemed trend-free direction was detected as N20° from the variogram map shown in Figure 6, and the experimental variogram and the model fitted are shown in Figure 8.

]]>

However, since the GPR data is densely sampled, only the samples that are 80-100 m away from any estimation grid are considered and are of greater importance. Therefore, the variogram model reaching the sill value at range value of 256 is considered as robust.

Two spherical models without a nugget variance provided the best fit to the experimental variogram computed by using 15 m as the lag distance and 0.5 **h **as the lag tolerance. The mathematical representation of the fitted model is shown in the following:

*Cross-validation*

The cross-validation technique was used to assess the accuracy of the variogram models fitted to the experimental variograms. It utilizes diagnostic statistics and the accuracy of the prediction is evaluated through various tools (Webster and Oliver, 2001). The criteria used to estimate the accuracy in this study are the mean absolute error (MAE) and mean squared deviation ratio (MSDR). The MAE should ideally be zero, which satisfies the unbiasedness condition. The MSDR is basically the ratio of the computed squared errors to the kriging variances, and the closer the MSDR is to unity, the better the model for kriging (Oliver 2010). These two criteria are calculated as follows:

where *N *is the number of data values (which is 735 for this study), Z**(u****„****) **is the true value, Z***(u****„****) **is the predicted value, and o2**(u) **is the kriging variance. The results of the cross-validation technique that was implemented to assess the accuracy of the variograms used for UK, IRF*-k, *and OK are given in Table IV.

]]>

The results of the diagnostic statistics show that, although the variogram modelled for OK has the minimum MAE and the closest MSDR value to unity, all of the computed MAE values are small, indicating that there is no significant difference among the error statistics. Based on these results, all of the modelled variograms are considered to be appropriate for spatial prediction.

*Prediction ofore unit base elevation*

The first step of the prediction of the elevation of the ore unit base is to define an estimation grid that is capable of covering the whole area containing the data-points as well as minimizing the extrapolation. In this case study, a two-dimensional grid having the following properties was created: the origin is *X _{0}*= -1872.00 m, y

_{0}= 10181.00 m, the dimensions of the mesh are

*dx*= 5.00 m,

*dy*= 5.00 m, the number of meshes is 75 along the

*x*direction and 163 along the

*y*direction, resulting in 12 225 meshes in total. This was then followed by constructing a polygon delineating the boundaries of the data, which defines the resource estimation area. The number of meshes within the polygon boundary is 7 089.

*Neighbourhood parameters*

A moving type of neighbourhood was used in all the aforementioned algorithms to estimate the ore unit base elevations. In order to make the data around the estimation grid evenly distributed, eight angular sectors around the estimation grid were defined with a minimum of four samples in each sector. Delfiner (1976) states that the number of angular sectors should at least be twice the number of unbiasedness conditions, and considering the linear trend, there should at least be six angular sectors around the data. Hence, the number of angular sectors selected satisfies the given rule of thumb. Based on the parameters selected, a minimum of 32 samples around the estimation grid were used in each search neighbourhood. The radius of the search window circle was determined to be 170 m. However, since GPR pick-points were densely sampled, 32 samples around any estimation grid node do not fall further than 80 m radius. In another words, points that are 80 m away from any estimation grid node are not taken into account for estimations through kriging techniques due to this neighbourhood selection. This selected neighbourhood was used for all the kriging methods.

*Predictions with UK, IRF-k, and OK techniques*

UK was used to estimate the ore unit base elevation variable by making use of the inferred underlying variogram of the residuals and the trend function, which describe the spatial relationships between the sample data. The spatial maps of the ore unit base elevation obtained from the UK estimation are given in Figure 9.

Contrary to UK, IRF-k does not require the trend and the underlying variogram of the residuals to be determined beforehand, since it has its own automatic structure identification algorithm allowing it to pick up the best set of parameters within a preselected set of models (Chiles and Delfiner, 2012). Hence, parameters describing the spatial relationship for IRF-k algorithm in this study were automatically detected. These parameters are the trend and the chosen optimal GC model. The spatial maps of the estimates from IRF*-k *method are given in Figure 10.

**Results and discussion**

The performance of the kriging estimators was tested by comparing the validation data with the kriged data. This was done by comparing the mean squared error (MSE) and MAE values. In order to compute these errors, kriged data, which was collocated with the validation data, was copied from each of the UK, IRF*-k, *and OK maps and used to calculate the error associated with each estimation algorithm. The results of MSE and MAE are shown in Table V.

The differences between the errors calculated from the different kriging algorithms are similar and the errors are not significant, considering the mean 15.01 m of the prediction data-set.

Scattergrams generated by plotting the kriged values obtained from the three kriged algorithms against the validation data yielded almost identical results. The coefficients of correlations were 0.990, 0.990, and 0.991 for UK, OK, and IRF*-k, *respectively.

*-k*slightly outperformed the other predictor algorithms, the estimation errors were not significant enough to conclude that the more sophisticated IRF

*-k*algorithm outperformed OK in this particular case study. The similarity in the results yielded by these techniques is thought to be due to factors such as the densely sampled GPR elevation variable and the selected neighbourhood parameters.

The elevation variables were densely sampled and a maximum of 32 samples were used in the estimation of the variable at an unknown grid node due to the neighbourhood definition. Since 32 samples were almost always within a 80100 m radius, the trend was probably not apparent within such a radius. Therefore, the effect of the trend was not experienced to any significant degree, since all the samples used in the estimation were from the immediate vicinity of the estimation grid.

**Conclusions**

The focus of this case study was to implement and compare the performances of different geostatistical estimators in the case of a trend apparent in a densely sampled GPR data-set. The UK and IRF*-k *methods were implemented to account for the trend seen in the data-set and OK was implemented by considering the spatial trend to be constant.

The performances were assessed by comparing the kriged values with the preselected validation data for each kriging algorithm. The results of the comparisons with the validation data have shown that IRF*-k *outperformed the other algorithms considering MAE and MSE criteria. However, the differences between the results were not sufficiently significant for one kriging algorithm to stand out among the others. For example, the MSE values obtained from the comparison with the validation data were 0.1267, 0.1322, and 0.1349 for the IRF-k, OK, and UK algorithms, respectively. Similarly, the plotted scattergrams demonstrated a similar outcome; the coefficient of correlations obtained from plotting the kriged values against the validation data were 0.990, 0.990, and 0.991 for UK, OK, and IRF-k, respectively. This similarity was mainly due to the large dataset and neighbourhood parameters chosen.

**References**

Armstrong, M. 1984. Problems with universal kriging. *Journal of the International Association/or Mathematical Geology, *vol. 16. pp. 101-108. doi:10.1007/Bf01036241 [ Links ]

ASTM D 6432-99. 2005. Standard guide for using the surface ground penetrating radar method for subsurface investigation. ASTM International, West Conshohocken, PA. [ Links ]

Atkinson, P.M. and Lloyd, C.D. 2007. Non-stationary variogram models for geostatistical sampling optimisation: An emprical investigation using elevation data. *Computers and Geosciences, *vol. 33. pp. 1285-1300. [ Links ]

Bardossy, G. and Aleva, G.J.J. 1990. Lateritic Bauxites. Elsevier, Amsterdam. [ Links ]

Buttafuoco, G. and Castrignano, A. 2005. Study of the spatio-temporal variation of soil moisture under forest using intrinsic random functions of order k. *Geoderma, *vol. 128. pp. 208-220. doi:10.1016/j.geoderma.2005.04.004 [ Links ]

Chihi, H., Galli, Α., Ravenne, C., Tesson, M., and de Marsily, G. 2000. Estimating the depth of stratigraphic units from marine seismic profiles using nonstationary geostatistics. *Natural Resources Research, *vol. 9. pp. 77-95. [ Links ]

Chiles, J-P. 1976. How to adapt kriging to non-classical problems: three case studies. *Advanced Geostatistics in the Mining Industry: Proceedings of the NATO Advanced Study Institute, *Istituto di Geologia Applicata of the University of Rome, Italy, 13-25 October 1975. Springer. pp. 69-89. [ Links ]

Chiles, J.P. and Delfiner, P. 2012. Geostatistics: Modeling Spatial Uncertainty). Wiley. [ Links ]

Christensen, R. 1990. The equivalence of predictions from universal kriging and intrinsic random-function kriging. *Mathematical Geology, *vol. 22. pp. 655-664. doi:10.1007/Bf00890514 [ Links ]

Cressie, N. 1993 Statistics for Spatial Data. Wiley-Interscience, New York. [ Links ]

Davis, J.L. and Annan, A.P. 1989. Ground-penetrating radar for high-resolution mapping of soil and rock stratigraphy. *Geological Prospecting, *vol. 37. pp. 531-551. [ Links ]

Delfiner, P. 1976. Linear estimation of non stationary spatial phenomena. *Advanced Geostatistics in the Mining Industry. *Springer. pp. 49-68. [ Links ]

Delhomme, J.P. 1978. Kriging in the hydrosciences. *Advances in Water Resources, *vol. 1. pp. 251-266. doi:10.1016/0309-1708(78)90039-8 [ Links ]

ERTEN, O. 2012. Profiling and mining control to mitigate dilution effect from SiO_{2} at the base of a bauxite deposit. PhD thesis, University of Queensland. [ Links ]

Erten, O., Kizil, M.S., Topal, E., and McAndrew, L. 2013. Spatial prediction of lateral variability of a laterite-type bauxite horizon using ancillary ground-penetrating radar data. *Natural Resources Research, *vol. 22. pp. 207-227. [ Links ]

Erten, O., McAndrew, L., Kizil, M.S., and Topal, E. 2015. Incorporating fine-scale ground-penetrating radar data into the mapping of lateral variability of a laterite-type bauxite horizon. *Mining Technology, *vol. 124. pp. 1-15. [ Links ]

Fallon, G.N., Fullagar, P.K., and Sheard, S.N. 1997. Application of geophysics in metalliferous mines. *Australian Journal of Earth Sciences, *vol. 44. pp. 391-409. [ Links ]

Farkhutdinov, Α., Goblet, P., de Fouquet, C., and Cherkasov, S. 2016. Α case study of the modeling of a hydrothermal reservoir: Khankala deposit of geothermal waters. *Geothermics, *vol. 59, Part Α. pp. 56-66. http://dx.doi.org/10.1016/j.geothermics.2015.10.005 [ Links ]

Francke, J. and Utsi, V. 2009. Advances in long-range GPR systems and their applications to mineral exploration, geotechnical and static correction problems. *First Break, *vol. 27. pp. 85-93. [ Links ]

Francke, J.C. and Nobes, D.C. 2000. Α prelimenary evaluation of GPR for nickel laterite exploration. *Proceedings of the Eighth International Conference on Ground Penetrating Radar, *Gold Coast, Queensland. Noon, D.A., Stickley, G.F., and Longstaff, D. (eds). SPIE, Bellingham, WA. pp. 7-12. [ Links ]

Francke, J.C. and Parkinson, G. 2000. The new role of geophysics in nickel laterite exploitation and development. *Proceedings of the Mining Millennium/PDAC Convention, *Toronto, Canada. Prospectores and Developers Association of Canada. [ Links ]

Goovaerts, P. 1997. Geostatistics for Natural Resource Evaluation. Oxford University Press, New York. [ Links ]

Goovaerts, P. 2000. Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. *Journal of Hydrology, *vol. 228. pp. 113-129. doi: 10.1016/S0022-1694(00)00144-X [ Links ]

Gow, N.N. and Lozej, G.P. 1993. Bauxite. *Geoscience Canada, *vol. 20. pp. 9-16. [ Links ]

ISATIS. 2016. Software manual, 9.0 edn. Avon Cedex, France. [ Links ]

Journel, A.G. and Huijbregts, C.J. 1978. Mining Geostatistics. Academic Press, London. [ Links ]

Journel, A.G. and Rossi, M.E. 1989. When do we need a trend model in kriging? *Mathematical Geology, *vol. 21. pp. 715-739. [ Links ]

Kitanidis, P.K. 1997. Introduction to Geostatistics: Applications to Hydrogeology. Cambridge University Press. [ Links ]

Leuangthong, O., Lyall, G., and Deutsch, C. 2002. Multivariate geostatistical simulation of a nickel laterite deposit. APCOM 2002: *Proceedings of the 30th International Symposium on the Application of Computers and Operations Research in the Mineral Industry, *University of Alaska. [ Links ]

Bandopadhyay, S. (ed.). Society for Mining, Metallurgy, and Exploration, Littleton, CO. pp. 261-273. [ Links ]

Loughnan, F.C. and Bayliss, P. 1961. The mineralogy of the bauxite deposits near Weipa, Queensland. *American Mineralogist, *vol. 46. pp. 209-217. [ Links ]

Matheron, G. 1963. Principles of geostatistics. *Economic Geology, *vol. 58. pp. 1246-1266. [ Links ]

Matheron, G. 1971. The theory of regionalised variables and its applications. Ecoledes Mines, Cahiers du Centre de Morphologie Mathematique, Paris. vol. 5. pp. 211. [ Links ]

Matheron, G. 1973. The intrinsic random functions and their applications *Advances in Applied Probability. *pp. 439-468. [ Links ]

McLennan, J.A., Ortiz, J.M., and Deutch, C.V. 2006. Geostatistical simulation of optimum mining elevations for nickel laterite deposits. *CIM Bulletin. *vol. 1, no. 6. pp. 1-9. [ Links ]

Morgan, C.M. 1992. Geology of the spheres at Weipa. *Travaux de Comité Internationale pour l'Etude des Bauxites, de Alumine et de l'Aluminium, *vol. 22. p. 61-74. [ Links ]

Myers, D.E. 1989. To be or not to be...stationary? That is the question. *Mathematical Geology, *vol. 21. pp. 343-362. [ Links ]

Neuman, S.P. and Jacobson, E.A. 1984. Analysis of nonintrinsic spatial variability by residual kriging with application to regional groundwater levels. *Journal of the International Association for Mathematical Geology, *vol. 16. pp. 499-521. doi: 10.1007/Bf01886329 [ Links ]

Odeh, I., McBratney, A., and Slater, B. 1997. Predicting soil properties from ancillary information: non-spatial models compared with geostatistical and combined methods. *Proceedings of Geostatistics Wollongong 96. *Vol. 1. Baafi, E.Y. and Scofield, N.A. (eds). Kluwer. pp. 1008-1019. [ Links ]

Odeh, I.O.A., McBratney, A.B., and Chittleborough, D.J. 1994. Spatial prediction of soil properties from landform attributes derived from a digitial elevation model. *Geoderma, *vol. 63. pp. 197-214. [ Links ]

Olea, R.A. 1974. Optimal contour mapping using universal kriging. *Journal of Geophysical Research, *vol. 79. pp. 695-702. doi: 10.1029/JB079i005p00695 [ Links ]

Oliver, M.A. 2010. Geostatistical Applications for Precision Agriculture. Springer. [ Links ]

Oliver, M.A. and Webster, R. 2015. Basic Steps in Geostatistics : the Variogram and Kriging. Springer. [ Links ]

Rossi, M.E. and Deutsch, C.V. 2014. Mineral Resource Estimation. Springer. [ Links ]

Singh, N.K. 2007. Ground penetrating radar (GPR) 'mineral base profiling and orebody optimisation'. *Proceedings of the 6th International Heavy Minerals Conference 'Back to Basics', *Hluhluwe, South Africa. Southern African Institute of Mining and Metallurgy, Johannesburg. pp. 179-182. [ Links ]

Wackernagel, H. 2002. Multivariate Geostatistics. Springer, Berlin. [ Links ]

WATTS, A. 1997. Exploring for nickel in the 90s, or 'Til depth us do part'. Geophysics and Geochemistry at the Millenium. *Proceedings of Exploration 97: Fourth Decennial International Conference on Mineral Exploration, *Toronto, Canada. Gubins, A.G. (ed.) Propspectore and Developers Association of Canada. pp. 1003-1013. [ Links ]

Webster, R. and Oliver, M.A. 2001 Geostatistics for environmental scientists. *Statistics in Practice. *Wiley, Chichester, UK. [ Links ]

Zimmerman, D., Pavlik, C., Ruggles, A., and Armstrong, M.P. 1999. An experimental comparison of ordinary and universal kriging and inverse distance weighting. *Mathematical Geology, *vol. 31. pp. 375-390. [ Links ] ♦

]]> Paper received Apr. 2016

Revised paper received Aug. 2017 ]]>