| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
|
|
|
|||||||
|
|||||||||
Australian Cotton Cooperative Research Centre, Dep. of Agricultural Chemistry and Soil Science, The Univ. of Sydney, NSW 2006, Australia
Corresponding author (johnt{at}acss.usyd.edu.au)
 |
ABSTRACT |
|---|
 TOP ABSTRACT INTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
Abbreviations: BLUE, best linear unbiased estimates • Co-K, cokriging • ECa, soil electrical conductivity • ECe, soil salinity • EM0,H, soil conductivity measurements made in horizontal mode • EM0,V, soil conductivity measurements made in vertical mode • ME, mean error • OK1, ordinary kriging of raw EM0,H data • OK2, ordinary kriging based on ECe estimates • RK, regression kriging • RMSE, root mean square error • SDR, standard deviation of ranks • 3-DK, three-dimensional kriging
|
INTRODUCTION |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
Cameron et al. (1981) and van der Lelij (1983) used an EM38 to describe the spatial variability of salinity at the field scale, producing choropleth maps showing areas of low, medium, and high average salinity. Unfortunately, the final maps contain inherent errors in their classification because of high estimation errors at unsampled locals. The method also assumed random and independent variation within and between mapping units (Chang et al., 1988). A more appropriate method, in terms of information transfer, is the use of isolinear maps where points of equal value are joined so that continuous representation is achieved with minimal information loss. This is significant, particularly in flat terrains, where landscape features or geological processes are subtle.
Within the last two decades, geostatistical methods, such as kriging, have been introduced into soil science to provide BLUE at unsampled locations (Burgess and Webster, 1980; Odeh et al., 1995). The major difference between geostatistics and classical statistics is that the former allows for the direct modeling of the inherent spatial data correlation. This is achieved through the initial calculation of a variogram, which acts as a quantified summary of all the available structural information of one or more random functions (Journel and Huijbregts, 1978). In comparing classical statistical and geostatistical methods, Hajrasuliha et al. (1980) found that the spatial dependent nature of the variance structure in some fields allowed geostatistical techniques to produce better predictions of soil salinity. Gallichand et al. (1992) produced similar results when comparing geostatistical methods with moving and weighted moving average techniques. The results suggested that not only were the resultant kriged maps more accurate, but more readily interpretable.
Lesch et al. (1992) mapped spatial variation of soil salinity on the field scale using geostatistical approaches. However, where there is a lack of spatial autocorrelation, a statistical calibration technique based on multiple-linear regression for predicting multiple depth, field scale spatial distribution of soil ECa is more appropriate (Lesch et al., 1995a).
To date, there has been no comparison of geostatistical methods to determine which may be an optimal approach. In this study, a comparison was made of a number of methods, including: ordinary kriging, cokriging, RK, and 3-DK. Similarly, no studies have been carried out to determine whether interpolation of raw ECa data was more accurate than first converting these values, via suitable calibrations, to soil electrical conductivity of a saturated soil paste extract ECe. Further comparisons were made from these two approaches using various methods of kriging. In all cases, the calibration approach of Triantafilis et al. (2000) was used to determine ECe at validation sites.
|
REVIEW OF STATISTICAL METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
 | (1) |
the error term with zero mean and constant variance. The regression coefficients, a and b, are not known a priori, but could be estimated by ordinary least-squares methods. Our interest here focuses on how Eq. [1] can be used to predict z at a location s0 where only xr (a row of predictor variables) has been determined.
Ordinary Kriging
Ordinary kriging is one of the most basic of kriging methods. It provides an estimate at an unobserved location of variable z, based on the weighted average of adjacent observed sites within a given area. The theory is derived from that of regionalized variables (Matheron, 1965, 1971) and can be briefly described by considering an intrinsic random function denoted by Z(si), where si represents all sample locations, i = 1,...,n. An estimate of the weighted average given by the ordinary kriging predictor at an unsampled site, z(s0), is defined by:
 | (2) |
i are the weights assigned to each of the observed sample sites. These weights sum to unity so that the predictor provides an unbiased estimation:
 | (3) |
The weights are calculated from the matrix equation
 | (4) |
