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

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 sparselyspaced borehole data, as the estimates would not reflect the actual variability in the surface. Therefore, it is imperative that the sparselyspaced borehole data be supplemented with appropriate geophysical information for more accurate resource estimations (Erten et al., 2015, 2013). Ground-penetrating radar (GPR) 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 lateritetype 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 ordertwo 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 Accounting for a spatial trend in finescale ground-penetrating radar data: a comparative case study

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 sparselyspaced borehole data, as the estimates would not reflect the actual variability in the surface.Therefore, it is imperative that the sparselyspaced borehole data be supplemented with appropriate geophysical information for more accurate resource estimations (Erten et al., 2015(Erten et al., , 2013)).Ground-penetrating radar (GPR) 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 lateritetype 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 ordertwo 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 Accounting for a special trend in fine-scale ground-penetrating radar data 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 nonsampled 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.UK, IRF-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.
The RF can be represented by the following model: [1] 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-order stationarity 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: [3] 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: [4] where (h) represents the variogram at lag h.The following relationship between the variogram (h) and covariance C(h) function exists for a RF which is second-order stationary: [5] where C(0) represents the variance 2 of the RF.
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). [6] 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 (u) are determined to satisfy this unbiasedness and minimized variance goals.The following OK system is used to obtain the kriging weights: [7] 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: In UK, the trend component is modelled as a smoothly varying deterministic function of u and is expressed as: [9] where m(u) represents the local mean, l , the unknown coefficients of the trend function, and f l (u) represent the known functions of the coordinates u and are called trial or base functions (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: [10] 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: [11] where the set of L constraints is termed the universality or unbiasedness conditions.The kriging system satisfying these requirements is then defined by: [12] where are the UK weights and l are the Lagrange parameters used to determine the coefficients minimizing the error variance.
The UK estimate variance is defined as: [13] The inference of the residual variogram 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: [14] 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 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 Accounting for a special trend in fine-scale ground-penetrating radar data 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(u) in Equation [9] are restricted to those which are only translation-invariant and pairwise orthogonal (Wackernagel, 2002).
Considering a set of weights applied to particular points u , a discrete measure u is defined as: [15] where u represents the Dirac measure at point u .Any linear combinations of the weights with RV at locations u are defined as: The estimation variance is defined by: [21] The region where the mine site is located is composed of Proterozoic and Palaeozoic basement in the eastern part.This basement comprises acid intrusives, extrusives, and metamorphics.Overlying the basement to the west, Mesozoic and Cenozoic sediments dominate.The sediments were intensively weathered, which played a crucial role in the formation of laterites rich in alumina and bauxites.The bauxite sits above an almost entirely kaolinized pallid zone (Morgan, 1992).
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).
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.
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 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 m 2 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.
The target variable to be regionalized in this study is the Accounting for a special trend in fine-scale ground-penetrating radar data Accounting for a special trend in fine-scale ground-penetrating radar data 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.
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  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.
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 l (u) and the covariance C R (h) of the residual component R(u) inferred from the residual variogram R(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.
The first step in determining the underlying variogram for UK was to estimate the coefficients of the drift function with a least-squares based estimator.Once the coefficients of the drift function are determined, the residuals are computed by subtracting the drift from the data.An experimental variogram of the residuals is then calculated and a model is fitted to the residuals experimental variogram.The bias associated with the experimental variogram is computed and an iteration is applied to compute the corrected experimental variogram.This iteration is carried out 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.
The second approach used to account for the nonstationarity was the IRF-k method.Nonstationary modelling with the IRFk 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 'Jacknife' 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).
In the OK case, the residual variogram R (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: 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 2 (u) is the kriging variance.The results of the crossvalidation 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.
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 twodimensional grid having the following properties was created: the origin is X 0 = -1872.00m, 0 = 10181.00m, 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.
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.
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.
In addition to these nonstationary methods, OK was implemented by neglecting the nonstationarity present in the data.In this method, the variogram model obtained along the trend-free direction was used as a structural input describing the spatial relationship.The maps of the results produced by OK are given in Figure 11.
Accounting for a special trend in fine-scale ground-penetrating radar data L  variable at an unknown grid node due to the neighbourhood definition.Since 32 samples were almost always within a 80-100 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.
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.
[16] The expression shown in Equation [16] is called the allowable linear combination of order k(ALC-k) in Equation [17] [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 IRFk is a symmetric function, K(h)=K(-h), and satisfies the condition: [18] An example to GC function is known the as polynomial generalized covariance function.The equation of the polynomial generalized covariance function is given as: [19] The conditions on b p are satisfied if b pp.The intrinsic kriging system minimizing the variance is expressed in terms of the GC and shown in the following:[20] , the variogram values in the N70° direction are higher only for the lag values L 178 VOLUME 118

Table I
.