Geostatistical Fusion of LiDAR and Hyperspectral Data to Improve the Spatial Prediction of Sand Content in a Forested Catchment of Southern Italy

1. Introduction

Many functions are carried out by soil, ensuring a fundamental contribution to human life and well-being, however, to understand the potential of soil in providing functions, it is necessary to characterize the variability of soil properties in space and time [1]. Soil texture is one of the most important properties that influence many soil functions and processes through both the absolute values of its properties and their spatial variability [2]. In particular, it controls variations in soil moisture, availability of water to crops, erosion and aggregation. In addition to direct measurements in the laboratory, different types of soil sensing technologies, from laboratory reflectance diffuse spectroscopy to proximal and remote sensing, can provide rapid and reliable information on soil properties variability at different spatial and temporal scales [3,4,5]. However, data measured in laboratory and/or with different soil sensing methods require their fusion and integration [6,7,8]. A proper integrated use or synthesis (fusion) of different spatial data results to be more informative than the one coming from individual sources and provides new knowledge, understanding, or explanations of the processes [9,10,11].

When mapping and modeling topsoil there has always been a compromise between the need to produce models of high spatial resolution but limited extent (at local scale) and the need to produce models of coarse spatial resolution but wide extent (as at regional or country scale) [12]. If extended spatial models could also be produced at high resolution, this would undoubtedly mean major advances in understanding the physical and biological processes that take place on the Earth's surface [13]. Morphometric attributes as elevation, slope gradient, aspect, local relief, slope curvature, have been commonly used to estimate the surface lithology and model the spatial variability of soils [14,15]. However, most of these studies used digital elevation models with a spatial resolution greater than 10 m. Although such a scale may be suitable for some environmental studies, it may be unsuitable for accurately describing the spatial variability of highly changing landscapes [16], such as the characteristic landscapes of the forested mountain of some Calabrian areas in southern Italy [17]. The climate, the topography, the vegetation cover, and land use are the main factors that influence soil properties and determine its heterogeneity at many scales [18].

One of the most important recent developments in remote sensing for lithological mapping of the soil surface is the increasing availability of high-resolution digital elevation models (DEM)s derived from LiDAR (Light Detection and Ranging) elevation data. The high density of measurement points obtained from a LiDAR sensor allows for precision mapping of land surface features and also makes it possible to clearly distinguish soil from vegetation [12]. LiDAR data are widely used to study landscape [19], topography [20], tree attributes ([21,22] phenotyping [23], and biomass estimation [17] but can be used to assess soil properties [24]. The use of LiDAR also enables to generate large-scale DEMs with sub-meter resolution in extensive heterogeneous mountain forested areas [25]. Recently, the DEM resolution derived by LiDAR data has increased reaching even 0.5 m with the use of airplanes [26] up to very accurate resolutions (~0.03 m) using the drone [27]. Airborne Laser Scanning (ALS) is the most highly accurate and efficient method to acquire 3D data from large areas and to generate DEMs, such as Digital Terrain Model (DTM) and Digital Surface model (DSM). Moreover, using geographic information systems (GIS) algorithms, these new high-resolution elevation data can be used to provide several primary and secondary topographic attributes.

It is widely accepted that a larger number of soil samples can more accurately describe and model spatial heterogeneity, however the limited budget for soil measurements generally restricts the number of samples when using standard laboratory techniques as they are generally expensive and time consuming [28,29]. As an alternative soil spectroscopy techniques, such as visible-near-infrared spectroscopy (Vis-NIR, 350-2500 nm), can be used to measure soil properties such as the different particle sizes (clay, silt and sand), organic and inorganic carbon content and iron content quickly, accurately and relatively inexpensively [30,31]. This allows the use of soil data from hyperspectral spectroscopic measurements, appropriately reduced in dimensions, as an additional source of information supplementing direct laboratory soil measurements, in order to improve the prediction of the variable of interest.

To produce soil maps associated with some measure of prediction uncertainty by efficiently combining (fusing) heterogeneous data, several techniques spanning from traditional spatial statistics, geostatistics to machine learning have been used in recent decades [32,33,34,35]. Actually, the different computing technologies lead to a single formulation of spatial modelling expressed by the following formula [36]:

$$Z(x)=\mu (x)+\epsilon \text{'}(x)+\epsilon \text{'}\text{'}$$

(1)

where:

x: location vector in one, two or three dimensions,

Z(x): the random variable Z at location x,

μ(x): deterministic component or trend (drift),

ɛ’ (x): stochastic component, i.e. spatially dependent residual from μ(x)

ɛ’’: non spatially correlated component (white noise or unexplained variability).

Such a formula makes it possible to model the spatial variability of any soil property on different spatial scales, separating the deterministic component (at long scale) from the stochastic one (at short or local scale). Non-geostatistical techniques (multi-linear regression or machine learning) have generally been used to describe the trend, while geostatistical techniques to quantify the stochastic component of variability. These two types of approaches can be kept separate or combined in hybrid techniques (non-stationary geostatistics) [37]. One of the latter is kriging with external drift [38,39], when the trend is externally modelled by auxiliary variables and trend and residuals are simultaneously estimated in a single system with a jointly calculated covariance function. However, the trend of the target variable or some of its transformation functions must be expressed by linear models.

The objectives of the study were: 1) to evaluate the ability of a multivariate approach based on non-stationary geostatistics to merge remotely sensed high-resolution LiDAR-derived topographic attributes with Vis-NIR diffuse reflectance spectroscopy and laboratory data to produce high-resolution maps of soil sand content and its estimation uncertainty in a forested catchment in southern Italy, and 2) compare this approach with the commonly used univariate approach of ordinary kriging.

