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.     

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.     

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.     

Scale-dependent correlation between topsoil copper and cobalt concentrations in Scotland

The soil of south-east Scotland is locally deficient in copper and cobalt. Measurements from nearly 3000 fields for which the soil association is known were analysed to study the coregionalization of the two elements and to assess the influence of parent material on the metals' concentrations. The experimental auto- and cross-variograms revealed distinct local (1.5 km) and regional (20 km) scales of spatial variation. A combination of indicator variograms of the soil associations had the same spatial structures, suggesting that parent material influences the concentrations of the metals. The coregionalization between copper and cobalt was modelled as a linear combination of three spatial structures. The resulting structural correlation coefficients showed the two elements to be fairly strongly positively correlated at the regional scale. Kriging allowed determination and mapping of each spatial component; these maps were then compared with the spatial distribution of soil associations in the region. An analysis of variance was performed before and after filtering out the nugget and short-range spatial components. Classification by soil association (parent material) accounted for a large proportion of the variance at the regional scale, suggesting that the parent material contributes substantially more to the trace element content of the soil than had been thought earlier.     

Changes in variance and correlation of soil properties with scale and location: analysis using an adapted maximal overlap discrete wavelet transform

The magnitude of variation in soil properties can change from place to place, and this lack of stationarity can preclude conventional geostatistical and spectral analysis. In contrast, wavelets and their scaling functions, which take non‐zero values only over short intervals and are therefore local, enable us to handle such variation. Wavelets can be used to analyse scale‐dependence and spatial changes in the correlation of two variables where the linear model of coregionalization is inadmissible. We have adapted wavelet methods to analyse soil properties with non‐stationary variation and covariation in fairly small sets of data, such as we can expect in soil survey, and we have applied them to measurements of pH and the contents of clay and calcium carbonate on a 3‐km transect in Central England. Places on the transect where significant changes in the variance of the soil properties occur were identified. The scale‐dependence of the correlations of soil properties was investigated by calculating wavelet correlations for each spatial scale. We identified where the covariance of the properties appeared to change and then computed the wavelet correlations on each side of the change point and compared them. The correlation of topsoil and subsoil clay content was found to be uniform along the transect at one important scale, although there were significant changes in the variance. In contrast, carbonate content and pH of the topsoil were correlated only in parts of the transect.     

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.     

Robust estimation of the pseudo cross‐variogram for cokriging soil properties

The pseudo cross‐variogram can be used for cokriging two or more soil properties when few or none of the sampling locations have values recorded for all of them. The usual estimator of the pseudo cross‐variogram is susceptible to the effects of extreme data (outliers). This will lead to overestimation of the error variance of predictions obtained by cokriging. A solution to this problem is to use robust estimators of the pseudo cross‐variogram, and three such estimators are proposed in this paper. The robust estimators were demonstrated on simulated data in the presence of different numbers of outlying data drawn from different contaminating distributions. The robust estimators were less sensitive to the outliers than the non‐robust one, but they had larger variances. Outliers tend to obscure the spatial structure of the cross‐correlation of the simulated variables as described by the non‐robust estimator. The several estimators of the pseudo cross‐variogram were applied to a multitemporal data set on soil water content. Since these were obtained non‐destructively, direct measurements of temporal change can be made. A prediction subset of the data was subsampled as if obtained by destructive analysis and the remainder used for validation. Estimators of the auto‐variogram and pseudo cross‐variogram were applied to the prediction data, then used to predict the change in water content at the validation sites by cokriging. The estimation variances of these predictions were best calculated with a robustly estimated model of coregionalization, although the validation set was too small to conclude that the non‐robust estimators were unsuitable in this instance.     

苏北海涂围垦区土壤水分空间变异性及其协同克立格估值

以苏北滨海滩涂围垦区为研究区域,应用地统计学方法研究了剖面不同深度土壤含水量的空间变异特征,以0～20 cm表土层含水量作为协同变量,分析了不同深度土壤含水量的交互半方差特征,探讨了土壤含水量的普通克立格和协同克立格估值方法.结果表明:各层土壤含水量都表现为弱变异强度;受结构性因素和随机性因素的共同作用,各土层含水量均呈中等空间自相关性;协同区域化变量间均表现为正相关,且协同区域化变量的空间结构性要优于单一变量;协同克立格对土壤水分局部变异的描述较普通克立格更为详细;与普通克立格法相比,协同克立格法估值产生的均方误差减小10.1%～21.3%,平均标准误差减小11.3%～17.7%,预测值和实测值间的相关系数提高15.9%～26.4%.该研究为滨海滩涂地区土壤水分分区管理和水盐优化调控提供了一定的理论参考.     

Exploring scale‐dependent correlation of soil properties by nested sampling

