A practical two‐step approach for mixed model‐based kriging,with an application to the prediction of soil organic carbon concentration

Soil scientists often use prediction models to obtain values at unsampled locations. The spatial variation in the soil is best captured by using the empirical best linear unbiased predictor (EBLUP) based on a restricted maximum likelihood (REML) approach that efficiently exploits available data on both mean trends and correlation structures. We proposed a practical two‐step implementation of the REML approach for model‐based kriging, exemplified by predicting soil organic carbon (SOC) concentrations in mineral soils in Estonia from the large‐scale digital soil map information and a previously established prediction model. The prediction model was a linear mixed model with soil type, physical clay content (particle size < 0.01 mm) and A‐horizon thickness as fixed effects and site, transect, plot, year, year‐transect random intercepts and site‐specific random slopes for clay content. We used only the site‐specific intercept EBLUPs for estimating spatial correlation parameters as they described most of the variation in the random effects (86.8%). Fitting an exponential correlation model to these EBLUPs resulted in an estimated range of 10.5 km and the estimated proportion of the variance from the nugget effect was 0.23. The results of a simulation study showed a downwards bias that decreased with sample size. The results were validated through an external dataset, resulting in root mean square errors (RMSE) of 1.06 and 1.07% for the two‐step approach for kriging and the model with only fixed effects (no kriging), respectively. These results indicate that using the two‐step approach for kriging may improve prediction.     

Optimal interpolation and isarithmic mapping of soil properties: I The semi‐variogram and punctual kriging

Kriging is a means of spatial prediction that can be used for soil properties. It is a form of weighted local averaging. It is optimal in the sense that it provides estimates of values at unrecorded places without bias and with minimum and known variance. Isarithmic maps made by kriging are alternatives to conventional soil maps where properties can be measured at close spacings. Kriging depends on first computing an accurate semi‐variogram, which measures the nature of spatial dependence for the property. Estimates of semi‐variance are then used to determine the weights applied to the data when computing the averages, and are presented in the kriging equations. The method is applied to three sets of data from detailed soil surveys in Central Wales and Norfolk. Sodium content at Plas Gogerddan was shown to vary isotropically with a linear semi‐variogram. Ordinary punctual kriging produced a map with intricate isarithms and fairly large estimation variance, attributed to a large nugget effect. Stoniness on the same land varied anisotropically with a linear semi‐variogram, and again the estimation error of punctual kriging was fairly large. At Hole Farm, Norfolk, the thickness of cover loam varied isotropically, but with a spherical semi‐variogram. Its parameters were estimated and used to krige point values and produce a map showing substantial short‐range variation.     

Estimation and Prediction in the Spatial Linear Model

Often in environmental monitoring studies interesting ecological factors will be observed at several locations repeatedly over time. Generally these space-time data are subject to a sequential spatial data analysis. In geostatistics, spatial data describing an environmental phenomenon like the pH value in precipitation at several locations are regarded as a realisation from a stochastic process. Component models are used to interpret the spatial variation of the process. Decomposing the spatial process into single components is based on the theory of linear models. Trend surface analysis is seen to be the geostatistical method for best linear unbiased estimation (BLUE) of the trend component, whereas universal kriging is equivalent to best linear unbiased prediction (BLUP) of the realisation of the spatial process. Furthermore trend surface analysis and universal kriging are shown to agree with the estimation of fixed effects and prediction of fixed and random effects in mixed linear models. Since estimation and prediction for spatial data result in different interpolations the differences are explained also graphically by example. The example uses acid-precipitation monitoring data. The extension of these spatial methods for application to space-time problems by combination with dynamic linear models is treated in the discussion.     

Prediction of spatial soil property information from ancillary sensor data using ordinary linear regression: Model derivations, residual assumptions and model validation tests

