Spatial estimation of surface soil texture using remote sensing data

Abstract

Understanding the spatial distribution and variability of soil texture is essential for land use planning and other activities related to agricultural management and environmental protection. This study was conducted to evaluate Landsat Enhanced Thematic Mapper (ETM) remote sensing data as auxiliary variables for spatial estimation of surface soil texture using a limited number of soil samples taken from a site located in the city of PingduShandong ProvinceChina. Three methods of evaluating variability in surface soil texture were evaluated: (1) multiple stepwise regression (MSR) based on the relationship between surface soil sandsilt and clay contents and remote sensing data; (2) kriging of surface soil sandsilt and clay contents; (3) cokriging with remote sensing data. Correlation analysis showed that surface soil sandsilt and clay contents were significantly correlated with Landsat ETM digital number (DN) of six bands (Bands 1–5 and Band 7)and the DN of Band 7 explained most of the variability in soil sandsilt and clay contents. The DN of Band 7 was selected as auxiliary data for the estimation of surface soil texture. The cross-validation results indicated that both MSR and kriged estimates had low reliability due to the variations in landscape and the low-density sampling in the study area. Cokriging with remote sensing data significantly improves estimates of surface soil texture compared with MSR and kriging.

Key words:

• soil texture
• remote sensing
• spatial estimation
• multiple stepwise regression
• geostatistics

INTRODUCTION

Soil texture is considered an important influential factor on many of the physical and chemical properties and behavior of soilsuch as water storagecation exchange capacitysoil fertilityinternal drainage and sorption characteristics (Rawls et al. 1982; Manrique et al. 1991; Gawlik et al. 1999; Makabe et al. 2009). In terms of organic matter contentsoil texture is often used for developing pedotransfer functions to predict soil hydraulic properties and cation exchange capacity (Wösten et al. 2001; Seybold et al. 2005). Furthermoresoil texture is also a pertinent property to quantify vulnerability to soil erosion (Le Bissonnais 1996; Warrington et al. 2009). Thereforehigh-resolution maps of soil texture are essential for hydrological and ecological modelingas well as other activities related to agricultural management and environmental protection (Hassink 1992; Oberthür et al. 1996; Zhao et al. 2009).

