Sample adequately to estimate variograms of soil properties

The variogram is central in the spatial analysis of soil, yet it is often estimated from few data, and its precision is unknown because confidence limits cannot be determined analytically from a single set of data. Approximate confidence intervals for the variogram of a soil property can be found numerically by simulating a large field of values using a plausible model and then taking many samples from it and computing the observed variogram of each sample. A sampling distribution of the variogram and its percentiles can then be obtained. When this is done for situations typical in soil and environmental surveys it seems that variograms computed on fewer than 50 data are of little value and that at least 100 data are needed. Our experiments suggest that for a normally distributed isotropic variable a variogram computed from a sample of 150 data might often be satisfactory, while one derived from 225 data will usually be reliable.     

基于最大似然法与矩法的黄土高原小流域土壤碳氮空间变异分析

以黄土高原寺底沟小流域为研究对象,根据不同土地利用方式采集46个样点的土壤样品,通过地统计方法对土壤有机碳和全氮的空间变异特征进行了分析。采用受限最大似然法(REML)和矩法(MOM)两种方法分别对变异函数进行了估计,通过交叉检验选择克里金预测效果较好的变异函数进行地统计插值。(1)与矩法(MOM)相比,在多数情况下受限最大似然法(REML)估计的变异函数进行克里金插值更加准确。(2)土层深度对土壤全氮空间变异影响较小,对土壤有机碳影响较大,表层土壤有机碳含量及变异程度明显高于下层土壤。(3)土地利用方式对土壤有机碳和全氮的空间分布有重要影响,灌木林和天然草地土壤有机碳和全氮水平最高,弃耕地其次,梯田、果园、人工草地最低,表明退耕还林对提高土壤碳氮水平有重要贡献。     

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.     

On spatial prediction of soil properties in the presence of a spatial trend: the empirical best linear unbiased predictor (E-BLUP) with REML

Geostatistical estimates of a soil property by kriging are equivalent to the best linear unbiased predictions (BLUPs). Universal kriging is BLUP with a fixed‐effect model that is some linear function of spatial co‐ordinates, or more generally a linear function of some other secondary predictor variable when it is called kriging with external drift. A problem in universal kriging is to find a spatial variance model for the random variation, since empirical variograms estimated from the data by method‐of‐moments will be affected by both the random variation and that variation represented by the fixed effects. The geostatistical model of spatial variation is a special case of the linear mixed model where our data are modelled as the additive combination of fixed effects (e.g. the unknown mean, coefficients of a trend model), random effects (the spatially dependent random variation in the geostatistical context) and independent random error (nugget variation in geostatistics). Statisticians use residual maximum likelihood (REML) to estimate variance parameters, i.e. to obtain the variogram in a geostatistical context. REML estimates are consistent (they converge in probability to the parameters that are estimated) with less bias than both maximum likelihood estimates and method‐of‐moment estimates obtained from residuals of a fitted trend. If the estimate of the random effects variance model is inserted into the BLUP we have the empirical BLUP or E‐BLUP. Despite representing the state of the art for prediction from a linear mixed model in statistics, the REML–E‐BLUP has not been widely used in soil science, and in most studies reported in the soils literature the variogram is estimated with methods that are seriously biased if the fixed‐effect structure is more complex than just an unknown constant mean (ordinary kriging). In this paper we describe the REML–E‐BLUP and illustrate the method with some data on soil water content that exhibit a pronounced spatial trend.     

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.     

Identifying management zones in agricultural fields using spatially constrained classification of soil and ancillary data

Site‐specific management requires accurate knowledge of the spatial variation in a range of soil properties within fields. This involves considerable sampling effort, which is costly. Ancillary data, such as crop yield, elevation and apparent electrical conductivity (EC_a) of the soil, can provide insight into the spatial variation of some soil properties. A multivariate classification with spatial constraint imposed by the variogram was used to classify data from two arable crop fields. The yield data comprised 5 years of crop yield, and the ancillary data 3 years of yield data, elevation and EC_a. Information on soil chemical and physical properties was provided by intensive surveys of the soil. Multivariate variograms computed from these data were used to constrain sites spatially within classes to increase their contiguity. The constrained classifications resulted in coherent classes, and those based on the ancillary data were similar to those from the soil properties. The ancillary data seemed to identify areas in the field where the soil is reasonably homogeneous. The results of targeted sampling showed that these classes could be used as a basis for management and to guide future sampling of the soil.     

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.     

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.     

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.     

县域农田土壤有机质空间变异及其影响因素分析

