COMPARISON OF MAPS OF SPATIAL VARIABILITY OF SOIL RESISTANCE TO PENETRATION CONSTRUCTED WITH AND WITHOUT COVARIABLES USING A SPATIAL LINEAR MODEL

A study about the spatial variability of data of soil resistance to penetration (RSP) was conducted at layers 0.0-0.1 m, 0.1-0.2 m and 0.2-0.3 m depth, using the statistical methods in univariate forms, i.e., using traditional geostatistics, forming thematic maps by ordinary kriging for each layer of the study. It was analyzed the RSP in layer 0.2-0.3 m depth through a spatial linear model (SLM), which considered the layers 0.0-0.1 m and 0.1-0.2 m in depth as covariable, obtaining an estimation model and a thematic map by universal kriging. The thematic maps of the RSP at layer 0.2-0.3 m depth, constructed by both methods, were compared using measures of accuracy obtained from the construction of the matrix of errors and confusion matrix. There are similarities between the thematic maps. All maps showed that the RSP is higher in the north region.


INTRODUCTION
The intensive cultivation of the soil and the use of machinery and heavy equipment lead to degradation of the physical conditions and, consequently, increased soil compaction.WEIRICH NETO et al. (2006) and LIMA et al. (2010) indicate the moto mechanization as one of the main causes of soil compaction, which may hinder the penetration of roots.In this context the soil resistance to penetration (RSP) is a good indicator of soil compaction for being property directly related to plant growth and quick and easy determination.A compacted soil adversely affects the growth of roots, and consequently, may decrease productivity in the area.COLLARES et al. (2008) found that the presence of compacted areas on the soil restricts root growth at that point and the roots are induced to grow in another direction.OLIVEIRA et al. (2011) studied the relation between soybean yield and soil chemical properties for the selection of statistical models and found that there is autocorrelation and cross-correlation between soil resistance to penetration in the first layers and soybean yield measured over the transversal space, and that in the scenario considered by the authors, the soil bulk density is not autocorrelated or have cross-correlation with soybean yield.
Although the literature presents critical values of RSP in which occurs mechanical impedance to root development of plants, there has been disagreements to these values, considering the soil type and species studied.In this context, CANARACHE (1994) proposed a model to show that the penetration resistance values vary under different conditions, such as combinations of textures, densities, and moisture, and from this model generally concluded that the values of RSP in the interval [1.10; 2.59] MPa have few limitations to root growth, while values in the range [2.6; 5.0] MPa indicate that there are some limitations.
The objective of this reseach was to study the spatial variability of soil resistance to penetration retreat in the layers of 0.0-0.1 m, 0.1-0.2m and 0.2-0.3m deep, using geostatistical methods in univariate analysis (traditional) and using a spatial linear model (SLM) of the RSP in the layer 0.2-0.3m depth as covariables with the RSP of the layers of 0.0-0.1 m, 0.1-0.2 m.To compare the thematic maps of the RSP layer of 0.2-0.3m of both methods, it was calculated measures of accuracy obtained from the error matrix and confusion matrix.