Geospatial measurements of ancillary sensor data, such as bulk soil electrical conductivity or remotely sensed imagery data, are commonly used to characterize spatial variation in soil or crop properties. Geostatistical techniques like kriging with external drift or regression kriging are often used to calibrate geospatial sensor data to specific soil or crop properties. More traditional statistical methods such as ordinary linear regression models are also commonly used. Unfortunately, some soil scientists see these as competing and unrelated modeling approaches and are unaware of their relationship. In this article we review the connection between the ordinary linear regression model and the more comprehensive geostatistical mixed linear model and describe when and under what conditions ordinary linear regression models represent valid spatial prediction models. The formulas for the ordinary linear regression model parameter estimates and best linear unbiased predictions are derived from the geostatistical mixed linear model under two different residual error assumptions; i.e., strictly uncorrelated (SU) residuals and effectively uncorrelated (EU) residuals. The theoretically optimal (best linear unbiased) and computable (linear unbiased) predictions and variance estimates derived under the EU error assumption are examined in detail. Statistical tests for detecting spatial correlation in LR model residuals are also reviewed, in addition to three LR model validation tests derived from classical linear modeling theory. Two case studies are presented that highlight and demonstrate the various parameter estimation, response variable prediction and model validation techniques discussed in this article.     

Factorial kriging analysis: a useful tool for exploring the structure of multivariate spatial soil information

Most studies of relations between soil properties fail to take account of their regionalized nature because of the lack of appropriate methods. This paper describes a geostatistical technique, factorial kriging analysis, that bridges the gap between classical multivariate analysis and a univariate geostatistical approach. The basic feature of the method is the fitting of a linear model of coregionalization, i.e. all experimental simple and cross-variograms are modelled with a linear combination of basic variogram functions. A particular variance-covariance matrix, the coregionalization matrix, can then be associated with each spatial scale defined by the range of the basic variogram function. Each coregionalization matrix describes relationships between variables at a given spatial scale. A principal component analysis of these matrices produces a set of components, the regionalized factors, that reflect the main features of the multivariate information for each spatial scale and whose scores are estimated by cokriging. The technique is described and illustrated with three case studies based on a simulated data set and soil survey data. The results are compared with those of the principal component analysis of the variance-covariance matrix and the variogram matrices.     

Comparing sampling needs for variograms of soil properties computed by the method of moments and residual maximum likelihood

It has been generally accepted that the method of moments (MoM) variogram, which has been widely applied in soil science, requires about 100 sites at an appropriate interval apart to describe the variation adequately. This sample size is often larger than can be afforded for soil surveys of agricultural fields or contaminated sites. Furthermore, it might be a much larger sample size than is needed where the scale of variation is large. A possible alternative in such situations is the residual maximum likelihood (REML) variogram because fewer data appear to be required. The REML method is parametric and is considered reliable where there is trend in the data because it is based on generalized increments that filter trend out and only the covariance parameters are estimated. Previous research has suggested that fewer data are needed to compute a reliable variogram using a maximum likelihood approach such as REML, however, the results can vary according to the nature of the spatial variation. There remain issues to examine: how many fewer data can be used, how should the sampling sites be distributed over the site of interest, and how do different degrees of spatial variation affect the data requirements? The soil of four field sites of different size, physiography, parent material and soil type was sampled intensively, and MoM and REML variograms were calculated for clay content. The data were then sub-sampled to give different sample sizes and distributions of sites and the variograms were computed again. The model parameters for the sets of variograms for each site were used for cross-validation. Predictions based on REML variograms were generally more accurate than those from MoM variograms with fewer than 100 sampling sites. A sample size of around 50 sites at an appropriate distance apart, possibly determined from variograms of ancillary data, appears adequate to compute REML variograms for kriging soil properties for precision agriculture and contaminated sites.     

Non-homogeneity of variance components from spatially nested sampling of the soil