研究县域农田空间变异特征可以为培肥地力,增加作物产量提供指导。本文运用地统计学和 GIS相结合的方法,研究了四川省双流县土壤有机质的空间变异特征及其影响因素。结果表明： ①研究区域土壤有机质含量处于中等偏高水平,平均值为 29.72 g/kg,变异系数为 30.11%,属中等变异强度。②有机质变异函数的理论最佳模型为球状模型,块金值与基台值之比为12.67%,表明有机质含量具有强烈的空间相关性,空间相关距离为 91.10 km,普通Kriging插值表明土壤有机质含量呈现北部向东南部减少的趋势。③影响有机质空间变异的主要因素为土壤类型、地貌类型等结构性因子,而土地利用方式、施肥等随机性因子也对有机质空间变异产生重要影响,其中秸秆还田是有机质含量普遍升高的原因。     

Spatial structures of copper,lead, and mercury contents in surface soil in the Shenzhen area

Eighty-three surface soil samples were collected from the Shenzhen area for determination of copper, lead, and mercury contents. The nature of spatial dependence of the measured results was quantified using variogram analysis. All variograms show well-defined structure with zero nugget and distinct sills and ranges and can be fitted by a spherical model. The range scale and the geometric anisotropy of the variograms suggest that the spatial structures of copper and lead relate closely to the distribution of parent material in the area. The variogram of mercury appears to be isotropic with a relatively small range, indicating significant influence of geographical distribution of paddy soil fields that have been severely polluted by agricultural practice.     

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.     

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

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

红壤区土壤有机碳时间变异及合理采样点数量研究

相对于土壤有机碳(soil organic carbon,SOC)空间变异性及合理采样点数量的研究,其时间变异性及揭示特定时段SOC变化所需采样点数量的研究较少。选择红壤丘陵区的江西省余江县为研究区,分析了1982—2007年SOC含量的时间变异特征,并估算了揭示该时段SOC变化所需土壤采样点数量。结果表明,1982—2007年SOC含量均值由14.18增至16.27 g kg-1,增幅为14.74%,其变异系数则由0.22上升为0.44。各土地利用方式中,水田和林地SOC含量分别增加了2.93和3.12g kg-1,而旱地则降低了2.55 g kg-1;同时各利用方式的SOC含量变异系数均出现较大幅度的提高。基于两时段的全部样点,在95%和90%置信区间上,计算得到揭示该时段全县SOC时间变异所需的采样点数量分别为186和147。基于各土地利用方式的SOC变化,计算得到水田、旱地和林地所需采样点数量分别为68、44和144(95%置信区间)及54、34和112(90%置信区间);揭示旱地SOC变化所需采样点数量应为水田的60%以上,而林地所需样点则为水田的2倍以上。该研究结果可为红壤区SOC时间变异性及其调查采样提供参考。     

On the Akaike Information Criterion for choosing models for variograms of soil properties

A problem in the application of geostatistics to soil is to find satisfactory models for variograms of soil properties. It is usually solved by fitting plausible models to the sample variogram by weighted least squares approximation. The residual sum of squares can always be diminished, and the fit improved in that sense, by adding parameters to the model. A satisfactory compromise between goodness of fit and parsimony can be achieved by applying the Akaike Information Criterion (AIC). For a given set of data the variable part of the AIC is estimated by where n is the number of experimental points on the variogram, R is the residual sum of squares and p is the number of parameters in the model. The model to choose is the one for which  is least.
The AIC is closely related to Akaike's earlier final prediction error and the Schwarz criterion. It is also equivalent to an F test when adding parameters in nested models.     

The analysis of ranked observations of soil structure using indicator geostatistics

Structure is an important physical feature of the soil that is associated with water movement, the soil atmosphere, microorganism activity and nutrient uptake. A soil without any obvious organisation of its components is known as apedal and this state can have marked effects on several soil processes. Accurate maps of topsoil and subsoil structure are desirable for a wide range of models that aim to predict erosion, solute transport, or flow of water through the soil. Also such maps would be useful to precision farmers when deciding how to apply nutrients and pesticides in a site-specific way, and to target subsoiling and soil structure stabilization procedures.