2. Materials and Methods

2.1. Study area

The study area is the Bonis catchment, situated in the upland landscape of the Sila Massif (Figure 1), Calabria Region, (southern Italy), which is an important site for investigating the effect of forest management on its hydrologic behavior. Due to its location, the installed instrumentation, and its forest cover, the Bonis catchment has been studied intensively in many years [17,40,41,42]. The catchment has an area of 140 hectares and its altitude varies from 1019 to 1341 m a.s.l., with a mean value of 1130 m a.s.l. (Figure 1). The forest cover is mainly characterized by Pinus nigra ssp. laricio (Poir.) Maire (about 95 hectares), whereas the remaining 45 hectares are formations of Castanea sativa Mill and riparian vegetation and a small portion of bare soil.

In the catchment, crop out Palaeozoic granitoid rocks, which are deeply fractured and weathered and often covered with thick layers of saprolite and/or colluvial deposits [42]. The landscape is predominantly characterized by rugged morphology with steep slopes, often cut by deep incisions [42]. Slope gradient varies from 0° to about 50° with a mean value of 21°. According to the Calabria Soil Map [43], the soils of the study area can be classified as Typic Xerumbrepts and Ultic Haploxeralf [44]. The dominant soil textural classes are sandy loam and sandy clay loam. The upper A-horizon often shows a very dark brown color due to the accumulation of organic matter (umbric epipedon). Regarding to the pedo-climate features, the study area is characterized by a mesic soil temperature regime and an udic soil moisture regime [43].

2.2. Remote Sensing (LiDAR) data acquisition and processing

In the study site, the Light detection and ranging (LiDAR) scanning system was used by Airborne Laser Scanning (ALS) to provide a high-quality point cloud and density.

The measurement was made with a Riegl LMS-Q560 laser scanning with frequency 400,000 Hz, FOV 60°, and inclination of 20°mounted on a helicopter operating at approximately 700 m above ground level. The LiDAR sensor has been set up to obtain the resulting point cloud about of 20 points per square meter to attain the accuracy of producing DEMs layers with 1m resolution. The World Geodetic System 1984 (WGS 84) and Universal Transversal Mercator (UTM) projection (zone 33N) coordinate system were assigned to the cloud points. The LiDAR data used in this research was acquired in 2017 as part of a national project [17].

To derive the DEM accurately, the raw point cloud acquired by the sensor must be processed to remove unwanted points from laser range scans and obtain the cleaned point cloud. The LiDAR data was processed using the commercial software LiDAR360 (Green Valley International Inc., California, USA) and the data has been archived in the LAS (LASer) format.

The workflow to generate the terrain models involves several steps, the first one is to check the data quality and remove the isolated points and noise through denoising. The Remove Outliers tool in LiDAR360 was used to remove the noise and improve the quality of the data. The tool analyses the average distance to all neighboring points and determines average and standard deviation for each point so to identify the outliers and remove them from the data. Once outliers have been removed from the point cloud the next was to classify ground points through Classify Ground Points tool in LiDAR360 and set up the value 2 to the ground points as defined by the American Society for Photogrammetry and Remote Sensing (ASPRS) for Standard LIDAR Point Classes. After this, the DEM tool has been used to create the raster of DTM with 1-m resolution using Inverse Distance Weighting (IDW) interpolation method.

2.3. Topographic attributes

The topographic attributes, used as auxiliary variables in the geostatistics analysis to produce high-resolution maps of soil sand content distribution, are the following: elevation, slope, aspect, curvature, length-slope (LS) factor, stream power index (SPI), and topographic wetness index (TWI). These topographic attributes were derived by the LiDAR-DTM, using the SAGA-GIS software [45].