.
Cokriging
Cokriging is an extension of kriging, and is a method of estimating one or more variables of interest using data from several variables by incorporating not only spatial correlation but also intervariable correlation. As such Co-K is the most versatile and rigorous statistical technique for spatial point estimation when both primary and secondary (covariate) attributes are available. It has been used widely in soil science (Vauclin et al., 1983; Trangmar et al., 1986; Yates and Warrick, 1987; Vaughan et al., 1995).
In defining Co-K, consider that Zj(s),..., Zm(s) are random functions denoting the values of the variables at location s in 1,..., n space, with m being the number of random variables (Myers, 1982). Cokriging is therefore defined as:
 | (5) |
If each component of z(s0) satisfies the Intrinsic Hypothesis (Journel and Huijbregts, 1978), then Eq. [5] is unbiased if
 | (6) |
 | (7) |
(si, sj) and (si, s0) are the cross-variograms, 
j are the weights associated with prediction, and is the Lagrange Multiplier for i = 1,..., n.
Regression Kriging
Regression kriging involves various combinations of linear regressions and kriging. The simplest model is based on normal regression followed by ordinary kriging with regression residuals (Odeh et al., 1995). The difference being that the ordinary kriging model is modified by replacing the variances in the diagonal of the A matrix with error terms that represent uncertainty. This is equivalent to "kriging with uncertain data" as defined by Ahmed and DeMarsily (1987). The method is a realization of an intrinsic random function Z(s), which is uncertain, the uncertainty giving rise to:
 | (8) |
i2. This uncertainty is due to errors in measurement or regression errors. These errors are assumed to be systematic or unbiased (E[e(si)] = 0), uncorrelated among themselves (Cov[e(si), e(si)] = 0), and uncorrelated with the variables (Cov[e(si), z(si)] = 0), i =1,..., n). The kriging equation is modified accordingly:
 | (9) |
Three-Dimensional Kriging
Three-dimensional kriging is an extension of ordinary kriging to a third dimension (i.e., depth). Despite its potential application in soil science, few examples of its use exist. Of these examples, Gallichand et al. (1992) were interested in the spatial interpolation of soil salinity and sodicity data. They found that spatial interpolation of soil salinity data with 3-DK was precise, robust, and led to the production of good quality maps. Similarly, Pan et al. (1992) found that 3-D Co-K, as compared with ordinary kriging, improved estimation significantly in terms of discrepancies between true and estimated values.
|
MATERIALS AND METHODS |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
651 ha. The northern fields (i.e., 16–20) are located around Galathera Creek. The southern fields (i.e., 23–25) are adjacent to a large water storage and were constructed using laser levelling.
|
EM38 Field Survey
The EM38 instrument has been used to measure soil ECa at the field scale, with sampling intervals ranging between 5 to 100 m (Cameron et al., 1981; Lesch et al., 1992; Hendrickx et al., 1992; Lesch et al., 1995b; Vaughan et al., 1995; Ceuppens et al., 1997). In order to determine an appropriate survey spacing, single east–west transects were laid in Fields 16, 17, 23, and 24. Soil ECa measurements were made in vertical (EM0,V) and horizontal (EM0,H) modes with the EM38. The instrument was placed on the ground and in the furrow at 10 m-intervals along these transects. This was carried out shortly after an irrigation event.
The variogram parameters of the ECa transect data are presented in Table 1. The variograms showed a range of 350 m. In order to adequately account for the spatial distribution, and for practicality, an ECa sampling spacing of 50 m was adopted. Surveying was undertaken shortly after an irrigation event, beginning from the southeast corner of each field. Measurements in each row were made just beyond the head ditch and 20 m from this point. Readings were then obtained every 50 m along the row until at a distance of 20 m from the end; another measurement was made before a final 1 m from the end. A measuring wheel was used to determine the distance. At the completion of each row, another transect was selected 50 m north of the previous transect. Where salinity was perceived to be greater (i.e., larger ECa reading), a sampling interval ranging from 12.5 to 25 m was used. To ensure geodatic accuracy at submeter level, a Global Positioning System (GPS NAVPRO5000, Magellan Systems, San Dimas, CA) equipped with an FM receiver (RDS 1000) was used to determine spatial coordinates in the Australian Map Grid (National Mapping Council of Australia, 1972).
|
|
|
EM38 Calibration
The methodology used to predict ECe at discrete depth increments using EM0,H measurements is described in more detail in Triantafilis et al. (2000). Briefly, it involved fitting a mixed random and fixed effects logistic model to each of the 26 calibration holes. The salinity profiles were modeled as smooth curves, with deterministic and random components, which can be expressed by the following equation (Eq. [13] of Triantafilis et al., 2000):
 | (10) |