Typically, soil structure is inferred from bulk density or penetrometer resistance measurements and more recently from soil resistivity and conductivity surveys. To measure the former is both time-consuming and costly, whereas observations made by the latter methods can be made automatically and swiftly using a vehicle-mounted penetrometer or resistivity and conductivity sensors. The results of each of these methods, however, are affected by other soil properties, in particular moisture content at the time of sampling, texture, and the presence of stones. Traditional methods of observing soil structure identify the type of ped and its degree of development. Methods of ranking such observations from good to poor for different soil textures have been developed. Indicator variograms can be computed for each category or rank of structure and these can be summed to give the sum of indicator variograms (SIV).

Observations of the topsoil and subsoil structure were made at four field sites where the soil had developed on different parent materials. The observations were ranked by four methods and indicator and the sum of indicator variograms were computed and modelled for each method of ranking. The individual indicators were then kriged with the parameters of the appropriate indicator variogram model to map the probability of encountering soil with the structure represented by that indicator. The model parameters of the SIVs for each ranking system were used with the data to krige the soil structure classes, and the results are compared with those for the individual indicators. The relations between maps of soil structure and selected wavebands from aerial photographs are examined as basis for planning surveys of soil structure.     

Optimized spatial sampling of soil for estimation of the variogram by maximum likelihood

Recent studies have attempted to optimize the configuration of sample sites for estimation of the variogram by the usual method-of-moments. This paper shows that objective functions can readily be defined for estimation by the method of maximum likelihood. In both cases an objective function can only be defined for a specified variogram so some prior knowledge about the spatial variation of the property of interest is necessary.This paper describes the principles of the method, using Spatial Simulated Annealing for optimization, and applies optimized sample designs to simulated data. For practical applications it seems that the most fruitful way of using the technique is for supplementing simple systematic designs that provide an initial estimate of the variogram.     

土壤空间变异研究中的半方差问题

简要回顾了土壤空间变异的研究。根据地质统计学理论和多年从事土壤空间变异研究的经验,对土壤空间变异研究的关键问题——半方差函数的基本假设、取样、模型选取及模型的检验进行了讨论,并对确定半方差函数模型应注意的问题提出建议。在保证取样样本容量的前提下,检查测定数据是否服从内蕴假设;注意提高每一个估算值的置信水平;尽量选择安全型模型作为半方差函数模型;对确定的半方差模型进行统计检验。由此可以求得较为客观合理的半方差模型。     

不同样点数量对土壤有机质空间变异表达的影响

以南京市六合区为研究区,通过完全随机和限制最小采样间距抽样分别设置5个样点系列,基于每个样点系列100次重复抽样的变异结构推断及空间预测误差结果,探讨了不同样点数量对土壤有机质(SOM)空间变异表达的影响。结果表明,两种抽样方式降低样点数量后推断的SOM含量的块金效应(C0/C0+C)均随样点数量减少而降低且限制最小采样间距抽样推断的C0/C0+C要低于完全随机抽样方法,说明适当的减少样点数量以便降低与SOM变异尺度不匹配的样点对变异结构推断的影响有助于提高SOM空间变异结构表达的可靠性。普通Kriging预测的SOM误差对比则表明,尽管两种抽样方式下空间预测的均方根误差(RMSE)随样点数量变化而波动,但均低于全部样点的预测误差;通过限制最小采样间距减少样点至250个时,SOM空间预测的RMSE最低,较全部样点预测误差降低了6%,因此,为了实现样点密度与SOM变异尺度相匹配,合理设置土壤采样点的间距及样点数量较单纯的增加采样点数量更为重要。     

Fluctuations in method‐of‐moments variograms caused by clustered sampling and their elimination by declustering and residual maximum likelihood estimation

Soil data accumulated in national and regional archives derive from many sources and tend to be concentrated in zones of particular interest. Experimental variograms computed from such data by the usual method of moments can appear highly erratic, and therefore models fitted to them are likely to be unreliable. We have explored two methods of avoiding the effects, one by computing declustering weights and incorporating them into the method of moments, the other using residual maximum likelihood. The methods are illustrated with data on bulk density, exchangeable magnesium, cation exchange capacity and organic carbon of 4182 samples of soil from numerous soil surveys in the whole of Australia and stored in the CSIRO's national archive. The experimental variograms of all four variables are erratic. Cell declustering produced much smoother sequences of estimates to which plausible models could be fitted with confidence. The residual maximum likelihood models closely matched those models over several hundred km. Finally values were simulated at the same sampling points from the residual maximum likelihood models, reproducing ‘spiky’ experimental variograms such as those computed from the data. The simulation showed that clustered design of sampling causes spiky artefacts. We conclude that where data are clustered experimental variograms should be computed with declustered weighting or variogram models be fitted by residual maximum likelihood.