The soil properties distribution can be influenced by elevation because it affects local microclimate by changing the pattern of temperature, precipitation, and soil moisture [46].The slope gradient is a morphometric attribute of primary importance in the processes governing both pedogenesis and soil erosion/deposition; in fact, it influences surface runoff, soil infiltration, and intensity of erosion processes [47,48,49]. The aspect plays a fundamental role in controlling soil moisture along the slopes through solar radiation and rainfall as well as in influencing many factors regulating soil formation processes and soil productivity [46,47]. Slope curvature provides information on slope shape; it is important for soil mapping because it influences local superficial water flowing in terms of convergence or divergence and of acceleration or deceleration of the flow across the surface. LS factor is a variable used in the RUSLE equation to consider the effect of topography on soil erosion [50]. It depends on slope steepness factor (S), and slope length (L); it influences surface runoff intensity and is considered as a sediment transport capacity index. SPI attribute is a measure of the erosive power of overland flow based on the assumption that discharge is proportional to the specific catchment area [51]. Generally, the higher value of SPI the higher the likely occurrence of water soil erosion. TWI attribute indicates the effects of local topography on the hydrological processes [47]. It is an index of the likelihood by a cell of collecting water and is considered a predictor of soil saturation. It exerts a great influence on hillslope processes such as soil and water redistribution and vertical infiltration. Finally, TRI is a morphometric indicator that describes whether the topography of an area is flat or undulated and represents the spatial variability of a landscape. It may be used to measure the variation in elevation between neighboring spatial pixels of a DEM [47].

2.4. Field soil sampling and laboratory analysis

For this study, a field survey was conducted to collect 135 topsoil (0 –0. 20 m depth) samples within the Bonis basin (Figure 1). At each sampling location surface litter was removed, and soil was sampled using a metallic core cylinder having a diameter of 0.075 m and a height of 0.20 m. The position of the soil sampling sites was acquired using a hand-held GNSS device, with a precision of about 1 m.

The soil samples were transferred to the laboratory in polyethylene bags, dried at 40°C in ventilated oven, homogenized and then sieved using a 2-mm mesh stainless steel sieve. Successively, the sieved samples were used for soil particle size measurement with hydrometer method after a pre-treatment with sodium hexametaphosphate as a dispersing agent [52]. According to USDA [53] classification for soil texture, sand, silt, and clay particle sizes were categorized as follows: 0.05–2.0 mm for sand, 0.002–0.05 mm for silt, and 0.002 mm for clay. In addition, the soil samples were finely sieved through a 0.25 mm mesh sieve, and SOC content was determined by dry combustion using a Shimadzu TOC-L analyzer with an SSM-5000A solid sample module (Shimadzu Corporation, Kyoto, Japan). Each measurement was made in triplicate, and the soil sample was re-measured if the standard deviation of the three replicates exceeded 2.5%.

2.5. Hyperspectral spectroscopy measurements

The visible near-infrared (Vis-NIR) reflectance spectroscopy analysis was performed in laboratory using an ASD FieldSpec IV spectroradiometer (Analytical Spectral Devices Inc., Boulder, Colorado, USA) with a spectral range of 350–2500 nm. The spectral measurements were performed in a dark room to reduce the effect of external light with the spectroradiometer held in the nadir position at a distance of 0.10 m from the soil sample. A 50 W halogen lamp with a zenith angle of 30° at a distance of 0.25 m from the sample was used [30]. Each soil sample was placed in a Petri dish (0.09 m in diameter and 0.01 m deep) and the surface of the soil was leveled with a spatula [42]. Before starting the soil samples measurements, dark current was removed and the spectroradiometer was calibrated with a white panel of known reflectance (Spectral on Diffuse Reflectance Panel). Generally, the calibration was repeated once every 20 minutes.

For each soil sample 50 scans were acquired, which were then averaged to obtain the soil spectrum. To reduce the noise level and external interference, the diffuse reflectance spectra were smoothed by means of a nine-point smoothing function [54], using the ViewSpecpro software (Malvern Panalytical Ltd, UK). Finally, each soil spectrum was resampled every 10 nm wavelength to reduce the total number of raw reflectance bands and avoid overfitting [55].

2.6. Preliminary statistical data analysis and non-stationary geostatistical approach

An exploratory data analysis was carried out for all study variables before applying the geostatistical approach. The main basic statistics were calculated (min, max, mean and median values, standard deviation, skewness, kurtosis) and those variables with a significant deviation from the Gaussian distribution (probably due to large outliers) were transformed to the standardized normal distribution. To do that, a set of Hermite polynomials truncated to the first 100 terms was used to estimate the Gaussian anamorphosis function of transformation [39]. At the end the gaussian transformed estimates were back-transformed to the raw values by using such estimated function [56].

Two datasets were required in order to apply the non-stationary geostatistical technique known as kriging with external drift (KED) [47,48], to estimate the soil sand content at the nodes of the interpolation grid (mesh = 1 m) based on the digital elevation model (DEM) obtained from the LiDAR data.