Spatial Prediction
In comparing all the spatial prediction models, a question in mind was whether interpolation of the raw ECa data should be done, or whether to convert the data into estimates of ECe at various depths before interpolation. The implications of this is whether the extra time needed to manipulate the data in the latter approach warrants the effort in terms of improved prediction. Some interpolation methods, such as ordinary kriging, involved two approaches: (i) kriging the raw ECa data before estimating ECe, and (ii) converting the ECa data into ECe prior to kriging. In addition, and despite the fact that the calibration approach of Triantafilis et al. (2000) only requires EM0,H data, another question which needs answering is whether Co-K of this variable and EM0,V improves prediction? Further still, does the inclusion of residuals in the regression-kriging model improve prediction?
In kriging the EM0,H data prior to the ECe prediction, it is assumed that the predicted ECe values are unbiased based on BLUE. This is not entirely true. This is because the logistic calibration model of Triantafilis et al. (2000) was designed to be used with known EM0,H survey data. Similarly, by converting the EM0,H data (using the logistic function) prior to kriging of ECe, some bias in the estimation is introduced. As a result, the true variograms of ECe can not be reasonably estimated. Considering these facts, the two approaches could be regarded as semiempirical. However, this does not detract from the applicability of the five different approaches that were carried out.
Five different geostatistical approaches were carried out.
Ordinary Kriging (OK1)
Variogram of EM0,H was calculated and data interpolated using ordinary kriging. Calibration approach of Triantafilis et al. (2000) (i.e., Eq. [10]) used to estimate ECe at successive 0.30 m depth increments.
Ordinary Kriging (OK2)
Predictions of ECe were made at successive 0.30 m depth increments at each ECa survey site (using Eq. [10]). Variograms of ECe were calculated and data interpolated using ordinary kriging.
Cokriging (Co-K)
Covariogram calculated for EM0,V and EM0,H and used for cokriging. Prediction of EM0,H at validation site was converted into estimates of ECe at 0.30 m increments (using Eq. [10]).
Regression Kriging (RK)
Data converted for OK2 used for RK. Residuals were calculated and included for each 0.30-m depth increment from regression analysis carried out on the predicted ECe (using Eq. [10]) and measured ECe (using Eq. [3]). The program used for interpolation was KRIS-2D (McBratney et al., 1991) as modified by Odeh et al. (1994). The variogram parameters used were those derived from the ISATIS (1994) platform for OK2.
Three-Dimensional Kriging (3-DK)
Data converted for OK2 used for 3-DK. A variogram was calculated in the third dimension. Three-dimensional kriging was carried out using the ISATIS (1994) platform. All the variograms and cross variograms were generated within the ISATIS (1994) platform.
Validation
As a means of statistically testing the methods of spatial prediction described above, two indices were calculated from the observed and predicted values as used by Voltz and Webster (1990) at each validation site. As there were l sites (i.e., 53 sites) belonging to the validation sample set, the indices determined from the observed values z(sj) and the predicted values z*(sj) were the mean error (ME), and the root mean square error (RMSE). The ME measured the bias of prediction and for unbiased methods should be close to 0. It was defined as:
 | (11) |
The RMSE was a measure of how precise the various prediction methods were. It should be as small as possible for unbiased and precise prediction. It was expressed as
 | (12) |
ME and RMSE were calculated by considering predictions across all depths.
In some cases, the noise-to-signal ratios of the data may be large due to the presence of outliers that affect the final result. Laslett et al. (1987) suggested that a more revealing measure of performance can be obtained by computing the rank of each prediction method for every data point i, for i = 1,...l, based on the point value of prediction square error of a method. Thus rij was the rank of the ith method at the jth data point, so that the rank of the ith method was:
 | (13) |
 | (14) |