Over the last two decadesboth remote sensing images of bare soil and spectral reflectance of soil samples have been successfully used for accurate and rapid estimation of soil properties (Henderson et al. 1989; Agbu et al. 1990; Sudduth and Hummel 1991; Ben-Dor and Banin 1995; Ben-Dor et al. 1997 1999; Chen et al. 2000; McCarty et al. 2002; Brown et al. 2006; Ben-Dor et al. 2008; Chen et al. 2008; D'Acqui et al. 2010). For exampleChen et al. (2000) used remotely sensed imagery to evaluate surface soil organic matter at landscape scale. Ben-Dor et al. (2008) combined active (using imaging spectroscopy) and passive (using ground penetrating radar and frequent domain electromagnetic) remote sensing methods for assessing soil salinity in the field. Studies also have indicated that soil texture correlates significantly with soil reflectance in the visible and near-infrared region (400–2500 nm) (Al-Abbas et al. 1972; Suliman and Post 1988; Zhang et al. 1992; Sullivan et al. 2005). Thereforesoil reflectance measurements can be used for estimating soil texture by statistical methods such as stepwise linear regression and partial least squares. The relationship between reflectance and soil texturehoweverdegrades when going from the laboratory scale to the field and large scalesdue to the variations in surface status including roughnesssoil moisture and nature of pebbles (Lagacherie et al. 2008). Although soil texture can be estimated using remote sensing data such as Landsat Enhanced Thematic Mapper (ETM) imagerythis approach may yield unacceptable performance over large areas.

Geostatistics is a powerful tool for characterizing and mapping the spatial distribution and variability of soil properties (Burgess and Webster 1980). Kriging is a basic geostatistical technique that provides the best linear unbiased estimation (BLUE) for a spatially dependent variable. Cokriging is another geostatistical method that extends kriging of a primary variable to secondary variables based on their correlations with the primary variable. It has been demonstrated that cokriging is superior to kriging for minimizing estimation variance (Vauclin et al. 1983; Istok et al. 1993; Wu et al. 2009). If auxiliary variables are not highly correlated to the primary variablecokriging is only slightly superior to (or at least equal to) ordinary kriging (Shouse et al. 1990; Triantafilis et al. 2001). Yates and Warrick (1987) suggested that cokriging is superior to ordinary kriging for predictions for 0.5 or larger correlation coefficients between variables and when the auxiliary variables are oversampled.

The spatial distribution of soil texture can be predicted by kriging interpolation. This method will produce satisfying results if enough sample points are available. Regretfullylaboratory measurements (e.g.pipette method) of the particle size distribution (PSD) are costly and time-consuming. Cokriging can also be used for predicting soil texture when it is highly correlated to easily measured auxiliary variables. Lagacherie et al. (2008) found that soil texture correlates with soil reflectance significantly over small to large scales with low variation in surface roughness. Compared with the measurement of the PSDremote sensing data is much easier to obtain for detailed descriptions of large areas. Thusremote sensing data such as Landsat ETM imagery may be useful auxiliary variables for the prediction of surface soil texture. While there have been several studies focusing on the spatial estimation of soil texture using laboratory soil reflectance data (Zhang et al. 1992) and satellite remote sensing data such as IKONOS (Sullivan et al. 2005) and AVHRR imageries (Odeh and McBratney 2000)to our knowledge little research has been conducted using Landsat ETM imagery. Thereforethe objectives of this study are to: (1) investigate the utilization of Landsat ETM remote sensing data as auxiliary variables for the estimation of surface soil textureand (2) evaluate the performance of three methods (multivariate regressionkriging and cokriging) in assessing the variability of surface soil texture.

MATERIALS AND METHODS

Descriptions of study area

This study was conducted in agricultural regions of Pingdu City. The study area (36°28′–37°02′ N119°31′–120°19′ E)covering 3166.54 km²is located in the eastern part of Shandong ProvinceChina (Fig. 1). Pingdu is located in a warm temperate climate zone with marine monsoonshot rainy summersand cold dry winters. The annual mean temperature is 12.3°C with an average annual sunshine of 2541 h. The study area is a typical dry-land farming area with average annual precipitation of 707 mmover 70% of which occurs from July to September. The landscape of the area is generally sloping. The elevation ranges from more than 150 m above sea level in the northeast to less than 50 m in the south and center of the study area (Fig. 2). Alsothe typical crop rotation for this region is winter wheat (Triticum aestivum L.) and summer maize (Zea mays L.) cropping.

Figure 1 The soil types, soil sampling sites (58), and sampling locations (344; 3-km grid) on a Landsat Enhanced Thematic Mapper (ETM) image in the study area.

Figure 2 The elevation change within the study area.

Parent materials in this area are mainly residuum and colluvium granite and gneissalluvialand lacustrine deposits. Soils formed in these parent materials are classified as LuvisolsFluvisols and Vertisols respectivelyaccording to the World Reference Base for Soil Resources (WRB) (Food and Agriculture Organization of the United Nations 1998). These three soil types cover more than 96% of the area. The Luvisols are mainly distributed in the northeastern part of the study area with high elevation and steeper slopewhile the Fluvisols are mainly distributed around rivers with low elevation and gentle slopeand Vertisols are also mainly distributed in the areas with low elevation and gentle slopebut far away from rivers (Fig. 1 and 2). Additionallyvery few Lixisols and Solonchaks are also distributed in this region.

Soil sampling and analysis

Based on the number of samples being proportional to each soil typea total of 58 samples were randomly taken from topsoil (0–15 cm) in November 2007 in agricultural regions of Pingdu Citycovering 3166.54 km² (Fig. 1). Twenty-seven samples of Luvisols16 samples of Fluvisols and 15 samples of Vertisols were collected. When samplingfive surface soil cores were collected with a hand auger at each sampling site. These five soil cores were fully mixedand then approximately 1 kg of the mixed soil was packed in a bag and brought back to the laboratory for analysis. The coordinates of all sampling locations were recorded using a global positioning system (GPS) receiver. Soil samples were air-driedpassed through a 2-mm sieve and then analyzed in the laboratory for the contents of sand (0.05–2 mm)silt (0.002–0.05 mm) and clay (0–0.002 mm) particles in the soil following the pipette method (Soil Science Society of China 2000).

To predict soil texture in the study area using remote sensing data as auxiliary variableswe selected a Landsat ETM image acquired on 7 November 2007. Landsat data are available for free download from the USGS Global Visualization Viewer website (U.S. Geological Survey 2010). Most values of the normalized difference vegetation index (NDVI)which is a commonly used vegetative index in remote sensing analysiswere smaller than 0.1 in the study area. This indicated that the agricultural fields were mostly bare at the time of Landsat data acquisition. To use remote sensing data as an auxiliary variable for the prediction of surface soil texturewe sampled 402 pixels; 344 pixels were sampled using a grid-based sampling scheme with 3-km grid spacingand the other 58 pixels were sampled based on the corresponding locations of the 58 soil samples described above (Fig. 1).

Multivariate regression

Multiple stepwise regression (MSR) was conducted for surface soil texture of the 58 soil samples (dependent variable) with digital numbers (DNs) of six bands (Bands 1–5 and Band 7) of the Landsat ETM image (independent variables). A backward method regression (Norusis 1994) was selected and the level for entry in the regression model was set at p < 0.10while a 0.05 significance level was applied to retain the variables in the model. The DNs of six bands were taken as independent variables and the contents of sandsilt and clay particles were taken as dependent variables in the linear regression model.

The multiple linear regression model was fitted using the regress function of MATLAB software (The MathWorks Inc.USA).

Kriging and cokriging

Ordinary kriging was used in this study. It is a typical geostatistical approach that relies on the observed soil textures and their corresponding spatial positions to predict the soil textures at unsampled locations. Cokriging is an extension of the kriging method that integrates the information carried by a secondary variable related to the primary variable being estimated. Such spatial variables are also called co-regionalized and are spatially dependent. In this studycokriging is expected to produce a more accurate estimation since the relationship between surface soil texture and remote sensing data (e.g.DN of Band 7) was considered in the interpolation. Both kriging and cokriging are viewed as unbiased prediction methods with minimum prediction variance. They respectively use the semivariogram and cross-semivariogram to quantify the spatial patterns of the dependent variableboth isotropically and anisotropically (Western et al. 2004). The estimator for the semivariogram and cross-semivariogram is:

(1)

where γ _uv is the semivariance (when u = v) with respect to random variable z_u h is the separation distanceand N(h) is the number of pairs of z_u (x_i ) and z_v (x_i ) in a given lagged distance interval of (h + dh). When u ≠ vγ _uv is the cross-semivariogramwhich is a function of h (Yates and Warrick 1987). In this studyanisotropy of variograms was not founded. Four variogram models (sphericalexponentiallinear and Gaussian) were used to describe the semivariograms or cross-semivariograms and the best-fitted models were selected with the largest coefficient of determination (R ²) between model predicted variances and the measured values of soil texture (Wang et al. 2002).

The semivariogram and cross-semivariogram can be used for kriging and cokriging analyses. The general equation of kriging estimator is:

(2)

where z*(x_p ) is the kriging estimate at location x_p z(x_i ) is the known value at location x_i and λ _i is the weight of sample z at x_i . In order to achieve unbiased estimations in krigingthe following set of equations should be solved simultaneously:

(3)

where γ(x_i x_j ) is the value of the variogram corresponding to a vector with origin in x_i and extremity in x_j and μ is the Lagrange coefficient. Alsothe general form of the cokriging equation can be formulated as:

(4)

subject to:

(5)

in which is the cokriging estimate at location x_p and a_i and b_j are weights of z (primary variable) at x_i and y (secondary variable) at x_j respectively. The weights are calculated by solving the cokriging equations which can be found in geostatistical textbooks (Journel and Huijbregts 1978; Isaaks and Srivastava 1989).

All the geostatistical computations were performed using GS+ version 9.0 (Gamma Design Software LLC.PlainwellMI). ArcGIS 9.3 (ESRIThe RedlandsCA) was used for mapping and spatial analysis.

Evaluation criteria

Due to the relatively small number of sites used in this studya cross-validation technique was used to evaluate the predictive power of kriging and cokriging (Lagacherie et al. 2008). We removed one sample from the data setadjusted the prediction model (e.g. kriging or cokriging) over the 57 remaining samplesand predicted the ‘leave out’ individual. We repeated this process until all samples had been removed individually. Estimated and measured soil texture were used to calculate the mean error (ME) and root mean square error (RMSE) which are defined based on Isaaks and Srivastava (1989):

(6)

(7)

where z(u_i) is the measured value of z at location u_i and z*(u_i) is the predicted value at the same location. The ME provides a measure of biasand the RMSE provides a measure of accuracy.

RESULTS

Correlation analysis

In order to detect the relationships between surface soil texture and Landsat ETM DNs of six bands (Bands 1–5 and Band 7)the linear correlation analysis was appliedand the results are shown in Table 1. A positive correlation exists between soil sand content and DNs of six bandsand silt and clay content have negative correlation with DNs of six bands. The DN of Band 7 has the strongest absolute correlation with surface soil texture.

Table 1  Pearson correlations between surface soil sand, silt and clay contents and digital numbers (DNs) of Landsat Enhanced Thematic Mapper (ETM) imagery

	DN 1	DN 2	DN 3	DN 4	DN 5	DN 7
Sand	0.223	0.297^a	0.436^b	0.489^b	0.503^b	0.568^b
Silt	−0.173	−0.268^a	−0.398^b	−0.412^b	−0.422^b	−0.453^b
Clay	−0.264^a	−0.366^b	−0.476^b	−0.522^b	−0.545^b	−0.597^b
^aCorrelation is significant at the 0.05 level^bCorrelation is significant at the 0.01 level. DN 1DN 2…DN 7 = The digital number of Band 1Band 2…Band 7 of Landsat ETM imagery.

Multiple linear regression analysis

The derivation of regression models for surface soil texture was performed for the 58 soil samples through multiple stepwise regression (MSR) using the DN values of six bands. Results indicated that only DN7 remained as a significant predict variable (p <0.05):

(8)

(9)

(10)

The coefficients of determination (R ²) of the model are 0.320.21 and 0.36 for sandsilt and clay respectivelywhich means that the regressions explain 32%21% and 36% of the variances in sandsilt and clayrespectively.

Geostatistical analysis

The quantile-quantile plots of surface soil texturewhich can be used to check whether this data set follows distributionare listed in Fig. 3. This figure shows most of the points align themselves along the 45-degree linethus indicating that surface soil texture meets the requirement of a normal distribution for geostatistical analysis. We derived the semivariance and cross-semivariance values of soil texture and DN of Band 7 and selected the best geostatistical model that fitted each variogram from the highest R ² values (Table 2).

Table 2  Semivariograms and cross-semivariograms for surface soil sand, silt and clay contents and digital number (DN) of Band 7 of Landsat Enhanced Thematic Mapper (ETM) imagery

	Model	Nugget	Sill	Nugget/Sill (%)	Range (km)	R ^2
Semivariograms
Sand	Spherical	4.70	125.20	3.75	10.17	0.28
Silt	Gaussian	31.40	71.25	44.07	12.51	0.45
Clay	Exponential	2.15	32.03	6.71	9.99	0.21
DN 7	Spherical	0.10	302.10	0.03	8.27	0.44
Cross-semivariograms
Sand-DN 7	Gaussian	0.10	101.00	0.10	10.83	0.44
Silt-DN 7	Spherical	−0.10	−51.00	0.20	22.04	0.48
Clay-DN 7	Spherical	−0.10	−53.07	0.19	10.17	0.23
DN 7digital number of Band 7 of Landsat ETM imagery.

Figure 3 The quantile-quantile plots of (a) soil sand, (b) silt and (c) clay content.

The results of the geostatistical analysis are presented in Table 2. The semivariograms of sandsilt and clay are fitted well by sphericalGaussian and exponential models respectivelyand a spherical model fits the semivariogram of DN of Band 7 well. The cross-semivariogram of sand is fitted well by a Gaussian modeland a spherical model fits the cross-semivariograms of silt and clay well. Because of the negative linear correlation between silt (r = –0.453) and clay and DN of Band 7 (r = –0.597) (Table 1)negative nugget and sill values for silt-DN7 and clay-DN 7 are listed in Table 2. The nugget/sill ratios of the kriging semivariogram of soil sandsilt and clay contents are 3.75%44.07% and 6.71%respectivelywhile the nugget/sill ratios of the cokriging cross-semivariogram of sandsilt and clay contents are 0.10%0.20% and 0.19%respectively. The range values of the kriging semivariogram of soil sandsilt and clay contents are 10.1712.51 and 9.99 kmrespectivelyand the range values of the cokriging cross-semivariogram of sandsilt and clay contents are 10.8322.04 and 10.17 kmrespectively.

The above features can be clearly observed from the semivariogram and cross-semivariogram plots of the soil texture (Fig. 4). Based on semivariogram and cross-semivariogram analysisspatial patterns of predicted soil texture can be derived by both kriging and cokriging with remote sensing data (Fig. 5).

Figure 4 (a–c) Semivariograms and (d–f) cross-semivariograms of soil sand, silt and clay content.

Figure 5 Predicted soil sand, silt and clay content by (a–;c) kriging and (d–;f) cokriging.

The performance of kriging and cokriging

The cross-validation results of kriging and cokriging are presented in Table 3. For comparisona statistical description of the measured sandsilt and clay contents is also shown in this table. The minimum and maximum values of the predictions by kriging are 38.2% and 58.3% for sand26.6% and 42.3% for silt and 12.8% and 19.9% for clayrespectively. The minimum and maximum values of the predictions by cokriging with DN7 are 34.0% and 65.9% for sand23.7% and 44.2% for silt and 11.5% and 23.0% for clayrespectivelywhereas the minimum and maximum values of the observations are 23.3% and 70.8% for sand12.0% and 50.7% for silt and 7.6% and 31.2% for clayrespectively. The mean and standard deviations (SDs) of the predictions by kriging are 49.0% and 4.9% for sand33.8% and 3.9% for silt and 17.3% and 1.8% for clayrespectivelyand the mean and SDs of the predictions by cokriging with DN7 are 48.8% and 7.1% for sand33.8% and 5.0% for silt and 17.2% and 2.9% for clayrespectively. By contrastthe mean and standard deviations (SDs) of the observations are 49.0% and 11.0% for sand33.9% and 8.4% for silt and 17.1% and 5.6% for clayrespectively.

Table 3  Results of cross-validation of kriging and cokriging methods for surface soil sandsilt and clay contents

	Minimum	Maximum	Mean	SD	ME	RMSE	R ^2
Sand (%)
Measured	23.3	70.8	49.0	11.0
Kriging	38.2	58.3	49.0	4.9	0.014	10.65	0.12
Cokriging	34.0	65.9	48.8	7.1	0.159	6.90	0.62
Silt (%)
Measured	12.0	50.7	33.9	8.4
Kriging	26.6	42.3	33.8	3.9	0.099	8.23	0.09
Cokriging	23.7	44.2	33.8	5.0	0.048	5.47	0.60
Clay (%)
Measured	7.6	31.2	17.1	5.6
Kriging	12.8	19.9	17.3	1.8	−0.126	5.55	0.06
Cokriging	11.5	23.0	17.2	2.9	−0.068	4.74	0.28
SDstandard deviation; MEmean error; RMSEroot mean square error.

The MEsRMSEs and R ² of kriging are 0.014%10.65% and 0.12 for sand0.099%8.23% and 0.09 for silt and –0.126%5.55% and 0.06 for clayrespectivelywhile the MEsRMSEs and R ² of cokriging with DN7 are 0.159%6.90% and 0.62 for sand0.048%5.47% and 0.60 for silt and –0.068%4.74% and 0.28 for clayrespectively. The predicted soil sandsilt and clay contents by cokriging are more similar to observations than those predicted by kriging (Fig. 6).

Figure 6 Comparison of (a–c) kriging estimates to (d–f) cokriging estimates of soil sandsilt and clay contents.

DISCUSSION

Relationship between remote sensing data and soil texture

Soil texture has a stronger absolute correlation with near-infrared bands (Bands 4–5 and Band 7) than with visible bands (Bands 1–3) of ETM imagery due to the spectral signatures of the texture typically dominating the near-infrared spectra (Stenberg et al. 2010). Among the three textural classessand and clay contents have better correlation with the DN values than silt content. This is consistent with previous studies (Islam et al. 2003; Wetterlind and Stenberg 2010). Wetterlind and Stenberg (2010) suggested that silt was more poorly predicted due to the fact that the silt fraction can be assumed to be a mixture of sand and clay minerals and was therefore difficult to distinguish in the visible and the near-infrared regions. In this studythe DN of Band 7 has a significant correlation with soil texture in a large area with different soil types. The main reason for this result may be that the Band 7 is referred to as short wave infrared (SWIR)which is sensitive to soil available water. The available water content depends greatly on the soil texture. In additionsoil texture also has a direct impact on SWIR reflectanceas incoming radiation is scattered differently by coarse particles as compared to fine particles. For exampleseveral studies have found that the spectral absorption band center of clay particles is around 2200 nm (in the range of wavelengths of Band 7) (Clark 1999; Brown et al. 2006; Lagacherie et al. 2008).

Performance of regression

The regression models for predicting surface soil texture in our study are worse than those reported by Coleman et al. (1991)Demattê and Garcia (1999)Broge et al. (2004)and Lagacherie et al. (2008)but slightly better than those of Sullivan et al. (2005) even though we only have a total of 58 soil samples. The main reason for the less robust performance of regression models in our study is probably that using ETM imagery to obtain remote sensing data cannot completely avoid the influences of atmospherictopographic and solar effects. In additionthe spatial variability of surface soil texture within a 30-m ETM pixel may also affect the accuracy of regression models. Thereforethe remote sensing data used in this study is less effective than the spectral reflectance obtained using field multiband radiometer (Coleman et al. 1991)Infra Red Intelligent Spectroradiometer (Demattê and Garcia 1999)hyper-spectrometer (Broge et al. 2004) and ASD pro FR Portable Spectro-radiometer (Lagacherie et al. 2008)which enable direct measurement of pure soil spectra. This suggests that the regression models have to be used with caution when extrapolating the relationship between soil texture and remote sensing data over large geographic regions.

Geostatistical characteristics of soil texture

The nugget values of all variables are greater than 0indicating that the semivariograms and cross-semivariograms of soil texture and DN of Band 7 exhibit nugget effects caused either by the measurement error or by the unresolved spatial variability of the property. The nugget/sill ratio is assumed to be a criterion to classify the spatial dependence of soil properties. Ratio values lower than 25% and higher than 75% corresponded to strong and weak spatial dependencerespectivelywhile the ratio values between 25 and 75% corresponded to moderate spatial dependence (Cambardella et al. 1994; Zhu and Lin 2010). The nugget/sill ratios of the kriging semivariograms of sand and clay are lower than 25%and the nugget/sill ratios of the kriging semivariogram of silt are between 25% and 75%. This indicates that sand and clay contents have a strong spatial dependenceand silt content has a moderate spatial dependence in the area. The estimates of nugget variance are unstable for soil texturehoweverbecause their semivariances at the smallest separation distance are almost the same as the ones at the second smallest separation distance (Fig. 4). This suggests that the large variability of texture is within the smallest sampling interval. This large variabilityhoweveris well described by the cokriging cross-semivariograms of soil texturewhere the absolute semivariances at the smallest separation distance are significantly less than those at the second smallest separation distance. Overall, low nugget/sill ratios of sandsilt and clay mean that the spatial variability of soil texture is mainly controlled by intrinsic factors such as soil parent materials (Cambardella et al. 1994).

The range is the separation distance beyond which no spatial correlation exists and estimates of covariance values remain essentially constant. Range can also reflect some information about spatial dependency of soil properties (Lόpez-Granados et al. 2002). A large range of the cokriging cross-semivariogram of soil silt content indicates that observed values of this soil variable are influenced by other values of this variable over greater distances than soil variables which have smaller ranges (Lόpez-Granados et al. 2002).

Performance of kriging and cokriging

The prediction maps of surface soil texture by kriging and cokriging exhibit similar distribution trends. Howeverthe prediction map of soil texture by cokriging has more classes than the prediction map by kriging (Fig. 5)indicating that the prediction maps by cokriging provide more spatially detailed information than those by kriging. Soil sand content is relatively higher in the northeastern parts of the study area. By incorporating the DN of Band 7 in the cokrigingmore spatially detailed DN7 (a spatial resolution of 30 m) information has been considered in the interpolation of soil texture. Thereforesoil texture maps generated from cokriging can better show the detailed spatial variation of soil texture that was determined by the spatial variation of soil parent materials than maps generated by kriging that only relied on the 58 soil samples. Sand contents in Luvisols formed in coarse-grained residuum and colluviums in the mountain area are higher than those in Vertisols formed in fine-textured lacustrine deposits in the flat and low elevation area. Thereforesoil sand content in the northeastern parts of the study area where most Luvisols distributes is much higher than that in the area with Vertisols.

Variations of predicted sandsilt and clay contents by kriging and cokriging are smaller than those of measured sandsilt and clay contentsindicating a sort of smoothing effect caused by these two interpolation methods (Wang et al. 2002). Howevercokriging produced a lower minimum output and greater maximum output than did kriging. Means of predicted sandsilt and clay contents by the two methods are comparable with those of measured sandsilt and clay contentsand the SDs of predicted sandsilt and clay contents are significantly less than those of the measured values. Howevercompared with krigingSDs of predicted sandsilt and clay contents by cokriging are more close to those of measured values. Statistical descriptions of kriging and cokriging interpolation results indicate that estimation by cokriging with remote sensing data provides better predictions than estimation by kriging.

The MEs of kriging and cokriging for sand and silt are greater than 0which suggests that the two methods generally underestimate the soil sand and silt contents. The MEs of both kriging and cokriging for clay are less than 0which virtually indicates overestimation. The RMSE provides a measure of accuracy. The RMSEs of kriging for sandclay and silt contents are greater than those of cokrigingindicating that cokriging with remote sensing data outperformed kriging for the prediction of surface soil texture. This result also indicates that remote sensing data such as Landsat ETM imagery have potential as powerful auxiliary variables for improving the estimates of surface soil texture. Our results are consistent with previous findingswhich show the superiority of cokriging soil properties with remote sensing data over kriging (Zhang et al. 1992; Bishop and McBratney 2001; Sullivan et al. 2005).

The coefficient of determination (R ²) between the measured and predicted values is another method for evaluating the accuracy. The R ² values of kriging are 0.120.09 and 0.06 for sandsilt and clay respectivelywhich means that predictions generated by kriging only explain 12%9% and 6% of the variability of sandsilt and clayrespectively. Compared to MSRwhich can be applied to explain 32%21% and 36% of the variance in sandsilt and clayrespectivelykriging produces worse predictions for surface soil texture. It suggests that both MSR and kriged estimates have low reliability over a large-sized area. The variations in landscape and the low-density sampling in the study area may be the two most important factors that affect the reliability of the predicted soil texture by MSR and kriging (Wu et al. 2009). Cokriged estimates explain 62%60% and 28% of the variability of sandsilt and clayrespectivelywhich provides evidence that remote sensing data as auxiliary data for cokriging soil sand and silt contents can significantly improve predictions despite different landscapes and the low density of soil sampling in the study area. Howevercokriging produces worse predictions for clay than MSR. This can probably be attributed to the fact that sampling densities are not enough to capture the variation of the cross-semivariogram of soil clay content (the R ² between predicted and measured cross-semivariogram values is only 0.23) (Table 2). Thusthe predicted sand and silt contents by cokriging were selected to determine the soil texture according to the USDA system. Fig. 7 shows the spatial distribution of surface soil texture predicted by cokrigingas well as the measured textures of the 58 soil samples. Comparison of measured and predicted textures indicates that the cokriging method has a high correct recognition rate of nearly 85% for soil textures of the 58 soil samples. Overall, cokriging is a robust method in the prediction of surface soil texture in the study area.

Figure 7 The spatial distribution of surface soil texture generated from cokriging-predicted sand and silt contents. Different black symbols on the map denote the measured textures of the 58 soil samples.

CONCLUSIONS

In our studysoil sand content is positively correlated with Landsat ETM DNs of six bands (Bands 1–5 and Band 7)and soil silt and clay contents are negatively correlated with DNs of six bands. The correlation coefficients between surface sandsilt and clay contents and DN of Band 7 are the largest in absolute value. The prediction of surface soil texture by cokriging with DN of Band 7 is an improvement over that by MSR and kriging. This indicates that remote sensing data such as Landsat ETM imagery have potential as powerful auxiliary variables for improving the estimates of surface soil texture.

The distribution of surface soil texture is controlled by soil parent materials in the study area. Soil sand content in the northeastern parts of the areawhere most are Luvisols formed from coarse-grained residuum and slope depositsis much higher than in the area of Vertisols formed from medium- to fine-textured lacustrine deposits. The spatial distribution maps of surface soil sand and silt contents will be transferred to those of surface soil hydraulic properties for use in hydrological and ecological modeling for agricultural and environmental management in the area.

ACKNOWLEDGMENTS

The study was financially supported by the National Natural Science Foundation of China (No. 41171183) and the National Science Foundation for Distinguished Young Scholars of China (No. 40725010).