The first of these is the co-regionalization data set (vector file of points) including the sand content (target variable) data and all the auxiliary variables (external variables for KED) measured/calculated at the 135 soil sampling locations and used to calculate the trend.

The second one is the interpolation data set (raster file) based on 1-m mesh DEM at whose nodes all the auxiliary variables have been calculated or estimated.

Preliminarily to the realization of the two aforementioned data sets, it was necessary to apply fusion techniques since the data to be processed were of different type and not in the same spatial location (not-collocated): 1) scattered sample data (at 135 locations): clay and soil organic carbon (SOC); 2) sample data at the same 135 locations but consisting of diffuse reflectance spectra discretized into 216 wavelengths intervals and to which principal component analysis (PCA) was applied in order to reduce their size to few principal components (PCs); 3) raster data at the 1- m DEM nodes (elevation and all derived topographic attributes).

For the realization of the collocated point data set, a migration of the raster data from the pixel to the nearest sample location was therefore necessary. Differently for the raster data set, the point data were required to be interpolated at the grid nodes, using ordinary kriging for clay content and SOC data, and ordinary cokriging [39] for the principal component (PCs) of spectroscopy data.

For a in-depth description of the basic geostatistical techniques, as Linear Model of Coregionalization (LMC), kriging and cokriging, the interested reader is referred to the numerous manuals on the subject [37,39,57]. Following there is a brief description of the main statistical and geostatistical procedures used in this work.

2.6.1. Principal Component Analysis

Principal component analysis (PCA) [58] is a method of analysis for multivariate data, which is widely used because of the simplicity of its algebra and straightforward interpretation (Wackernagel, 2003). Moreover, it is frequently applied in remote sensing to reduce the dimensionality of radiometric variables recorded by multi or hyperspectral radiometers. The principal components (PCs) are calculated as the eigenvectors of the correlation matrix [59] and are not then directly observable variables for which they require to be interpreted from a scientific-rational perspective. Mathematically, they are expressed as a linear combination of the observed variables and are of the same number. Actually, only a part of them is necessary and namely those with eigenvalue greater than 1 determine the number of retained PCs (Kaiser's criterion) [60]. The results can be subjected to rotation by a variety of methods, those employing orthogonal rotations are preferred as they preserve the independence of the PCs. In this study, the VARIMAX procedure was used and the values of the loading coefficients were multiplied by 100 and rounded to the nearest integer and those with absolute values greater than.68.73 (the root mean square of all the values multiplied by 100 in the matrix of PC loadings) were plotted in a graph.

PCA was performed with the SAS/FACTOR procedure of the statistical software package SAS/STAT [61].

2.6.2. Kriging with external drift

Non-stationary geostatistics is applied in cases where the target variable exhibits systematic variations (trend or drift) so that estimation of the variogram model may be problematic or in some cases impossible. The variogram is a mathematical model which describes the spatial continuity of the attributes of interest and measures the average dissimilarity between observations depending on separation distance and direction [57]. Its main features are summarized by two parameters: sill and range. Generally, after an initial increase, the variogram will reach a maximum (sill variance) at a finite lag distance (range), which is the variogram value over which pairs of values are spatially correlated. Therefore, the range is the distance at which the spatial correlation becomes zero [57]. Moreover, the variogram, can show a discontinuity at its origin, called nugget effect, which characterizes the very short-scale variability within the shortest sampling interval and the errors of measurement [57]. In the non-stationary case, the sill is increasing without any finite distance within the study area. Non-stationary geostatistical techniques (universal kriging, UK, and kriging with external drift, KED) [39] are based on the formula (1) given in the introduction: a generic random variable is constituted by a deterministic component plus a random component. Universal kriging represents the trend as a linear combination of generally polynomial functions of the coordinates while KED as a linear combination of independent external variables. However, combinations of the two trend types are also possible.

Matheron [62] introduced the kriging of Intrinsic Random Functions of order k (IRF-k), which essentially consists of two steps: 1) trend identification, and 2) determination of the spatial structure of residuals through the definition of generalized spatial covariance models using higher order increments, compared with linear increments of ordinary kriging. In this way, the trend can be completely filtered out and thus stationarity again achieved [39]. Reference is made to the text cited before [62] and to other publications on the subject [41,63] for a complete description of the theory. Only the main operational steps for producing the spatial maps of the estimated variable of interest and estimation uncertainty will be outlined below.

Unlike universal kriging, in IRF-k theory it is not necessary to know the coefficients of the trend but only its order, which is generally assumed is less or equal to 2.

Structural analysis according to this theory consists of two steps:

1) Determination of the order k of the trend;

2) Calculation of the generalized covariance function K(h) of the module of the distance vector (h) and fitting of an authorized parametric model to it.