A plot of mean rank vs. SDR could be used to determine the optimal approach. The method that exhibited a low mean rank in addition to a low SDR would be optimal in terms of prediction performance. As a result, these measures cancelled out the effects of any large outliers in the data that could not be accounted for in the overall calculation of RMSE and ME, which provide general measures of prediction (Odeh et al., 1994, 1995).
|
RESULTS |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
|
The variogram parameters for the ECa and calibrated ECe data, using the approach of Triantafilis et al. (2000), are shown in Table 1. The ECa variograms were, in general, similar to those calculated for the preliminary-transect data, although the variance was greater. The reason for this is that the ECa survey data was obtained at a 50-m interval, as compared with a 10-m interval along east–west transects. As such, the latter incorporated more of the variation in the x–y plane than did the transect survey. This might explain why the range was larger with the areal survey variograms than with the transect variograms. In general, the variograms of ECe at the various depths shows increasing nugget variance and sill with increasing depth, although as a proportion of the sill, the spatially-dependent component also showed an increase. High prediction errors may result, especially in the surface layers. By comparison, the smaller values of the spatially-dependent components suggest that slightly larger prediction errors were incorporated into the kriging system when the data was first transformed into ECe estimates prior to interpolation.
The ECa data was generally found to be quasi-stationary in the x–y plane. This was similarly the case for each of the ECe variograms calculated at each 0.3 m depth. As a result, vertical–spatial isotropy was assumed for each of the methods of interpolation, except in the three-dimensional case. For each of the two-dimensional interpolations, spherical models were fitted to the experimental variograms. In calculating the variogram in the vertical or z direction for the calibrated EM0,H survey data, it was apparent that a drift exists. This is due to an increasing trend of salinity to a depth of 2.0 m. This drift was used for 3-DK (ISATIS, 1994).
Statistical Validation
In total, n = 3363 data points were used for the purpose of spatial prediction of the ECa survey data. Another smaller number of validation points (l = 53), independent of the model data, were used to test the geostatistical methods. First the RMSE and ME were used to compare the methods as shown in Fig. 4
. With respect to the ME, all methods were close to zero (i.e., |ME| < 1). This was not an unusual result, considering the unbiased nature of the geostatistical methods. It should be noted that the negative ME suggested that the methods of prediction used here was overestimating the ECe within the profile (i.e., observed < predicted: see Eq. [11]). This overestimation is also shown in Table 3, where average measured ECe is 3.33 dS m-1, as compared with predicted ECe, ranging between 3.49 and 3.52 dS m-1.
|
|
|
|
In order to determine whether differences in precision and bias existed between the methods of interpolation in different fields, we compared the RMSE and ME from Fields 18 and 20 independently (Fig. 6) . Figure 6a shows that in Field 18, the bias was positive (i.e., observed > predicted), and hence soil ECe was underestimated. Conversely, Fig. 6b shows that in Field 20 the bias was negative (i.e., observed < predicted), and hence overestimating soil ECe. These results were consistent with the fact that the EM38 survey was undertaken in Field 18 when soil moisture was below field capacity, and in Field 20 when the soil condition was slightly greater than field capacity, as compared with the calibrated soil samples. Field 18 was surveyed soon after a heavy rainfall and during a fallow season. This means that when the calibration samples were taken (i.e., several days after an irrigation event), the soil condition, with respect to field capacity, was suboptimal. Conversely, Field 20 was surveyed at the beginning of an irrigation season and 3 d after a heavy rainfall event, and therefore would have been above field capacity. As a consequence, the estimates are therefore slightly greater than the observed soil ECe. With respect to precision, there was little difference between the methods of RK, Co-K, OK1, and OK2 in Field 18. As with the overall precision of 3D-K, in Field 18 it was the worst performer. Conversely, in Field 20, 3D-K was the most precise, although there was little difference between this approach and that of OK1, OK2, and Co-K. The worst performing method was RK.
|
|
0.6 m km-1 , the water flowing in these creeks on the plains lose energy (Ward, 1999). This is evidenced by the meandering nature of Galathera Creek on the plains proper. Just beyond the study area and to the west, the channel itself becomes ill-defined as the water is dispensed, fanning out over a small area. As a result, dissolved salts have also been dispensed and subsequently accumulated in these areas, and in particular at depths between 0.9 to 2.0 m along old drainage lines and previously lower-lying areas. This is the case in the northeast-facing arc of higher ECe shown in Fields 19 and 20 (Fig. 7a), and the higher areas of ECe in the north and western part of Field 18 (Fig. 7b), respectively. |
DISCUSSION AND CONCLUSIONS |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
By comparison, each of the other approaches (i.e., OK2, 3-DK, and RK) proceeded only after ECa was first converted to ECe at 0.30-m depth increments. Ordinary kriging (OK2) of the ECe estimates required calculation and modeling of a separate variogram for each of the various depths of interest (0–0.3, 0.3–0.6 ... 1.8–2.0). A total of 21 parameters were required, three for each depth. Therefore, seven separate interpolations were necessary to predict ECe at each validation site. To overcome this problem, 3-DK of the converted ECe data used in OK2 was used. This required a third variogram to be calculated in the vertical direction, so that 6 parameters were calculated altogether (i.e., 3 in the x–y plane and 3 in the vertical plane). Unfortunately, the approach only improved prediction in Field 20.
The RK approach, which required regression residuals to be incorporated within the kriging system at each ECe depth, performed best overall, in terms of precision and bias. The cost vs. benefit balance of the methods was therefore more important. The user must consider whether the cost required in preparing the data to predict ECe using RK is worth the extra precision. In all, 35 parameters were calculated (the 3 by 7 variogram parameters calculated for 7 depth increments, as well as the two additional regression parameters, a and b, of each line (i.e., 14). Today, these issues are not as relevant as they may have been a decade ago, given the increased computational speeds of computers and the software that runs them.
In summary, the results obtained here are equivocal based on these conclusions. The reason for this may be a result of an insufficient number of validation samples being taken in one particular field, and in others, not enough were generated across a broad range of soil ECe values. In addition, some of the work carried out here was at times conducted at less than ideal soil moisture conditions. Large errors may have been the result if moisture content was well below field capacity due to different soil or irrigation management. This was shown to be significant when the EM38 instrument is used in the horizontal mode of operation, as discussed by Vaughan et al. (1995). Caution is therefore required in timing field surveys. However, and judging by the overall low RMSE and ME, the results appeared to be free of any significant cultural effects or variability in moisture condition. This was consistent with the results of Hendrickx et al. (1992), who found errors from this effect in irrigated soil profiles was minimal.
Where the problem of moisture content or site selection is considered to be potentially significant, we suggest the approach of Lesch et al. (1995b) be adopted. In this case, both soil calibration and validation sites will be taken in the field of interest. Most of the issues of differences in moisture conditions caused by seasonal effects, such as soil or crop management within the field, would be minimized. Of course, the cost vs. benefit balance of such a detailed approach still needs to be explored.
|
ACKNOWLEDGMENTS |
|---|
Received for publication May 9, 2000.
|
REFERENCES |
|---|
|
TOPABSTRACTINTRODUCTIONREVIEW OF STATISTICAL METHODSMATERIALS AND METHODSRESULTSDISCUSSION AND CONCLUSIONSREFERENCES |
|---|
This article has been cited by other articles:
|
|
 |
|
  T.-L. Liu, K.-W. Juang, and D.-Y. Lee Interpolating Soil Properties Using Kriging Combined with Categorical Information of Soil Maps Soil Sci. Soc. Am. J., May 23, 2006; 70(4): 1200 - 1209. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 Y. Ouyang, J. E. Zhang, and L.-T. Ou Temporal and Spatial Distributions of Sediment Total Organic Carbon in an Estuary River J. Environ. Qual., January 3, 2006; 35(1): 93 - 100. [Abstract] [Full Text] [PDF]
|
|
|
|
 |
|
 A. B. Tarr, K. J. Moore, and P. M. Dixon Spectral Reflectance as a Covariate for Estimating Pasture Productivity and Composition Crop Sci., May 6, 2005; 45(3): 996 - 1003. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
Y. Ouyang, L.-T. Ou, and G. C. Sigua Characterization of the Pesticide Chlordane in Estuarine River Sediments J. Environ. Qual., March 1, 2005; 34(2): 544 - 551. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
S. Ersahin Comparing Ordinary Kriging and Cokriging to Estimate Infiltration Rate Soil Sci. Soc. Am. J., November 1, 2003; 67(6): 1848 - 1855. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
Y. Ouyang, P. Nkedi-Kizza, R. S. Mansell, and J. Y. Ren Spatial Distribution of DDT in Sediments from Estuarine Rivers of Central Florida J. Environ. Qual., September 1, 2003; 32(5): 1710 - 1716. [Abstract] [Full Text] [PDF]
|
|
|
|
|
|
T. G. Mueller and F. J. Pierce Soil Carbon Maps: Enhancing Spatial Estimates with Simple Terrain Attributes at Multiple Scales Soil Sci. Soc. Am. J., January 1, 2003; 67(1): 258 - 267. [Abstract] [Full Text] [PDF]
|
|
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