Spatially nested sampling and the associated nested analysis of variance by spatial scale is a well-established methodology for the exploratory investigation of soil variation over multiple, disparate scales. The variance components that can be estimated this way can be accumulated to approximate the variogram. This allows us to identify the important scales of variation, and the general form of the spatial dependence, in order to plan more detailed sampling by design-based or model-based methods. Implicit in the standard analyses of nested sample data is the assumption of homogeneity in the variance, i.e. that all variations from sub-station means at some scale represent a random variable of uniform variance. If this assumption fails then the comparable assumption of stationarity in the variance, which is an important assumption in geostatistics, will also be implausible. However, data from nested sampling may be analysed with a linear mixed model in which the variance components are parameters which can be estimated by residual maximum likelihood (REML). Within this framework it is possible to propose an alternative variance parameterization in which the variance depends on some auxiliary variable, and so is not generally homogeneous. In this paper we demonstrate this approach, using data from nested sampling of chemical and biogeochemical soil properties across a region in central England, and use land use as our auxiliary variable to model non-homogeneous variance components. We show how the REML analysis allows us to make inferences about the need for a non-homogeneous model. Variances of soil pH and cation exchange capacity at different scales differ between these land uses, but a homogeneous variance model is preferable to such non-homogeneous models for the variance of soil urease activity at standard concentrations of urea.     

A stochastic-geometric model of soil variation

Stochastic models of soil variation are used in geostatistical analysis, but in general they bear no relation to our mechanistic understanding of the processes in soil that cause its properties to vary spatially. It is proposed that we require a suitable stochastic model in which space is partitioned into discrete domains as a first step towards random spatial models that incorporate our understanding of processes in soil. Even though the soil is essentially continuous in its spatial variation, there are components of soil variation (e.g. differences between parent materials) which are discontinuous. This paper shows how variogram models can be derived directly from the Poisson Voronoi Tessellation (PVT), a stochastic-geometric partition of d -dimensional space. The PVT variogram models, for d = 2 and 3, were fitted to variograms estimated from data over disparate scales, including computerized tomographic images of soil aggregates (pixels of a few tens of micrometres long) and the land systems of Swaziland. In all cases, PVT variogram models fitted better than the conventional geostatistical ones. The good performance of PVT variogram models at these disparate scales encourages further work on tessellation models for soil variation. In principle such models could incorporate information on underlying factors of soil formation such as the spatial distribution of individual plants, the origin and growth of microbial colonies, spatial processes in soil chemistry (such as reaction–diffusion processes) and geometrical information on boundaries between geological strata or contrasting plant communities. PVT models may therefore be one component of a random model of soil variation which reflects our understanding of soil-forming processes, and so have a stronger scientific basis than the models that are now in standard use.     

Model‐based analysis using REML for inference from systematically sampled data on soil

The general linear model encompasses statistical methods such as regression and analysis of variance (anova ) which are commonly used by soil scientists. The standard ordinary least squares (OLS) method for estimating the parameters of the general linear model is a design‐based method that requires that the data have been collected according to an appropriate randomized sample design. Soil data are often obtained by systematic sampling on transects or grids, so OLS methods are not appropriate. Parameters of the general linear model can be estimated from systematically sampled data by model‐based methods. Parameters of a model of the covariance structure of the error are estimated, then used to estimate the remaining parameters of the model with known variance. Residual maximum likelihood (REML) is the best way to estimate the variance parameters since it is unbiased. We present the REML solution to this problem. We then demonstrate how REML can be used to estimate parameters for regression and anova ‐type models using data from two systematic surveys of soil. We compare an efficient, gradient‐based implementation of REML (ASReml) with an implementation that uses simulated annealing. In general the results were very similar; where they differed the error covariance model had a spherical variogram function which can have local optima in its likelihood function. The simulated annealing results were better than the gradient method in this case because simulated annealing is good at escaping local optima.     

Distribution of Cd in the vicinity of a metal smelter: Interpolation of soil Cd concentrations with regard to regulative limits