The value of nested sampling for exploring the spatial structure of univariate variation of the soil has been demonstrated in several studies and applied to practical problems. This paper shows how the method can be extended to the multivariate case. While the extension is simple in theory, in practice the direct estimation of covariance components by equating mean‐square matrices with their expectation will often lead to estimates that are not positive semidefinite. This paper discusses solutions to this problem for balanced and unbalanced sample designs. In the balanced case there is a residual maximum likelihood (REML) estimator that will find estimates of covariance components that maximize an overall likelihood on the condition that all components are positive semidefinite (p.s.d.). This is possible because the condition is met if the differences of successive mean‐square matrices are positive semidefinite, and this constraint can be incorporated into an algorithm. This does not hold for unbalanced designs. In this paper the problem was solved for unbalanced designs by scaling covariance components that were not p.s.d. to the nearest p.s.d. matrix according to a Euclidean distance. These methods were applied to data from three surveys, two with balanced and one with unbalanced sampling. Different patterns of scale‐dependence of the correlation of soil properties were found. For example, at Ginninderra Experimental Station in Australia the soil water content and bulk density were correlated significantly, with the correlation increasing with distance to 56 m, but at longer distances the properties were not significantly correlated. By contrast, the pH of the soil and the available P content showed correlation that increased with distance. The implications of these results for planning more detailed sampling, both for prediction and for investigation of processes, are discussed.     

Multiscale spatial variability of CO2 emissions and correlations with physico-chemical soil properties

Spatial variability of greenhouse gas (GHG) emissions from agricultural lands is not well known although it has a great impact on the accuracy of GHG budget.The objectives of this study were to assess the spatial variability of CO₂ emission fluxes (CO_2-flux) and correlate these emissions with soil physico-chemical properties at two spatial scales and at different depths using a new geostatistical approach (coregionalization analysis with a drift, CRAD) that performs multiscale spatial analysis.Two agricultural sites with sandy and loamy soils were instrumented at 108 geo-referred sampling points and at two depths during spring 2007 where soil surface CO_2-flux and soil physico-chemical parameters were measured. The CO_2-flux presented spatial patterns characterized by different scales (i.e., non-spatial, small spatial and large spatial scale components), each describing a different fraction of its variability. About a quarter of CO_2-flux variability at the first site and one fifth at the other site was attributed to the non-spatial component. Strongest correlations were obtained between CO_2-flux and soil temperature, water saturation (S_w), elevation, electrical conductivity, soil bulk density, and the C/N ratio, but with differences between sites. Correlations were much stronger at large scale. Analyzing correlations between CO_2-flux and soil properties without discriminating for scales can miss important scale-dependent processes controlling soil gas emissions. Scales at which these processes vary should therefore be taken into account.     

Spatial Variability in Slash Linear Modeling with Finite Second Moment

This article studies the dependence of spatial linear models using a slash distribution with a finite second moment. The parameters of the model are estimated with maximum likelihood by using the EM algorithm. To avoid identifiability problems, the cross-validation, the Trace and the maximum log-likelihood value are used to choose the parameter for adjusting the kurtosis of the slash distribution and the selection of the model to explain the spatial dependence. We present diagnostic techniques of global and local influences for exploring the sensibility of estimators and the presence of possible influential observations. A simulation study is developed to determine the performance of the methodology. The results showed the effectiveness of the choice criteria of the parameter for adjusting the kurtosis and for the selection of the spatial dependence model. It has also showed that the slash distribution provides an increased robustness to the presence of influential observations. As an illustration, the proposed model and its diagnostics are used to analyze an aquifer data. The spatial prediction with and without the influential observations were compared. The results show that the contours of the interpolation maps and prediction standard error maps showed low changes when we removed the influential observations. Thus, this model is a robust alternative in the spatial linear modeling for dependent random variables. Supplementary materials accompanying this paper appear online.     

Exploring the variation in soil saturated hydraulic conductivity under a tropical rainforest using the wavelet transform

Saturated hydraulic conductivity (K_s) of the soil is a key variable in the water cycle. For the humid tropics, information about spatial scales of K_s and their relation to soil types deduced from soil map units is of interest, as soil maps are often the only available data source for modelling. We examined the influence of soil map units on the mean and variation in K_s along a transect in a tropical rainforest using undisturbed soil cores at 0–6 and 6–12 cm depth. The K_s means were estimated with a linear mixed model fitted by residual maximum likelihood (REML), and the spatial variation in K_s was investigated with the maximum overlap discrete wavelet packet transform (MODWPT). The mean values of K_s did not differ between soil map units. The best wavelet packet basis for K_s at 0–6 cm showed stationarity at high frequencies, suggesting uniform small‐scale influences such as bioturbation. There were substantial contributions to wavelet packet variance over the range of spatial frequencies and a pronounced low frequency peak corresponding approximately to the scale of soil map units. However, in the relevant frequency intervals no significant changes in wavelet packet variance were detected. We conclude that near‐surface K_s is not dominated by static, soil‐inherent properties for the examined range of soils. Several indicators from the wavelet packet analysis hint at the more dominant dynamic influence of biotic processes, which should be kept in mind when modelling soil hydraulic properties on the basis of soil maps.     

区域尺度下黄土高原土壤全钾含量的空间模拟