The order of the trend (k) corresponds to the degree of the polynomial used to describe the large-scale variation (i.e., at a scale longer than the size of the study area). For determining the degree of trend, the least-squares errors corresponding to the various polynomials of degree 0, 1 and 2 are calculated for each sampling point. The degree corresponding to the minimum error is assigned rank 1 and that with the maximum error rank 3. The ranks for each degree are finally averaged over all sample points and the degree corresponding to the smallest average rank is assumed to be the optimal degree for the trend.

In this procedure the neighborhood is split into two rings: the closest samples to the target sample belong to the ring numbered 1, the other samples to ring number 2. Fitting is based on a cross-validation procedure. All the data from ring 1 are used to fit the functions corresponding to the different trend hypotheses. Each datum of ring 2 is used to evaluate the quality of the fit. Then the roles of both rings are inverted. The best fit corresponds to the minimal average variance of the cross-validation errors.

A convenient model for generalized covariance function is the polynomial model given by the linear combination of a given set of generic basic structures. All possible combinations can be reduced in practice to a combination of four basic models with terms arranged according to increasing regularity (Chilès and Delfiner, 2012):

$$K(h)=C_0\delta (h)-b_0\left|h\right|+b_S{\left|h\right|}^2\log \left|h\right|+b_1{\left|h\right|}^3$$

(2)

where δ(h)=0 for h=0, otherwise δ(h)=1. The coefficients, C₀, b₀, b_s and b₁ in a two-dimensional space R², must satisfy the following set of inequalities for K(h) to be a valid generalized covariance of IRF-k (authorized model):

$$C_0\ge 0,b_0\ge 0,b_1\ge 0,b_S\ge -\frac{3}{2}\sqrt{b_0{b}_1}$$

(3)

Unlike the variogram, the generalized covariance function cannot be estimated directly, but depends on knowing the order of the trend in the IRF-k theory.

The procedure of selection of the generalized covariance function is based on a cross-validation technique performed using the two rings of samples as previously defined when determining the optimal degree of the trend. The criterion used is to compare the ratio (jackknife estimator) between the experimental and the theoretical variances: the closer this ratio to 1 the better the result.

In the case of the external drift the technique consists of replacing the large scale trend function, previously modelled as a low order polynomial, by a linear combination of a few deterministic functions of external variables (auxiliary variables). However, also a combination of the two types of trend is allowed.

Finally, kriging applied to an IRF-k provides an optimal solution to the case in which it is necessary to filter out from the estimation error any linear external function (trend) assumed to be known at each point of the interpolation grid. This is mathematically equivalent to imposing that the kriging estimator respects not only the polynomial trend of the coordinates, but also the spatial variation of the auxiliary function(s) of 'external drift'. It is therefore necessary to add a new constraint to the set of weights ${\lambda }_{\alpha }$, given by:

$${\displaystyle \sum _{\alpha =1}^N{\lambda }_{\alpha }FF\left(x_{\alpha }\right)}=FF\left(x_0\right)$$

(4)

where $FF\left(x_{\alpha }\right)$ and $FF\left(x_0\right)$ represent the values of the external drift function at the N points x_α, where the experimental data exist, and at all points x₀, where the value of the original IRF-k is to be estimated, respectively.

2.7. Mapping methods comparison

The comparison between KED and ordinary kriging (OK) was carried out in cross-validation using the two statistics suggested by Carroll and Cressie [64]: the mean error (ME) and the root mean squared standardized error by standard deviation of kriging (RMSSE):

$$ME=\frac{1}{n}{\displaystyle \sum _{\alpha =1}^n\left[z^{*}(x_{\alpha })-z(x_{\alpha })\right]}$$

(5)

$$RMSSE=\sqrt{\frac{{\displaystyle \sum _{\alpha =1}^n\frac{{[z^{*}(x_{\alpha })-z(x_{\alpha })]}^2}{{\sigma }^2(x_{\alpha })}}}{n}}$$

(6)

where $z^{*}(x_{\alpha })$ is the estimated value at location x, $z(x_{\alpha })$ the measured value at the same location, and $\sigma (x_{\alpha })$ is the standard deviation of kriging (OK or KED) at the same location (x). ME was used to assess the unbiasedness of the estimation and its optimal value should be about zero; RMSSE was used to assess estimation accuracy and should be about 1.

For the comparison, also the correlation coefficient between observation and estimation as well as the one between standardized error and estimation were calculated.

3. Results and discussion

The basic statistics of the soil attributes and elevation measured at 135 points are reported in Table 1. As can be seen, sand, silt and clay exhibited sufficiently symmetrical distributions, while SOC showed a great positive asymmetry. Therefore, it was decided to transform SOC through a Gaussian anamorphosis process before interpolating it with ordinary kriging.