The German soil protection regulation (BBodSchV) requires the investigation and evaluation of sites with known, or suspected contamination. The purpose of this study is the application of geostatistical methods to locate hazardous zones within such a site and to estimate the amount and uncertainty of the contaminant load in these zones. The study site is an area around a metal smelter in the city of Nordenham, Germany, where among other heavy metals, Cd was released to the environment by dust emissions for many decades. In an earlier study soil cores were taken in the area and analyzed for Cd using various extraction methods. After translation of data to results corresponding to a single extraction method using linear regression analysis, Cd concentrations were mapped by ordinary and lognormal kriging. Crossvalidation showed that both methods perform similarly. However, neither ordinary nor lognormal kriging were able to account for the uncertainty of the kriged estimates. We repeated ordinary kriging with a relative variogram having a unit sill. The estimated relative kriging variance was scaled locally. This method considerably improved the estimation of uncertainty. Subsequently, we estimated Cd contents for the land use dependent size of support as specified in the BBodSchV. The kriged Cd estimates as well as their uncertainty were evaluated with regard to limits set by the BBodSchV. Parts of the area which may be declared safe based on merely the kriged estimates, can actually exceed a sanction or test limit by a chance of up to 50 % when uncertainty is also considered. Within the BBodSchV a recommended limit should therefore always be accompanied by a tolerable uncertainty that it may be exceeded on a given support (e.g. 5 %).     

Analysis of Genebank Evaluation Data by using Geostatistical Methods

Genebanks often characterize accessions based on evaluation trials. This paper evaluates geostatistical methods as a tool to increase the utility of evaluation data. These methods were selected to overcome limitations resulting from a relative lack of replication and the scarcity of standards or check varieties. The data employed in the present study comprise nine characteristics of spring and winter barley, evaluated mostly as ratings. Ratings with quasi-metric scales were transformed by using the folded exponential transformation. To estimate the genetic component of the total effect, we compared two methods: Method 1 whereby a variogram is fitted by non-linear regression, and subsequently the implied spatial correlation is embedded into a mixed model analysis, which estimates the genetic effect by Best Linear Unbiased Prediction (BLUP); and Method 2 where each data value is re-estimated by kriging to correct for spatial effects and then the corrected data are submitted to a mixed model analysis. For practical application we propose Method 1 (though occasionally we met convergence problems): Fit the short range of the empirical variogram, visually choose the suitable covariance model. Use this and the initial values from non-linear regression fit with the mixed model, fixing the spatial parts at their starting values from non-linear regression, and estimate genetic effects by BLUP by using the fitted mixed model. To improve performance, we recommend that more standard or check varieties be used and, wherever possible, replace rating scales with metric scales or free-percentage scales (without categories).     

Kriging a soil variable with a simple nonstationary variance model

Kriging is a standard tool in the environmental sciences for spatial prediction from limited sample data, subject to the assumption of intrinsic stationarity, made about the underlying spatially correlated random function. It is generally well understood how the assumption of stationarity in the mean can be relaxed within the linear mixed model framework, using residual maximum likelihood to estimate variance parameters for the random effects. The Best Linear Unbiased Predictor (BLUP) is equivalent to the kriging predictor in these circumstances. However, nonstationarity in the variance is a harder problem to solve. Stationarity assumptions are necessary if the spatial covariance of a random process is to be estimated from the single realization which nature provides. However, they are not always plausible for variables arising from processes in complex landscapes across contrasting topography and geology.     

A coherent geostatistical approach for combining choropleth map and field data in the spatial interpolation of soil properties

Information available for mapping continuous soil attributes often includes point field data and choropleth maps (e.g. soil or geological maps) that model the spatial distribution of soil attributes as the juxtaposition of polygons (areas) with constant values. This paper presents two approaches to incorporate both point and areal data in the spatial interpolation of continuous soil attributes. In the first instance, area-to-point kriging is used to map the variability within soil units while ensuring the coherence of the prediction so that the average of disaggregated estimates is equal to the original areal datum. The resulting estimates are then used as local means in residual kriging. The second approach proceeds in one step and capitalizes on: 1) a general formulation of kriging that allows the combination of both point and areal data through the use of area-to-area, area-to-point, and point-to-point covariances in the kriging system, 2) the availability of GIS to discretize polygons of irregular shape and size, and 3) knowledge of the point-support variogram model that can be inferred directly from point measurements, thereby eliminating the need for deconvolution procedures. The two approaches are illustrated using the geological map and heavy metal concentrations recorded in the topsoil of the Swiss Jura. Sensitivity analysis indicates that the new procedures improve prediction over ordinary kriging and traditional residual kriging based on the assumption that the local mean is constant within each mapping unit.     

采样点数量对长三角典型地区土壤肥力指标空间变异解析的影响