为了揭示高原区域地区尺度上土壤全钾的空间异质性及其影响因素,该文采用状态空间方程和传统线性回归模型对该区土壤全钾含量的空间分布进行了模拟,并分析了其与土壤体积质量、黏粒含量、粉粒含量、土壤酸度、降水、气温和海拔高度等因素之间的关系。结果表明,以上变量在30～50km的采样间距下均表现出较好的空间自相关性,其中土壤体积质量、黏粒含量、粉粒含量、降水和气温与土壤全钾之间存在显著的交互相关关系,可用于土壤全钾的状态空间模拟。不同因素组合下的状态空间方程均比使用相同变量的线性回归方程能更好的模拟土壤全钾含量的空间分布。使用土壤体积质量和黏粒含量的双因素状态空间方程模拟效果最好,决定系数R2为0.978,均方根误差(RMSE)为0.049。状态空间模拟在大尺度区域的应用表现出较好的效果,为研究该区其他土壤属性的空间异质性提供了参考。     

Choosing functions for semi-variograms of soil properties and fitting them to sampling estimates

The semi-variogram is central to geostatistics and the single most important tool in geo-statistical applications to soil. Mathematical functions for semi-variograms must be conditional negative semi-definite, and there are only a few families of simple function that meet this demand. These include the transitive models with finite a priori variance deriving from moving average processes. The spherical and exponential schemes are the most often encountered members. The other major group is that of unbounded models in which the variance appears to increase without limit. The linear model is the most common in this group. If more complex models are needed they can be formed by combining two or more simple models. The usual estimator of the semi-variance is often considered inefficient and to be sensitive to departures from normality in the data. It is compared with a robust estimator and shown to be generally preferable in being unbiased and having confidence intervals that are no wider. For routine analysis, fitting models to sample semi-variograms by weighted least squares approximation, with weights proportional to the expected semi-variance, is preferred to the more elaborate and computationally demanding statistical procedures of generalized least squares and maximum likelihood. The Akaike information criterion is recommended for selecting the best model from several plausible ones to describe the observed variation in soil, though for kriging it may be desirable to validate the chosen model. Examples of models fitted to soil semi-variograms are shown and compared.     

复杂景观区土壤有机质预测模型的尺度效应

为研究复杂景观区土壤有机质预测模型的尺度效应,探讨不同空间尺度数据综合利用的问题,本文运用回归Kriging方法对河南省登封市土壤有机质进行预测,分析了不同空间尺度数据在建模过程中的作用和影响.结果表明:土壤有机质关于高程因子的趋势属于宏观趋势,以大尺度数据拟合该趋势值效果最优;小尺度数据不适合用于拟合土壤有机质的趋势值,但揭示了小尺度残差值的空间变异细节,增强了大尺度残差值的空间结构性,能够有效提高土壤有机质预测精度.因此,景观复杂区土壤有机质预测中,应基于大尺度数据模拟趋势值,大尺度数据和小尺度数据相结合拟合残差值的空间变异函数以预测残差值,最后趋势值加上残差值得到土壤有机质预测值.     

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.     

Estimation of soil salinity in semi‐arid land using a geostatistical model

Soil salinity is one of the great problems in arid and semi‐arid environments. The estimation and prediction of spatial soil salinity may be considered as a stochastic process, observed at irregular locations in space. Environmental variables usually show spatial dependence among observations which is an important drawback to traditional statistical methods. Geostatistical techniques that analyse and describe the spatial dependence and quantify the scale and intensity of the spatial variation, provides spatial information for local estimation of soil salinity. In this paper we propose a Gaussian Spatial Linear Mixed Model (GSLMM), which involves a non‐parametric term, accounting for a deterministic trend given by exogenous variables, and a parametric component defining the purely spatial random variation, possibly due to latent spatial processes. We focus here on the analysis of the relationship between soil electrical conductivity as a parameter related directly to soil salinity as well as sodium content (sodicity) to identify spatial variations in these parameters. This kind of methodology is demonstrated as a useful tool for environmental land management. Copyright © 2007 John Wiley & Sons, Ltd.     

Estimating temporal change in soil monitoring: I. Statistical theory

Detecting small temporal change of spatially varying soil properties demands precise estimation. Design– and model–based methods are compared for estimating temporal change of soil properties over finite areas. Analytical expressions for the estimators and their variances arc derived for the two approaches, and formulae for the expectations of the variances under the random–process model are developed. Among the randomized designs simple, stratified, and systematic random sampling using the arithmetic mean as estimator have been studied. Pairing the sampling positions on the different occasions increases the precision of design–based estimation if the observations are positively cross–correlated. The relative precisions of the means of stratified and systematic samples depends on the spatial correlation. Neither is more precise than the other in all circumstances. The stratified design provides an unbiased estimator for the sampling error, which is not available from systematic samples. Theoretically, the geostatistical global estimator is more precise than the estimates derived from any of the classical designs when many realizations arc repeatedly sampled at random. In practice, with only a single realization of the process, this is no longer relevant. Moreover, errors in estimating the variograms add to the total error of the method. It seems that only by sampling from large auto–correlated random fields can the precisions of the methods be compared in practice.     

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.