Regarding the diffuse reflectance spectra measurements, in Table 2 are reported the results of PCA. Only the principal components (eigenvectors) with an eigenvalue greater than 1 are reported in Table 2. It is to note that the first PC explains more than 85%, due to the high correlation among the reflectance at the different bands, while the second explains only less than 10% of the total variance. Since each of the remaining four PCs explains just or less than 2% of the total variance, it was decided, in order to simplify the subsequent analyses and facilitate the interpretation of the radiometric indices, to retain only the first two PCs which cumulatively explain almost 95%.

The composition of the first two rotated PCs can be derived from Figure 2.

Practically, all bands (410-1760 nm) from the visible (Vis, 380-760 nm) to the near infrared (NIR, 760-1500 nm) and part of the short-wave infrared (SWIR, 1500-3000 nm), except for a narrow range centered on about the 1400 nm band, weighed positively and significantly on PC1 (Figure 2). This can then be interpreted as an indicator of the average diffuse reflectance of the soil in the indicated radiometric range and of its color. Furthermore, remembering that absorbance is inversely proportional to reflectance, the PC1 can also be assumed to be an inverse indicator of iron oxide content [65].

On the PC2, mainly the bands from 1760 to 2500 nm weigh positively and significantly (Figure 2), comprising the major absorption bands for OC (primarily 1772 nm; 1871 nm; and to a lesser extent 660 nm; 2070 nm; 2177 nm; 2246 nm; 2351 nm; 2483 nm) and for clay (mainly 2201 nm and to a lesser extent 1877 nm; 1904 nm; 2177 nm; 2192 nm; 2220 nm; 2492 nm) [65]. The PC2 can therefore be interpreted as an inverse indicator of SOC and clay contents.

3.1. Raster data

Table 3 shows the basic statistics of elevation and topographic attributes obtained from LiDAR data and calculated at the nodes of a one-meter mesh grid. As can be seen, the study area with a maximum height difference of about 300 m is characterized by a high variability in its topography and consequently a high local relief. There are some very steep areas with slopes of about 72 degrees, which are practically inaccessible, where it was impossible to take soil samples. LS and SPI also show very high outliers and there are areas with negative curvature which indicate a surface is concave, and others with positive curvature (convex surface) with a prevalence of those with negative curvature (negative skewness).

The high variability of the topographic attributes can be explained by the changing morphology of the catchment under study (Figure 1) with an increasing elevation gradient from west to east and southeast. Furthermore, the area is characterized by a dense drainage network with numerous tributaries of the main stream, which are affect by of deep and narrow incisions.

In order to apply the KED to the sand content data, all auxiliary variables had to have been calculated or estimated at the nodes of the 1 m DEM. Therefore clay, SOC, PC1 and PC2 kriging estimates were added to the topographic attributes.

An isotropic model consisting of a nugget effect and an exponential model with a range of 394.5 m and therefore an effective range of about 1185 m (approximately the maximum size of the basin in direction 132° from the north) was fitted to the experimental variogram of the clay content (Figure 3).

This denotes a great continuity in the variation of clay content as is evident from its kriging map (Figure 4a) in which the highest contents occurred in the lowest parts whereas the lowest ones in the highest parts in the East and Southeast.

Different is the behavior shown by SOC, which was previously transformed to a Gaussian variable as mentioned before. An isotropic model including a nugget effect and a spherical model with a range of 137 m was fitted to the experimental variogram of its Gaussian transform (Figure 5).

The variability of SOC is thus characterized by less continuity than that of clay with spatial autocorrelation over a shorter range (137 m). The implication of the above is also evident in the back transformed map of SOC's kriging estimates (Figure 4b). It is characterized by a pattern of numerous small areas of limited size. However, two macro-areas can still be distinguished: one above the median line (132° N) of the basin with higher values and the part south of this line with lower values.

An isotropic LMC consisting of three basic structures: a nugget effect, a spherical model with a short-range of 159 m and a spherical model with a long-range of 716 m was fitted to the set of direct and cross experimental variograms of PC1 and PC2 (Figure 6).

As can be deduced from Figure 6, both in PC1 and PC2, the largest spatial component is unfortunately represented by the discontinuity at the origin of the graphs (intercept). That is called nugget effect and results from insufficient density of sampling with large unsampled areas. It is also shown that in the cross-variogram, the model reaches zero within a distance of about 150 m (very low sill). This means that the extraction and rotation of the PCs produced not only statistical but also quite spatial independence of the two components at a relatively short distance. This is also verified by the high distance from the intrinsic correlation curve (dashed line), which represents the maximum possible spatial correlation between the two variables (Wackernagel, 2003).

Figure 7a shows the map of PC1, which has a somewhat reverse pattern to that of SOC. This is consistent with the previously given interpretation of PC1 as a sort of spatial indicator of average soil reflectivity. It is known in fact that soils with a higher organic matter content are generally darker and therefore less reflective [65].