通过对江苏省如皋市888个采样点的不重复随机抽样,探讨了采样点数量对土壤肥力指标空间变异解析的影响。从半方差函数估计的可靠性角度考虑,在长江冲积物形成的土壤上,针对县级农业管理和生态环境规划的土壤肥力指标调查采样,采集250个样点较为合适。另外,仅使用交互验证和独立验证评价半方差函数估计的可靠性及确定采样点数量是不完善的,而平均克里格方差理论上随采样点数量的增加而单调递减,可以作为不同采样点数量条件下,综合评价半方差函数估计可靠性及确定采样点数量的补充指标。     

A comparison of some robust estimators of the variogram for use in soil survey

The standard estimator of the variogram is sensitive to outlying data, a few of which can cause overestimation of the variogram. This will result in incorrect variances when estimating the value of a soil property by kriging or when designing a sampling grid to map the property to a required precision. Several robust estimators of the variogram, based on location and scale estimation, have been proposed as improvements. They seem to be suitable for analysis of soil data in circumstances where the standard estimator is likely to be affected by outliers. Robust estimators are based on assumptions about the distribution of the data which will not always hold and which need not be made in kriging or in estimating the variogram by the standard estimator. The estimators are reviewed. Simulation studies show that the robust estimators vary in their susceptibility to moderate skew in the underlying distribution, but that the effects of outliers are generally greater. The estimators are applied to some soil data, and the resulting variograms used for ordinary kriging at sites in a separate validation data set. In most cases the variograms derived from the standard estimator gave kriging variances which appeared to overestimate the mean squared error of prediction (MSEP). Kriging with variograms based on robust estimators sometimes gave kriging variances which underestimated the MSEP or did not differ significantly from it. Estimates of kriging variance and the MSEP derived from the validation data were generally close to estimates from cross‐validation on the prediction set used to derive the variograms. This indicates that variogram models derived from different estimators could be compared by cross‐validation.     

Estimating variograms of soil properties by the method-of-moments and maximum likelihood

Variograms of soil properties are usually obtained by estimating the variogram for distinct lag classes by the method‐of‐moments and fitting an appropriate model to the estimates. An alternative is to fit a model by maximum likelihood to data on the assumption that they are a realization of a multivariate Gaussian process. This paper compares the two using both simulation and real data. The method‐of‐moments and maximum likelihood were used to estimate the variograms of data simulated from stationary Gaussian processes. In one example, where the simulated field was sampled at different intensities, maximum likelihood estimation was consistently more efficient than the method‐of‐moments, but this result was not general and the relative performance of the methods depends on the form of the variogram. Where the nugget variance was relatively small and the correlation range of the data was large the method‐of‐moments was at an advantage and likewise in the presence of data from a contaminating distribution. When fields were simulated with positive skew this affected the results of both the method‐of‐moments and maximum likelihood. The two methods were used to estimate variograms from actual metal concentrations in topsoil in the Swiss Jura, and the variograms were used for kriging. Both estimators were susceptible to sampling problems which resulted in over‐ or underestimation of the variance of three of the metals by kriging. For four other metals the results for kriging using the variogram obtained by maximum likelihood were consistently closer to the theoretical expectation than the results for kriging with the variogram obtained by the method‐of‐moments, although the differences between the results using the two approaches were not significantly different from each other or from expectation. Soil scientists should use both procedures in their analysis and compare the results.     

应用土壤质地预测干旱区葡萄园土壤饱和导水率空间分布