MATERIAL AND METHODS
Data collection was accomplished in the year 2007/2008 in a commercial producer of grains with 101.47 hectares located in the city of Cascavel, state of Paraná, Brazil, with approximate geographic coordinates of Latitude 24º57'21"S and Longitude 53º34'32"W, and an average altitude of 650 m.The climate presents as mild mesothermal and super humid, Cfa -Köeppen type climate, with moderate temperatures, well-distributed rainfall and hot summer.In winter, the average temperature is below 16°C, subject to frost, and in summer the maximum exceeds 30ºC and the average annual temperature is 21°C.
The planting system used was no-till and the soybean cultivar VMax NK412113 was planted in the north, and in the rest of the area the soybean cultivar CD213 was planted.In soil, which is classified as Hapludox, clayey, and presents historical succession of oat yield in the winter and soybean in the summer, it was performed a systematic sampling centered on pairs of nearby points (lattice plus close pairs), with maximum distance of 141 m between points and in some random places the sampling was performed with distances of 75 and 50 m between points.63 sample elements were used in each layer of the study (0.0-0.1 m, 0.1-0.2m and 0.2-0.3m deep) to 101.47 ha, and the location of sampling points was performed with a GPS GEOEXPLORE 3, a Universal Transverse Mercatur (UTM) spatial coordinate system, with an accuracy of five meters.
The RSP in layers studied was measured with a penetrometer SC-60 SOILCONTROL brand, with shaft of 600 mm and 9.53 mm in diameter, equipped with a cone at the tip of 129.3 mm 2 base area, 12.83 mm in diameter and 30 degrees of angle vertex, four replications per location were made, taken at random in each plot, setting an average value of soil mechanical resistance to penetration.The layers 0.0-0.1 m, 0.1-0.2m and 0.2-0.3m depth had mean humidity of 12.74%, 14.70% and 15.17%, respectively, considering that for the three layers the coefficient of variation did not exceed 13%, characterizing homogeneity of observations of humidity in the region (GOMES, 2009).
To model the spatial structure with RSP, it was considered a Gaussian stochastic process {Z(s),sS},  a two-dimensional Euclidean space.Suppose that the data, Z(s 1 ),...,Z(s n ), from this process are recorded in known spatial locations, s i (i = 1,...,n), and generated by the model Z(s i ) = µ(s i ) + e(s i ).In this model, the deterministic µ(s i ) and stochastic e(s i ) terms may depend upon the spatial location where Z(s i ) was obtained.It is assumed that the stochastic e(.) error has zero mean, E[e(s i )] = 0, and the variation between points in the space is determined by a function covariance COV[e(s i ), e(s u )] = C(s i , s u ), and for some known functions of s, such as x 1 (s),...,x p (s), the mean of the process is µ(s (  , and is called spatial linear models (SLM), in which β 1 ,...,β p are unknown parameters and to be estimated.
Equivalently, in matrix notation, it is understood that the SLM is: In which the vector of random errors ε has E(ε) = 0 (zero vector) and covariance matrix Σ = [(ζ iu )], with ζ iu = C(s i , s u ).It is assumed that the covariance matrix Σ is not singular, that X is a nxp matrix of full rank, β = (β 1 ,...,β p ) T and Z follows a normal distribution n-variate with mean vector Xβ and covariance matrix Σ, i.e., Z ~ N n (Xβ, Σ).
Given the parametric form of the covariance matrix, In which, φ 1 is the nugget effect or error variance; φ 2 is the contribution or dispersion variance (sill); R(φ3) is a nxn matrix which is a function of φ 3 , R(φ 3 )=[(r iu )] is symmetrical with its diagonal elements r ii = 1, to i=1,...,n; φ 3 is a function of range (a) of the model and I n is a nxn identity matrix.The parametric form of the covariance matrix is isotropic for various processes (GUEDES et al., 2008), where the covariance C(s i ,s u ) is defined according to the covariance function C(h iu )=φ 2 r iu , in which h iu =||s is u || is the Euclidean distance between the points s i and s u .In the covariance functions C(h iu ), the variance of Z is C(0) = φ 1 + φ 2 , and semi variance may be defined as To identify the structure of spatial dependency between the sampling elements, it was used the classic Matheron semivariogram, as defined in Equation (3) Choosing an appropriate model is to obtain estimators of the parameter vector δ =(β T ,θ T ) T , in which β =(β 1 ,...,β p ) T and θ =(φ 1 , φ 2 , φ 3 ) T by the methods of parameter estimation for maximum likelihood (ML) and maximum restricted likelihood (MRL).To estimate the parameter vector  , it was chosen a vector  ˆ=(  ˆT, ˆT) T which maximizes the likelihood function in the field, in which  is the parameter space.Considering the stochastic process Z =(Z(s 1 ),...,Z(s n )) T , in which Z ~ N n (Xβ, Σ), the method of ML estimation of δ consists in maximizing the logarithm of the likelihood function The MRL method used to estimate the parameters of the covariance matrix Σ, consists in obtaining less biased estimators than the estimators obtained by ML.The method of MRL estimation of θ which consists in maximizing the logarithm of the restricted likelihood function (BORSSOI et al., 2009), To choose the model that best fits the data; it was used the cross-validation technique and the maximum value of the logarithm of the likelihood function (FARACO et al., 2008, JOHANN et al., 2010;BORSSOI et al., 2011).After the selection, it was used the vector of parameters  ˆ for the universal kriging interpolation.The comparison between the univariate map (traditional method), called the reference map, and the map built by SLM, called model map, of the RSP in the layer 0.2-0.3m depth was performed using measurements obtained from the error matrix (Table 1) which has as measurement unit the pixel.Each matrix element represents the total number of pixels or the area belonging to the class i of the model map and to the class j of the reference map, being N the total number of pixels or the total area.TABLE 1. Error matrix of the reference map in relation to the model map.

Reference Map
The main diagonal is the case in which the pixels (or areas) had the same classification in the two maps, while the elements outside the main diagonal represent the classifications that do not coincide.The global accuracy (GA), , is a measurement used to measure the similarity between the reference map and the model map and according to ANDERSON et al. (2001), the minimum level of accuracy is 0.85.It is possible to build a range of (1-α) % confidence for the form GA, is the nominal value of a random variable with standard normal distribution (FOODY, 2009).
The Kappa index (K) (BAZZI et al., 2008, FOODY, 2009;FOODY, 2010) has been used to measure the accuracy of thematic classifications, recommended as an appropriate measure of accuracy for using all elements of the error matrix.Provides a measure of agreement between the values of the reference map and the values from the model map, being defined by Equation ( 4 is the number of pixels (or area) of the row i, x is the number of pixels (or area) of the row i and column i.According to the classification of KRIPENDORFF (1980), K is classified with low accuracy if K < 0.67, average accuracy if 0.67  K < 0.80, and high accuracy if K  0.80.The variance of the Kappa index is obtained using Equation (5), In which, MA & REDMOND (1995), proposed an index that shows less subjective and easier to understand and use, the concordance index Tau (T), also known as Kappa modified, defined by T=(θ 1 -p i )/(1-p i ), being pi the a priori probability for each class i.When the a priori probabilities for the classes are the same, there is p i =1/r, in which r represents the number of classes of the error matrix.The index Tau can follow the same classification of Kappa.The variance of Tau index is calculated with the same formula of variation of Kappa, and the difference is the inclusion of a priori probability.Another way to compare thematic maps is by setting up the confusion matrix by class, which is obtained by the elements of the error matrix, presented in Table 2.
In Table 2, a i = ii x is the amount of pixels (or area) of a class i of the reference map which was correctly classified as belonging to class i on model map (true positive - is the amount of pixels (or area) that does not belong to class i of the reference map and was classified as belonging to class i in the model map (false positive -FP i ); is the number of pixels (or area) of a class i of the reference map and which belongs to a different class in the model map (false negative -FN i ) and d i = N -(a i +b i +c i ) is the amount of pixels (or area) that does not belong to class i in the reference map and was classified as not belonging to class i on model map (true negative -TN i ).TP i and TN i values are correct predictions in each class.The values FP i and FN i are considered estimation errors, FP i known as errors of commission or overestimation, and FN i known as errors of omission, for i = 1, ..., r.
Table 3 shows measurements of the confusion matrix per class.The sensitivity index (S i ) is a measure indicating the probability that a unit area on the model map is classified as belonging to class i (i = 1,..., r) if it really belongs to the class i of the reference map.The specificity index (E i ) indicates the probability that a unit area does not belong to the class i of the reference map or to the class i of the model map.The errors of commission (FPR i ) indicate the proportion of area units that do not belong to the class i of the reference map and that belong to class i of the model map.Errors of omission are considered serious as it indicates the proportion of area that belongs to the class i on reference map and has been classified in another class in the model map.TABLE 3. Measures obtained from the confusion matrix for the i-th class (i =1,..., r).
The negative estimation power (NEP i ) indicates the probability of a unit area not to be classified by both reference map and model map since it was classified as absent by the model map.Positive estimation power (PEP i ) indicates the probability of a unit area to be classified by both reference map and model map since it was classified as present by the model map.The conditional Kappa index (CKU i and CKP i ) is presented as an additional means of incorporation of hit by hypothetical chance in the evaluation of the accuracy per class, and the conditional Tau index (CTU i and CTP i ) is an alternative to this, for the Kappa may provide wrong conditional estimates for not knowing the distribution of reference data in advance.The values of UA i and PA i , called user accuracy and producer accuracy, are given respectively by Table 5 presents measurements of global confusion matrix.The interpretation for the global indexes of total sensitivity (S) and total specificity (E) is analogous to the versions by class.The values of Matthews correlation coefficient (MCC) are within the range [-1, 1], and the value 1 corresponds to a perfect prediction, the value 0 corresponds to a random prediction and value -1 corresponds to an inverse prediction.The average mutual information AMI (FINN, 1993) calculates the conditional probability that an area on a map belongs to a particular class given the class of this area on the second map, thereby providing means to measure the similarity between maps with different themes.
For data analysis, it was used the software R (R DEVELOPMENT CORE TEAM, 2009) and its modules geoR (DIGGLE & RIBEIRO JR, 2007).

RESULTS AND DISCUSSION
Table 6 presents descriptive statistics of the data of RSP [MPa] in the different layers.According to the average values, it is verified that there is little restriction in the growth of roots (CANARACHE, 1994).According to the coefficient of variation (CV) there is an average data homogeneity of RSP [MPa] in three layers, for 10% < CV < 20% (GOMES, 2009).According to the coefficient of asymmetry and kurtosis, the RSP in the three layers have characteristics of normal distribution, according to criteria shown in JONES (1969).In this study no evidence sample was discarded.Table 7 presents the results of univariate geostatistical analysis with the models adjusted and the estimated parameters of the RSP for the layers under study.According to the cross-validation criteria and the maximum value of the logarithm of the likelihood function (MLL), for the layer of 0.0 -0.1 m depth the model which best fit was the Gaussian with parameters estimated by ML method.To the layers 0.1-0.2mand 0.2 -0.3 m, the model which best fit was the Gaussian with parameters estimated by MRL.The ranges (a) obtained represent the distance at which the sampling points are correlated, ranging from 1274.35 m to 1381.60 m, not exceeding the maximum distance of 1710.26 m.Analyzing the coefficient of relative nugget effect (RNE), according to the scale of CAMBARDELLA et. al (1994), we found a moderate spatial dependence (25% < RNE  75%) for the RSP in the layers of 0.0 -0.1 m and 0.1 -0.2 m deep and weak spatial dependence (RNE > 75%) to the RSP in the layer of 0.2 -0.3 m deep.Table 8 presents the results of cross validation for the adjusted models for the three layers in the study.Note that the values of the mean error (ME) and reduced error (RE) are close to zero.The standard deviation of the mean error (S ME ) and absolute error (AE) are the smallest among the models studied.And yet, the standard deviation of the reduced error (S RE ) is very close to one, as expected (FARACO et al, 2008).Figure 2 presents the thematic maps of RSP built by ordinary kriging.There was a decrease in soil compaction in accordance with increasing depth.Figure (2a) shows the thematic map of the RSP in the layer of 0.0 -0.1 m deep, where 17.86% of the area has RSP greater than 2.88 MPa, 18.36% between 2.61 and 2.88 MPa, 42.36% between 2.34 and 2.60 MPa, and 21.42% between 2.07 and 2.33 MPa. Figure (2b) shows the map of the RSP in the layer of 0.1 -0.2 m deep, in which 30.07% of the area has RSP between 2.34 and 2.60 MPa, 65.09% between 2.07 and 2.33 MPa, and 4.84% less than 2.07 MPa. Figure (2c) shows the map of RSP in the layer of 0.2 -0.3 m deep, and 26.69% of the area has RSP between 2.07 and 2.33 MPa, and the remaining 73.31% presents RSP less than 2.07 MPa.Although there are questions about values of soil resistance to penetration considered critical to crop development, CARVALHO et al. (2006) reported that RSP values ranging between 1.3 and 2.9 MPa did not limit grain yield of the bean in an Oxisol.Considering the conditions of this study and the criteria for CANARACHE (1994), it appears that for the three layers, the RSP is higher in the North region, and 36.22% of the soil has limited the roots growth up to 0.1 m deep, featuring a compacted soil.It should be noted that the layers samples showed little change in water content of the soil.
In the analysis using spatial linear model (SLM) the interest was to study the spatial variability of the mean ) ( X 2 (RSP in the layers from 0 -0.1 m) and X 3 (RSP in the layers 0.1 -0.2 m).The model that best fit was the Gaussian by the ML according to the cross-validation criteria and the maximum value of the logarithm of the likelihood function, LML = 16.66.The SLM obtained for the mean of the RSP in the layer of 0.2 -0.3 m deep was: , with spatial dependence parameters obtained of the covariance matrix Σ, in which 1  =0.0326; 2  =0.0019 e 3  =12.694, with a range estimative of 21.99 on relative nugget effect (RNE) equal to 94.5%, characterizing weak spatial dependence.Table 9 shows the values obtained for the cross-validation criteria for the adjusted model with covariables, which followed the same criteria for analyzing the values shown in Table 8.The thematic map of Figure 3 shows the RSP to be higher in the North region, but the soil has no limitation to root growth.Table 10 presents the error matrix associated with the RSP maps in the layer 0.2 -0.3 m deep, in which it was considered as a reference map the univariate map obtained by ordinary kriging and as a model map the map obtained by the SLM with universal kriging.In parenthesis is the value of area in ha.
It is noteworthy that both maps were generated with 16,204 pixels, which means that each pixel represents an area of approximately 62.62 m 2 .
Table 11 presents the estimates of the accuracy measures and their corresponding 95% confidence interval for the comparison of the maps.The index of global accuracy (GA) obtained was 0.98, indicating an acceptable precision according to ANDERSON et al (2001).According to the classification of KRIPENDORFF (1980), Kappa (0.94) and Tau (0.95) indexes are regarded as high accuracy.Thus, it can be seen that there are few differences in the classification of the thematic maps studied.Table 13 presents measures of accuracy calculated by class of RSP in the layer of 0.2 -0.3 m deep.In class [1.47 -2.06] MPa, the sensitivity index S 1 = 0.973 is the probability that a pixel belonging to this class in the two kinds of maps constructed, the specificity index E 1 = 0.992 is the probability a pixel does not belong to this class in both types of maps constructed.The omission error (FNR 1 ) indicates 2.7% of the pixels belonging to this class of the reference map, belongs to a class different from the model map.The commission error (FPR 1 ) indicates 0.8% of the pixels that do not belong to this class in the reference map, belongs to this class in the model map.The positive estimation power (PEP 1 ) indicates 0.997 and is the probability of a pixel to be classified by both the reference map and the model map since it was classified as present by the model map.The negative estimation power (NEP 1 ) is 0.923 indicating the probability that a pixel not to be classified by neither the reference map nor the model map since it has been classified as absent by the model map.
In class [2.07 -2.34] MPa, the sensitivity index S 2 = 0.992, is the probability of a pixel to belong to this class of both types of maps constructed, the specificity index E 2 = 0.973 is the probability a pixel to not belong to this class in both types of maps constructed.The omission error (FNR 2 ) indicates 0.8% of the pixels belonging to this class of the reference map, belongs to another class different from the model map.The commission error (FPR 2 ) indicates 2.7% of pixels that do not belong to this class in the reference map, belong to this class in the model map.The index PEP 2 indicates 0.923 is the probability that a pixel is classified by both the reference map and the model map since it was classified as present by the model map.The index PEN 2 indicates that 0.997 is the probability that a pixel not to be classified by both the reference map and the model map since it has been classified as absent by the model map.The global confusion matrix for the RSP in the layer of 0.2-0.3m deep presented 15,843 pixels classified as (+) (+), 15,843 as (-) (-) and 361 as (+) (-) and (-) (+).The global indexes of sensitivity (S), specificity (E), Matthews correlation coefficient (MCC) and Average Mutual Information (AMI) were 0.978, 0.978, 0.955 and 0.455, respectively.Thus, the sensitivity (S), which is the conditional probability that a pixel in the reference map belongs to a model map is 0.978.The specificity index (E) of 0.978 indicates the probability of not belonging to the reference map and to the model map, the Matthews correlation coefficient (MCC = 0.955) shows that there is a direct proportional correlation in the classification of thematic maps.The index of average mutual information AMI = 0.455 indicates the similarity between the reference map and the model map.

CONCLUSIONS
By geostatistical methods, it was possible to detect that the RSP in three layers showed higher degrees of compression in the North region.Through the thematic maps constructed, it was possible to determine the areas where it is necessary to apply techniques to reduce the problem of soil compaction.The modeling by SLM of RSP in the layer 0.2-0.3m deep is very important to estimate the RSP in terms of previous layers in sectors not sampled in the experimental area.The accuracies measures presented allowed realization of the comparison between the univariate thematic map using ordinary kriging and the thematic map using the SLM with universal kriging.The index values showed the existence of similarity between the constructions of both maps.

Figure 1
Figure1presents the theoretical semivariograms which were best adapted to the experimental semivariograms for each layer, according to the criteria for validation of models studied.


of the RSP in the layer 0.2 -0.3 m deep as a function of covariables

FIGURE 3 .
FIGURE 3. Thematic map of RSP [MPa] at 0.2 -0.3 m according to the preceding two layers.

Table 3 ,
it is possible to compare for each class the map reference with the model map.For a global comparison of the maps, JENNESS & WYNNE (2005) present the global confusion matrix, shown in Table4.

TABLE 5 .
Measures derived from the global confusion matrix.

TABLE 6 .
Descriptive statistics of the RSP [Mpa] in the three layers of depth.

TABLE 7 .
Adjusted spatial models and the estimated parameters of RSP[MPa].
 ))x100 relative nugget effect; LML: maximum value of the logarithm of the likelihood function; ML: maximum likelihood; MRL: maximum restricted likelihood; Gaus: Gaussian model.

TABLE 8 .
Cross-validation criterion for the adjusted models.
ME: mean error; RE: reduced error; S ME : standard deviation of mean error; S RE : standard deviation of reduced error; AE: absolute error

TABLE 9 .
Cross-validation criterion for the adjusted models with covariables.RE: reduced error; S ME : standard deviation of mean error; S RE : standard deviation of reduced error; AE: absolute error ME: mean error;

TABLE 10 .
RSP error matrix (0.2 -0.3 m) of the number of pixels and area (ha) by class.

TABLE 11 .
RSP accuracy measures in the layer 0.2-0.3m deep.

TABLE 13 .
Accuracy measures by class of RSP [MPa] in layer 0.2-0.3 m.