Figure 7b shows the map of PC2 with a large zone of higher values in the north. Interpretation becomes more difficult here, as the area at lower elevation (Figure 1) includes both zones with higher clay contents (north-west corner) (less reflective) and zones (south and southeast) with low clay and organic matter contents (more reflective). Actually, the diffuse reflection by soil in the range 1700-2500 nm depends on the presence of other components beyond the texture and organic carbon. Furthermore, selective absorption also depends on the particular composition of the clay and organic matter [66].

3.2. Coregionalization data set

The correlation matrix between sand content and the variables assumed as potential auxiliary variables is reported in Table 4. As can be seen, apart from the expected high negative correlation with clay, and the smaller positive correlation with elevation, the other correlations are rather low. It was therefore decided to exclude SPI and curvature as possible auxiliary variables, nevertheless to retain SOC and TWI to account for possible local influences of organic carbon and water content on estimation of sand content, because, generally, the different contents of sand are affect the soil moisture and the rate of decomposition and stabilization of the organic carbon [67]. On the other hand, spatial correlations might also be nonlinear and therefore not evaluated by Pearson's correlation coefficient.

3.3. Kriging with external drift

3.3.1. Trend estimation

For the trend calculation, various combinations (T=17) of internal (a linear function of x, y coordinates) and external drift (different linear combinations of the ten auxiliary variables selected: 1. clay; 2. SOC; 3. PC1; 4. PC2; 5. elevation; 6. slope; 7. aspect; 8. LS; 9. TRI; 10. TWI) were compared. For the internal drift, it was preferred to consider only a linear function of the coordinates (deterministic type variation), treating most of the variation in sand content due to spatial coordinates as stochastic and then described by the generalized covariance function. As can be seen from the examination of Table 5a in which the various trend models are sorted according to mean rank, the best model with the lowest mean rank and also the lowest mean square error was the one consisting of a constant coefficient plus clay content as auxiliary variable.

The automatic structure identification is made in two steps (Table 5a and b). In Table 5a, each line indicates the performance of an order k. The results based on polynomials fitted are reported under the heading Ring 1 and Ring 2. They correspond to the separation of data in two subset (ring 1 and 2) to improve the performance of jack-knife procedure: data in the ring 1 are used for the estimation of points of rings 2, and vice versa [39]. The second steps (Table 5b) is the identification of the generalized covariance, which is the equivalent of the variogram in the non-stationary case.

In the results, it is surprising is that the trend model with the inclusion of the first topographical attribute (elevation), although significantly correlated with sand content (Table 4) occupies the 5th position with lower mean absolute error but significantly higher mean square error and mean rank. The inclusion of the other topographic attributes worsened even further the performance of the model.

The interpretation of these results might be that the particle size composition of the soil in this catchment is essentially related to other intrinsic properties, while is little sensitive to relief features.

3.3.2. Generalized Covariance Function (GCf) identification

With regard to the estimation of the generalized covariance function, the model that realized the jackknife estimator closest to 1 is the one consisting of the nugget effect plus the linear component with scale equal to 200 m (Table 5b). The spatially uncorrelated component is approximately one order of magnitude larger than the spatially correlated one.

Synthesizing the two previous results on trend estimation and identification of GCf, it can therefore be said that according to the KED approach, all intrinsic spatial variation is stochastic with a prevalence of the spatially uncorrelated component (nugget effect).

In consideration of the fact that measurements of texture and in general of all soil properties are demanding of both man-hours and costs, it was decided to assess how much the sand model degraded by considering only topographic attributes as auxiliary variables. These are generally more available (less expensive) and at a much higher sampling density than observations of soil properties.

Eight trend models were then compared (Table 6), by subsequently adding all the other seven topographic attributes (TRI, aspect, slope, LS; curvature, SPI, TWI) as auxiliary variables to the elevation.

The best model in terms of both minimum mean rank and minimum mean square error is the one that uses elevation as the only auxiliary variable (Table 6a). This is indeed the only topographic attribute that showed a significant correlation with sand content, while the correlations with the other attributes were rather low, probably due to the uncertainty of their estimates and the high morphological variability due to denudation processes (landslides and water erosion) that affected the site [40,68,69]. These processes caused a reworking and a spatial redistribution of surficial soil.

As far as the generalized covariance function is concerned, it consisted again of two structures, nugget effect and linear component with a 200-m scale (Table 6). However, there was a marked prevalence of the spatially uncorrelated component, which was well about 2 orders of magnitude higher than the structured component. The exclusion of clay as an auxiliary variable in the trend therefore also resulted in an increase of randomness regarding the stochastic component of the sand estimation variance.

Assuming ordinary kriging as the univariate reference approach, an isotropic model containing three basic structures was fitted to the experimental variogram of the sand content: a nugget effect, a spherical model with a range of 150.99 m and a spherical model with a range of 513 m (Figure 8).