田间表层土壤饱和导水率的空间变异性是影响灌溉水分入渗和土壤水分再分布的主要因素之一,研究土壤饱和导水率的空间变化规律,有助于定量估计土壤水分的空间分布和设计农田的精准灌溉管理制度。为了探究应用其他土壤性质如质地、容重、有机质预测土壤饱和导水率空间分布的可行性,试验在7.6 hm2的葡萄园内,采用均匀网格25 m×25 m与随机取样相结合的方式,测定了表层(0～10 cm)土壤饱和导水率、粘粒、粉粒、砂粒、容重和有机质含量,借助经典统计学和地统计学,分析了表层土壤饱和导水率的空间分布规律、与土壤属性的空间相关性,并对普通克里格法、回归法和回归克里格法预测土壤饱和导水率空间分布的结果进行了对比。结果表明:1)土壤饱和导水率具有较强的变异性,平均值为1.64 cm/d,变异系数为1.17;2)表层土壤饱和导水率60%的空间变化是由随机性或小于取样尺度的空间变异造成;3)土壤饱和导水率与粘粒、粉粒、砂粒和有机质含量具有一定空间相关性,而与土壤容重几乎没有空间相关性;4)在中值区以土壤属性辅助的回归克里格法对土壤饱和导水率的预测精度较好,在低值和高值区其与普通克里格法表现类似。研究结果将为更好地描述土壤饱和导水率空间变异结构及更准确地预测其空间分布提供参考。     

An error budget for digital soil mapping of cation exchange capacity using proximally sensed electromagnetic induction and remotely sensed γ‐ray spectrometer data

The cation exchange capacity (CEC) of soil is widely used for agricultural assessment as a measure of fertility and an indicator of structural stability; however, its measurement is time‐consuming. Although geostatistical methods have been used, a large number of samples must be collected. Using pedometric methods and incorporating easy‐to‐measure ancillary data into models have improved the efficiency of spatial prediction of soil CEC. However, mapping uncertainty has not been evaluated. In this study, we use an error budget procedure to quantify the relative contributions that model, input and covariate error make to prediction error of a digital map of CEC using gamma‐ray (γ‐ray) spectrometry and apparent electrical conductivity (EC_a) data. The error budget uses empirical best linear unbiased prediction (E‐BLUP) and conditional simulation to produce numerous realizations of the data and their underlying errors. Linear mixed models (LMMs) estimated by residual maximum likelihood (REML) are used to create the prediction models. The combined error of model [5.07 cmol(+)/kg] and input error [12.88 cmol(+)/kg] is ~12.93 cmol(+)/kg, which is twice as large as the standard deviation of CEC [6.8 cmol(+)/kg]. The individual covariate errors caused by the γ‐ray [9.64 cmol(+)/kg] and EM error [8.55 cmol(+)/kg] were large. Preprocessing techniques to improve the quality of the γ‐ray data could be considered, whereas the EM error could be reduced by the use of a smaller sampling interval in particular near the edges of the study area and at pedoderm boundaries.     

Two robust estimators of the cross‐variogram for multivariate geostatistical analysis of soil properties

If we wish to describe the coregionalization of two or more soil properties for estimation by cokriging then we must estimate and model their auto‐ and cross‐variogram(s). The conventional estimates of these variograms, obtained by the method‐of‐moments, are unduly affected by outlying data which inflate the variograms and so also the estimates of the error variance of cokriging predictions. Robust estimators are less affected. Robust estimators of the auto‐variogram and the pseudo cross‐variogram have previously been proposed and used successfully, but the multivariate problem of estimating the cross‐variogram robustly has not yet been tackled. Two robust estimators of the cross‐variogram are proposed. These use covariance estimators with good robustness properties. The robust estimators of the cross‐variogram proved more resistant to outliers than did the method‐of‐moments estimator when applied to simulated fields which were then contaminated. Organic carbon and water content of the soil was measured at 256 sites on a transect and the method‐of‐moments estimator, and the two robust estimators, were used to estimate the auto‐variograms and cross‐variogram from a prediction subset of 156 sites. The data on organic carbon included a few outliers. The method‐of‐moments estimator returned larger values of the auto‐ and cross‐variograms than did either robust estimator. The organic carbon content at the 100 validation sites on the transect was estimated by cokriging from the prediction data plus a set of variograms fitted to the method‐of‐moments estimates and two sets of variograms fitted to the robust estimates. The ratio of the actual squared prediction error to the cokriging estimate of the error variance was computed at each validation site. These results showed that cokriging using variograms obtained by the method‐of‐moments estimator overestimated the error variance of the predictions. By contrast, cokriging with the robustly estimated variograms gave reliable estimates of the error variance of the predictions.