The variogram model appears to be well structured, upper bounded and with the spatial component at shorter range approximately twice as large as the uncorrelated spatial component and the structured component at longer range.

3.4. Comparison among the three approaches

The results of the cross-validation of the three models for estimating sand content are shown in Table 7: Model1, whose trend contains only the clay content as auxiliary variable; Model2 in which the trend has only the topographic attribute of elevation, and Model 3 with no trend.

As expected, Model1 was the best in terms of observation-estimate correlation, lowest bias, best accuracy and lowest correlation between standardized errors and estimates. This result was already predicted by comparing the ten trend models with different combinations of auxiliary variables including both soil attributes and topographic attributes (Table 5 and Table 6). However, Table 7 allows to assess the actual advantage obtained by supplementing the sand content observations with elevation, obtained from high spatial resolution LiDAR data, as the auxiliary variable. Model2 compared to Model3 has a higher observation-estimate correlation albeit small, a lower bias and a lower standardized error-estimate correlation. The reason for this rather modest result after the inclusion of elevation as an auxiliary variable is probably due to the particular contextual characteristics. However, this should not dissuade from continuing to investigate the real advantages of using high spatial resolution LiDAR data to improve the prediction of soil attributes. Soil measurements are generally expensive and the collection of soil samples at some locations may be difficult or even impossible as occurred in this case due to the extremely steep slope.

The previous results made it possible to compare the three models from a statistical point of view, essentially in terms of smoothing effect, bias, and accuracy. The purpose is now to continue the comparison in terms of mapping or graphical restitution of the sand content estimates and their corresponding uncertainties.

Maps of the sand content estimates for the three models are shown in Figure 9.

It is firstly to be said that the three maps are coherent in reporting the main structures of spatial dependence. However, it is possible to highlight clear differences among the maps. The one for Model1 (clay as the auxiliary variable, Figure 9a) shows a high degree of spatial continuity and, and mostly reproduces the inverse relationship with the clay due to the strong negative correlation between the two variables. The Model2 map (trend with elevation as the auxiliary variable, Figure 9b), although quite similar to the previous one, shows greater short-range variability due to the finer spatial resolution of elevation compared to the sampling scale of clay. However, it is worth remembering that clay was interpolated with ordinary kriging at the DEM grid nodes to be able to apply KED. Finally, the Model3 map (no trend, Figure 9c) shows less spatial continuity compared with the one of the one of Model1 and is characterized by more numerous micro-structures of limited extent, probably due to the random nature of the spatial variation of sand content. In Figure 10 the uncertainty of the estimates, as expressed by the standard deviation of estimation, is compared for the three models.

Being OK and KED linear interpolators, the estimation standard deviation depends on the sampling scheme and the mathematical model of the variogram for Model3, whereas the mathematical model of the generalized covariance function for Model1 and 2. Consistently, the 3 models show the lowest values in the central part, where sampling is densest, and the highest values at the edges of the area. However, compared to Model3 (Figure 10c), which faithfully reproduces the sampling pattern with minima at the sampling locations, Model1 (Figure 10a) and 2 (Figure 10b) reveal more spatial continuity with a greater accentuation of short-range variability for Model2.

What these results highlight is that the choice of auxiliary variables can have a significant impact on mapping of estimation and its uncertainty of the variable of primary interest. Even if the same macro spatial structures are reproduced by the various models, there can be appreciable local differences and only a careful validation process can help in choosing the optimal model.

4. Conclusions

The main objective of this work was to evaluate whether combining observations of soil properties with LiDAR-derived topographic attributes can be considered an efficient and straightforward approach for mapping soil texture (sand content) both at fine resolution and at large spatial scale. Through LIDAR, we have topographical knowledge of a vast area compared to manual measurements, and the results highlighted the importance of applying this tool for accurate soil measurements. Various linear trend models using different auxiliary variables were compared. Statistically, the best one used clay content as the auxiliary variable. Actually, the improvement obtained by using only LiDAR data (elevation) compared to the univariate (no trend) model was rather marginal.

Nevertheless, the terrain morphology detected by Lidar can affect the soil moisture and nutrients available and play an important role in forest dynamics, and growth rates.

However, it is widely recognized that there are various advantages, including economic ones, that would result from the use of LiDAR data in predicting soil lithology (Prince et al., 2020), and that LiDAR elevation data are now becoming the new standard for the production of high-resolution DEMs worldwide. What was said before should therefore lead to the development of new and novel methods. We believe that with the increased availability of data, an improvement in digital soil mapping (DSM) can be achieved by the application of Machine Learning (ML) techniques to the development of trend models, leaving geostatistics to assess and model the stochastic component of spatial variation. Anyway, ML models should be more interpretable and explainable [70] and not just more accurate, so that important scientific insights can be extracted from soil data and such models can serve as an effective source of information and knowledge.