An Approach for Spatial Statistical Modelling Remote Sensing Data of Land Cover by Fusing Data of Different Types

1. Introduction

The use of remote sensing in mapping and assessment of soil resources is constantly increasing, as is the possibility of achieving ever greater spatial resolution. There is a rich literature on the applications of remote sensing in soil ranging from general assessment of soil resources and its main properties for different purposes [1,2,3,4,5]. Interest covers the application of remote sensing in agriculture [6] and especially in precision agriculture [6,7].

The application of remote sensing requires the use of explicit spatial statistical techniques whose appropriateness needs to be proven even though they are widely used, and any shortcomings of them might require new spatial statistical modelling approaches. Geostatistics [8,9] is a suitable approach to spatial modelling remote sensing data and its use is almost missing although the link between remote sensing and geostatistical methods was established in the late 1980s [10,11,12] to quantify image spatial structure [13,14,15], integrating spatial statistics and remote sensing [12], image classification [16]. Remote Sensing is a great source of spatial (and space–time) environmental data and geostatistical methods can be to right approach to many issues arising from their analysis.

The set of measurements of reflectance obtained from a scene is associated with a spatial position in the scene and, therefore, as a regionalized variable function of spatial position and neighboring points in the scene are consequently more similar spectrally than distant points [17,18]. That is called spatial autocorrelation, and it disappears with distance. The mean distance at which this spatial autocorrelation disappears is related to the size, spacing, and shape of the objects in the scene [17,18]. The spatial structure resulting also depends on the spatial resolution of the sensor with respect to the objects imaged. Therefore, to use data from different remote sensors, it is necessary to consider three basic concepts related to a remote sensing image: support, spatial resolution and pixel [19,20]. The support is an important geostatistical concept representing the space on which a measurement or observation is made, and it is characterized by geometrical shape, size, and spatial orientation [19,21]. It should not be confused with spatial resolution [22] ,which represents the average sampling distance in sparse sampling, whereas support is the size of the grid mesh or pixel [20].

Using remote sensing data from different satellites or sources it is required to define a data fusion of heterogeneous spatial data for prediction which considers the different supports. The support of measurement in RS is commonly represented by a 2D Gaussian function (or similar function) referred to as the ‘point spread function’ (PSF), with its tails extending far beyond the limits of a pixel [23]. Little spatial statistics research in remote sensing has been devoted to the impact of support and PSF shape on the prediction and its uncertainty. It is worth underlining that the main target of spatial statistics in RS is to extract information from remotely sensed data. Information can be defined as the difference between data, and therefore in a single waveband remotely sensed image information consists in the differences between neighboring pixels [24]. This concept is sufficiently general [23] and allows the variance of a Gaussian data distribution to be used as a diagnostic tool because it describes what is expected on average. In the case of spatial data, its extension to the concept of a spatial covariance function is particularly useful. It is possible to calculate the empirical spatial covariance and to fit a model to it by using an appropriate method of inference. Common permissible (authorized) functions include the Matérn family of models, such as the popular exponential covariance model. The exponential model has two parameters, the so-called sill variance and range (or practical range). To a certain extent, the sill variance of the exponential model can be thought of as the spatial equivalent of the sample variance. The range on the other hand has no equivalent in the data distribution in traditional statistics. However, the range parameter is very useful because it tells us a lot about the potential information content of an image. The range in fact informs about the scale(s) of spatial variation present in the image [19,25,26,27,28,29,30] and represents the upper limit of spatial correlation between observations. Also, the shape of the spatial covariance function is quite informative. For example, the exponential model, which is asymptotic towards the sill variance, can be assumed representing a set of scales of variation, each with its own information-to-redundancy ratio. Yet in recent times there has been a lack of interest in RS in dealing with spatial statistical functions and in particular spatial covariance functions, neglecting all the potential information about the nature of the variable of interest that they can provide us with geostatistics [24].

Remote sensing images are essentially tables of data commonly represented by digital numbers. However, when the main interest is extracting information from these data, it is essential to define statistical models that can characterize the spatial variation and thus predict a variable of interest at points where sampling or observation has not been done. Indeed, various techniques exist to do this, such as geostatistics and its recent developments in multi-point geostatistics [31,32], spatially explicit regression models that introduce spatial dependence between pixels [33] and more recently Machine Learning techniques [34,35]. This paper is mostly focused on the application of geostatistics to RS data.

Geostatistical model can in fact be used for (i) handling missing data (e.g., generated by the presence of clouds or cloud shadow [36,37], or by a sensor's operation failures [38]; (ii) upscaling from a set of pixels to a bigger (regular or irregular) block [39]; (iii) downscaling from a coarse pixel to a finer spatial resolution [40,41] ; and (iv) fusion of images of a certain spatial resolution with other images of a different spatial resolution [42]. Both upscaling and downscaling involve the problem of change of support which initially was addressed to eliminate the strictures of the pixel in RS. Pixel is considered as the average reflectance of light from a limited support, therefore decomposing the original support into a set of finer supports is indeed a fairly common problem, widely reported but not easy. This can be done if additional information is available in a multivariate approach, which means fusing multi-source data with different support. Multivariate geostatistics provides a way to solve the downscaling problem in RS, which is particularly crucial in the presence of mixed pixels and in the application of site-specific management in agriculture. A key concept in geostatistical analysis is the spatial covariance function, or its related function, the variogram, which characterizes the spatial variation and parameterizes a random function (RF). Under the condition of second-order stationarity, remotely sensed images can be treated as a realization from a RF and then it is possible to apply geostatistical techniques in fusing them with different kinds of data considering their respective supports. Indeed, there is currently a scarcity of attempts in RS that apply downscaling or upscaling considering change of support with its implications on the variance of estimation or better on the spatial dependence function [43].

The choice of spatial resolution is an extremely crucial problem in RS, because the pixel size could be too large relative to the scales of variation of the process which remains the researcher’s main interest rather than the image itself, i.e. an understanding of reality with its mechanisms and functions.

Geostatistics offers the possibility to fit mathematical-statistical models defined on punctual support to the data on finite support such that other mappings can be generated readily (e.g., to coarser support, to a finer support or to a support with complex geometry such as Management Zones (MZ) in Precision Agriculture). It is quite important, because one of the most significant challenges in environmental data science today is the inability to allow datasets, that were obtained with different supports and with other measurement characteristics, to speak to each other [44]. Interoperability among different data sets requires that data on one support can be changed to data on another support reciprocally. This demand can be addressed by an appropriate multi-source data fusion approach although, unfortunately, such a vision is not common, especially in RS. Nevertheless, it is a crucial issue that needs to be addressed because environmental data comes in a wide variety of shapes and sizes (property definitions, sampling frameworks, measurement characteristics, error characteristics). Dealing jointly with this multiplicity of heterogeneous data properly and drawing as much information as possible is not a trivial problem and should appeal to spatial statisticians.

Moreover, geostatistical RF characterizations can be helpful in revealing the scales of spatial variation in spatial data and the potential information content or conversely, the redundancy in data. They can also provide some insights into the data which can have different uses, such as stating whether a particular sampling or method is appropriate for a task: using the right model for the right task by considering one's own purposes carefully.

The objective of this paper is to propose an approach for spatial statistical modelling in remote sensing of land cover by fusing data of different types. Agricultural land use composed of different land covers as soil and trees arranged spatial patterns was considered. The results of this study may provide useful insights and foster discussion on conceptually important topics in relation to the development and application of spatial statistics in RS.

2. Materials and Methods

2.1. Study Area and Soil Data

The study area (Figure 1) is an olive grove (’Ogliarola Salentina’ cultivar) located in south-eastern Italy (Brindisi, 40°31′10.57′′ N; 17°39′42.86′′ E). The field is flat and gravel free and soil is classified as Calcic Haploxeralf [45]. “Calcic” regards a calcic horizon with its upper boundary within 1 m of the mineral soil surface, and “Haploxeralf” involves a loamy granulometry and a relatively thin argillic or kandic horizon.

Composite topsoil samples (0–0.20 m) were collected in July 2019 at 61 locations around each olive trees. Three soil samples were collected around the trees (at about 1 m) and were mixed. More details on the study area can be found in Belmonte et al. [46].

Soil samples were analyzed in the laboratory for different properties, but here only textural fractions were included in the study. The different soil textural fractions (%) were classified according to the particle size limits used in the USDA scheme [47]: total sand (2.0-0.05 mm), silt (0.05–0.002 mm), and clay (<0.002 mm). Clay, silt and coarse sand (2.0–0.2 mm) were determined by pipette method [48]. Total sand represents the sum of coarse and fine sand (0.2-0.05 mm); the latter was calculated as the complement to 100 of the sum of the measured soil particles fractions.

2.2. Hyperspectral Analysis

VIS-NIR-SWIR reflectance of topsoil (0–0.20 m) was measured in laboratory on a fraction of the soil samples collected for granulometry analysis with an ASD FieldSpec IV spectroradiometer (Analytical Spectral Devices Inc., Boulder, Colorado, USA) over a spectral range of 350–2500 nm [49]. The spectra were averaged over 10 nm so reducing the number of wavelengths from 2151 to 215 to smooth them and avoid the risk of over-fitting [50].

The following spectral pre-processing methods were carried out to reduce the multiplicative and additive effects, preserving the useful information for prediction models [51]: 1) multiple scatter correction (MSC) and standard normal variate (SNV); 2) smoothing/denoising with Savitzky–Golay polynomials (SG). Multiplicative scatter correction (MSC) removes additive and multiplicative components [51,52], mostly due to scatter effect so reducing the number of components needed in the calibration model and improving the linearity [51]. The Standard Normal Variate (SNV) method carries out a normalization of the spectra that consists in subtracting each spectrum by its own mean and dividing it by its own standard deviation [52]. After SNV, each spectrum will have a mean of 0 and a standard deviation of 1: the result is like that of MSC. SNV is frequently used to compensate for changes in surface roughness of the material [53].

Savitzky–Golay (SG) algorithm [54] reduces the random noise of the measurements. For a given spectrum measured at N points and a filter of width w, it calculates a polynomial fit of any order in each filter window as the filter is moved across the signal. The size of the filter consists of (2n + 1) points, where n is the half-width of the smoothing window (w). The filter calculation is complete when the filter window moves the signal one-at-a-time. A window size (w) of 11 (w = 2n+ 1) and polynomial order of two were applied [55].

Spectral data transformation and pre-processing were performed with ParLeS software [56].

Different combinations of the pre-processing methods were compared to the use of raw data of reflectance by assessing the performance of prediction models [57].

2.3. Prediction Models

To estimate soil granulometric fractions, Partial Least Squares Regression (PLSR) [58] was applied. PLSR is based on decomposition of two sets of variables: the matrix X of predictors (matrix n x NX, where n is the number of observations (61) and NX is the number of spectral data at different wavelengths (215)), and the vector Y of the response variable (sand or clay or silt). PLSR selects successive orthogonal factors (or latent variables) that maximize the covariance between predictors (X, reflectance) and response variable (Y, sand or clay or silt) and explains most of the shared variation. For more details on PLSR method. Leave one out cross validation was performed to determine the number of latent variables that minimized the predicted residual sum of squares (PRESS). Each calibration set is created by taking all the samples except one, the validation set being the sample left out. The performance of PLSR prediction model was evaluated by means of three statistics: (1) the coefficient of determination R² of the regression: observation vs. prediction; (2) root mean square error of prediction (RMSEP), and 3) the ratio of performance to interquartile distance (RPIQ) [59] defined as follows:

$$RPIQ={\displaystyle \frac{IQ}{RMSEP}}$$

(1)

where:

$$RMSEP=\sqrt{{\displaystyle \frac{{\displaystyle \sum _{i=1}^n\left(y_i-y_i^{\text{'}}\right)}}{n}}}$$

(2)

y_i are the observed values, y_i’ the predicted values, and n is the number of samples with i=1, 2,…, n.

IQ is the difference between the third and first quartiles of y: $IQ=Q_3-Q_1$

RMSEP has the same unit of measurement as the dependent variable therefore RPIQ is unitless and an inverse measure of error related to population spread, the greater than 1.5-2 the better. R² measures the precision of each model and the higher the value the better the fit. After cross-validation, analysis of residuals was performed with Shapiro-Wilk and Kolmogorov-Smirnov tests, to evaluate if the distribution was normal. When the probability level for both statistics was greater than 0.10, the distribution was assumed normal [55].

To identify the best wavelength ranges for prediction, the analysis was performed both on the full range and on specific ranges by analyzing 100 nm intervals at a time. The ranges, causing a model with the highest R², were chosen. Thus, for clay were selected the wavelengths between 360-440 nm and 1560-2500 nm, while for sand and silt 450-1350, 1430-1550 and 2060-2400 nm.

A total of 108 models were carried out, based on: 1) two types of data: reflectance and absorbance; 2) 9 combinations of the pre-processing methods on reflectance and absorbance data; 3) 3 granulometric fractions, clay, sand and silt; 4) 2 wavelength ranges: full (350-2500 nm) range and selected ranges as described above.

2.4. Remote Sensing Products

Reflectance unmanned aerial vehicle (UAV) acquisition data

The high spatial UAV image was acquired on July 2019 with a DJI Mavic Pro drone whose custom payload carried a multispectral sensor (Parrot Sequoia) that captures data (16 Megapixel) with four bands (see Table 1) at 70 m height and with a theoretical resulting Ground Sampling Distance (GSD) of 0.07 m/pixel. The mission flight plan is managed by the Mission Planner software which ensures 80% front and side overlap of the frames. The reconstruction of the scene is performed by version 4.4.12 of Pix4d Mapper Pro software by applying Structure from Motion (SfM) [60] techniques on a total of 1152 images acquired (288 frames for each multispectral band). A final orthophoto image for each band was generated from a quoted point cloud.

Reflectance Planet acquisition data

In recent years new commercial satellite platforms, such as PlanetScope, have become available. PlanetScope is operated by Planet Labs [61]. This platform is a CubeSats constellation that consists of over 180 satellites deployed in sun synchronous orbit with a daily revisit time, which offers more opportunities for cloud-free images. The first Planet satellites were launched in 2015. The first generation satellites have an optical multispectral sensor with four spectral bands at red, green, blue, and near-infrared wavelengths, with a spatial resolution of 3.7 m at nadir; while the new generation of PlanetScope satellites has the potential to produce imagery with eight spectral bands.

In this study the reflectance of PS2 product Plane, orthophoto image, with four bands (see Table 1) was used. This product is into unit16 data type, to represent the data into range [0, 1] it is necessary to divide by 10000, as specified into metadata.

The Planet image was downloaded from Planet website [62].

Reflectance Sentinel 2 acquisition data

The Sentinel 2 image was downloaded from Copernicus Space Data Ecosystem (CDSE) Portal [63]. The S2-MSI satellite data, as surface reflectance product, were used for this work. This satellite was launched in 2015 by the European Space Agency (ESA) with a variable spatial resolution from 10 to 60 m and thirteen spectral bands. In detail we worked with Level-2A data that provides atmospherically corrected bottom of the atmosphere reflectance images.

In this work all UAV bands, all Planet bands and only two Sentinel 2 bands (11 and 12 one), see Table 1, were used as auxiliary variables in the process of soil property prediction. All the images have been chosen so that their acquisition date was as close as possible to the samples date.

2.4.1. Filtering of Soil Data from Remote Sensing Product

−

UAV and Planet data

The processing focuses on retrieval soil information involved supervised classification to select thematic information. The process was using the following principal steps present in the commercial ENVI (ENvironment for Visualizing Images) 5.1 software (Harris Geospatial):

−

Training Samples

To identify thematic information, the set of ground cover types (i.e., information classes) into which the image would be segmented was first determined. Prototype pixels, or training samples, were then assigned to these selected classes. At the beginning, for this study, three classes were chosen based on different but sufficiently homogeneous regions of the same image: soil, canopy, and shadow. To account for class variability and minimize the influence of transition areas where mixed pixels may occur [64], enough training pixels (at least 150 pixels per class) were selected. This ensured a representative sample of each class.

−

Maximum likelihood

The training data were applied to estimate the parameters for the classifier algorithm. The decision rule applied was the maximum likelihood method, utilizing a Gaussian model. This process resulted in the generation of a thematic map, with the desired classes clearly delineated.

−

Morphological filter

To keep intact the geometric structure of image objects, a morphological filter based on the closing operator [65] was applied. This operator modifies the geometric features of the signal by removing small holes within polygons. The morphological convolution is performed on the signal using a structuring element, which is a set of simpler shapes and sizes [66].

−

Combine classes

This step is used to build a binary image merging two classes (canopy and shadow) into one class labeled ‘soil’, identified with 1 respect to class ‘not-soil’, identified with 0.

−

Conversion raster to vector

The binary images on each UAV and Planet satellite band [49] were converted from a labelled raster image into vector data format (such as ESRI Shapefile format), employing the tool (r.tovector) of the GRASS (Geographic Resources Analysis Support System) module integrated in the open-source software: QGIS (Quantum Geographic Information System) version 3.38.0 – Grenoble.

−

Soil information

To retrieval soil information the ‘raster extraction’ tool available in QGIS was applied. By it each raster image is clipped by mask layer formerly created into shapefile format.

For Sentinel-2 bands 11 and 12, soil retrieval is performed exclusively using a geostatistical method, as described in Section 3.2, due to their too coarse spatial resolution (20 m).

2.5. Geostatistical Multi-Source Data Fusion

Spatial data conversion between different spatial scales is a common challenge in environmental applications because of spatial mismatch between different data sources, such as various types of observation sensors in Remote and proximal sensing and /or sample collection methods. One of spatial support include the pixel size in satellite imagery. Spatial support is a very important concept in geostatistics because it determines the scale at which spatial variability of input variables can be analyzed and the one of the predictions.

Change of spatial support refers to changing from one type of spatial support to another during the spatial prediction. Using the terminology by [67], upscaling consists in aggregating point data into an area and downscaling in disaggregating an area data into points.

Accurate downscaled information is crucial for making decisions in precision agriculture or environmental management. Geostatistical spatial scaling for remote sensing downscaling was used by [40]. In the case of co-registered multivariate satellite sensor images, [68] proposed area-to-point downscaling cokriging, where the pixel size to be predicted is smaller than the one of the input images.

A practical problem of interest in remote sensing is to transform a coarse spatial resolution image into another fine spatial resolution image by introducing spatial ‘detail’ but in a way that is coherent with the spectral information of the original image. A particular solution is provided by area-to-point kriging [69], but in most cases more than one experimental variable may be represented on different supports and the support of the variable to be predicted may not be a point. This is often the case in remote sensing, where a satellite sensor may provide different images for several spectral wavebands with potentially different spatial resolutions (pixel sizes). In such cases, a problem of digital image analysis is to combine the images to increase the spatial resolution (i.e. decrease the pixel size) of the coarser spatial resolution images while preserving their spectral content. The objective may be twofold: improving the visual aspect of the image and facilitating image interpretation for better e.g. soil characterization. However, most of the existing algorithms for image fusion are aimed at obtaining a visually pleasant image but do not consider the spatial autocorrelation of the waveband and then its spatial structure either at coarse or fine resolution. Moreover, in order to jointly analyze (fuse) multi-source data, both sample and grid data such as RS images of different resolution , geostatistics uses a multivariate interpolator (cokriging) which has several advantages compared with traditional approaches in RS, already underlined by [70]: 1) provides a prediction which is unbiased and minimizes the prediction variance; 2) explicitly takes into account the support (pixel size in RS) of each variable via its variogram; 3) takes into account the correlation structure via the set of direct and cross variograms; 4) can explicitly take into account the point spread function (PSF) of the satellite sensor. PSF makes it possible to describe the response of an optically focused (image) system to a point source or object; and 5) provides an estimate of the uncertainty of the prediction.

To be able to use point cokriging applied to the entire multivariate data set, including the sample soil variables, satellite images (Planet and Sentinel2) and UAV images, it was necessary to reduce in advance all variables to the same point or reference support.

Assuming as the reference (point) support that of the soil samples (i.e. a 1m x 1m block), the satellite images were deconvolved to the same block through deconvolution by kriging, while the UAV images were convolved to the same block through block cokriging.

Below the whole approach used to fuse the multi-source data of our study is illustrated and the techniques used were briefly described. The interested reader may refer to the numerous bibliographical references given for further details.

2.5.1. Deconvolution of Images by Kriging

In image analysis using satellite or UAV images the available data are not on a point support but are convoluted by a weighting function on an extended area (pixel). This function depends on the underlying physical process and the characteristics of the sensor, and can be established from theoretical considerations, or can be measured with appropriate experiments [71]. Various approaches enable images to be deconvoluted. All these approaches are very common in the image analysis literature [72]. Deconvolution of images is generally treated by Fourier transformation. However, this approach is known to be not applicable when considering convoluted data with noise or when a regular grid of data is not available and therefore missing data have to be interpolated. It is for this reason that in this work another approach was used based on the process of estimating point data by deconvoluting kriging [71].

An image in Rⁿ Y(x) at the finer pixel x can be modified through a convolution by a weighting function $\stackrel{⌣}{p}(x)$, which is normalized, so that its integral over spatial coordinates is equal to 1, and $\stackrel{⌣}{p}(x)=\stackrel{⌣}{p}(-x)$. Experimentally, the observed data Z at coarser pixels x_α are given by:

$$Z\left(X_{\alpha }\right)=Y\left(x\right)\times \check{p}\left(x\right)$$

(3)

Therefore, knowing the weighting function $\stackrel{⌣}{p}(x)$, the problem is to estimate the image $Y(x)$ at each pixel x from the data $Z(x_{\alpha })$ known at coarser pixels $x_{\alpha }$. This can be done by estimation with kriging. The weighing or convolution functions depend on physical conditions, however the most common are: uniform, exponential, gaussian and sinc functions [71]. This approach has been developed in the frame of geostatistics [73] and considers information on the spatial structure of the underlying phenomena. It essentially consists in the estimation of the unknown data by a local linear regression. Applications to RS image analysis have already been proposed [68,74] but they are not very common, especially in recent years when Machine Learning techniques are preferred. Here we consider the common case of a single datum per pixel x, such as the recorded reflectance $Z\left(x\right)$ at a given wavelength. The variables $Y\left(x\right)$ and $Z\left(x\right)$ are assumed to have the same mean and to be realizations of an intrinsic random function, whose increments $Z\left(x+h\right)-Z\left(x\right)$ are assumed as a second order stationary random function (Webster and Oliver, 2006). In geostatistics variogram is the cornerstone of the prediction, which for the random function Z is defined as follows:

$${\gamma }_z\left(x\right)={\displaystyle \frac{1}{2}}E\left[{\left(Z\left(x+h\right)-Z\left(x\right)\right)}^2\right]$$

(4)

where E is mathematical expectation. To estimate the underlying variable $Y$ from the observed variable $Z$, an estimator $Y^{*}\left(x\right)$, i.e. a linear combination of data in a neighborhood of x (pixels $x_{\alpha }$), was used:

$$Y^{*}\left(x\right)=\sum _{\alpha }{\lambda }^{\alpha }Z\left(x_{\alpha }\right)$$

(5)

The weights in equation (${\lambda }^{\alpha }$) are the unique solution of the linear kriging system, obtained from unbiases and minimal error variance conditions for the estimate:

$$\begin{gathered} \sum _{\alpha }{\lambda }^{\alpha }{\gamma }_Z\left(x_{\alpha }-x_{\beta }\right)+\mu ={\gamma }_{YZ}\left(x-x_{\beta }\right) \\ \sum _{\alpha }{\lambda }^{\alpha }=\left\{1\right\} \end{gathered}$$

(6)

where ${\gamma }_{YZ}$ is a cross-variogram:

$${\gamma }_{YZ}\left(h\right)={\displaystyle \frac{1}{2}}E\left[\left(Y\left(x+h\right)-Y\left(x\right)\right)\left(Z\left(x+h\right)-Z\left(x\right)\right)\right]$$

(7)

and the variance of the estimate is given by $\mu +\sum _{\alpha }{\lambda }^{\alpha }{\gamma }_{YZ}\left(x_{\alpha }-x\right)$, where$\mu$ is Lagrange multiplier.

The practical implementation of deconvolution by kriging requires one to solve the cokriging system (Equation 6) after performing a so-called structural analysis to estimate the underlying variogram ${\gamma }_Y$. The procedure consists of different steps:

1.

Calculation of the experimental variogram ${\gamma }_Z$ from data Z. This variogram should show a nugget effect in the presence of noise, followed by very regular behavior for small distances h, due to the convolution function $\check{p}\left(x\right)$;

2.

Choice of a variogram model for ${\gamma }_Y$, possibly including different scales. Theoretical models are inferred from a priori knowledge of the structure and from the behavior of the experimental variogram ${\gamma }_Z$ and setting of supports of the variables $Y$ and $Z$

3.

Calculation (mostly by numerical interpolation) of a theoretical model of ${\gamma }_Z$ from ${\gamma }_Y$, and $\check{p}\left(x\right)$.It is worth pointing out that if the knowledge of $\check{p}\left(x\right)$ is necessary, otherwise no correct deconvolution can be expected;

4.

Comparison between experimental and the model variogram${\gamma }_Z$. If necessary, perform correction of the model and repeat steps 2-3-4;

5.

Choice of the optimal neighborhood and numerical solution of the cokriging system (Equation 6);

6. 

Restoration of the image from the system of weights and calculation of the quality of the deconvolution. This procedure is implemented in the software package ISATIS Release 18.5 (www.geovariances.com) and was developed by [71].

2.5.2. Coregionalization Data Set

The multivariate data set of the case study included both sample variables (point vectorial variables) and gridded areal variables (raster variables). The target variables or variables of interest were the soil texture variables (clay, silt and sand). However, only clay and sand were subjected to the multivariate geostatistical analysis, while silt was obtained as a complement to 100 of the sum of the other two texture components. This was made necessary by the fact that not all texture components were independently measured, as above described. Therefore, the three variables were linearly dependent, thus making the resolution of the cokriging system impossible because of the singularity of the coefficient matrix. The auxiliary variables were both sample (clay and silt predictions using hyperspectral spectra) and raster (all the RS bands of UAV and Planet), while for Sentinel 2 only the R11 and R12 bands were used as they were [49].

To apply data fusion according to the techniques of multivariate geostatistics, it is necessary, after having reduced all study variables to the same support (1 m x 1 m), to construct a coregionalization data set to have enough pairs of observations to calculate cross-variograms. These calculations require using the migration procedure: move the closest pixel of the grid variables to each soil sample, if the distance does not exceed a threshold, as small as possible (1 m here, i.e. the size of one pixel), which allows obtaining sufficient observation pairs. Sixty-one was the number of observations for the multivariate coregionalization data set.

2.5.2.1. Linear Model of Coregionalization (LMC)

The multivariate spatial correlations were modelled using the linear model of coregionalization (LMC) [75], which assumes all the n variables as produced by the same independent physical processes, acting on N_S different spatial scales. Therefore, the n(n + 1)/2 direct and cross-variograms (a direct variogram for each variable and a cross-variogram for the pair of variables) are modelled as a linear combination of the same series of N_S basic functions of variogram $g^u\left(h\right)$:

$${\gamma }_{ij}\left(h\right)=\sum _{u=1}^{N_s}{b}_{ij}^u{g}^u\left(h\right) i,j=1,\dots , n$$

(8)

where ${\gamma }_{ij}\left(h\right)$ is the model of cross-variogram between variables Z_i and Z_j; the apex u indicates the spatial scale; h the lag, and $b_{ij}^u$ the coefficients of the function $g^u\left(h\right)$, namely the partial sills of the variograms at each scale u [76]. Putting these coefficients together for each scale u, the corresponding coregionalization matrix is obtained, which can be interpreted as a variance-covariance matrix at the scale u.

2.5.2.2. Multicollocated Cokriging with Images Deconvoluted by Kriging

To fuse the sample variables with gridded variables for the final prediction of soil properties (clay and sand contents), the variant multi-collocated point cokriging was preferred rather than ordinary cokriging [77,78,79]when modelling sparse primary variables with secondary information much more finely sampled. As using all secondary exhaustive information contained within the interpolation neighborhood may lead to an intractable solution, due to too much information (big data), only the secondary data at the target location (where the estimation of the primary variable is required) and at points where the primary variable is available are employed for interpolation. This simplified version has generally produced more reliable and stable results than more known collocated cokriging, that uses only the secondary variable collocated with the position where the primary variable is sought [80].

Since the secondary datum collocated with the target position tends to screen out the influence of other secondary data placed at a larger distance, the loss of accuracy in considering the reduced neighborhood is indeed mostly negligible. However, the loss of information is well compensated for by the considerable reduction in computer time requirements.

To assess the influence of the auxiliary variables on the prediction of clay and sand content, the previous approach was compared with the univariate kriging prediction of the same variables. To perform this comparison, cross-validation [81], whereby one observation is temporally removed from the data set and re-estimated from the remaining data.

For the comparison the three statistics suggested by [82] were used:

$$ME={\displaystyle \frac{1}{N}}\sum _{i=1}^N\left(Z\left(x_i\right)-Z^{*}\left(x_i\right)\right)$$

(9)

$$RMSSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{I=1}^N\left({\displaystyle \frac{Z\left(x_i\right)-Z^{*}\left(x_i\right)}{{\sigma }_i}}\right)}$$

(10)

$$RMSE=\sqrt{{\displaystyle \frac{1}{N}}\sum _{I=1}^N{\left( Z\left(x_i\right)-Z^{*}\left(x_i\right) \right)}^2}$$

(11)

where Z($x_i$) is the estimate of soil property at site $x_i$_, $Z^{*}\left(x_i\right)$ is the observation at the same site, ${\sigma }_i$is the root mean squared prediction error of Y(x)₎. Mean error (ME) was used to assess the unbiasedness of the predictor and its optimal value should be approximately zero; Root Mean Squared Standardized Error (RMSSE) assesses the accuracy of mean squared prediction error and should be approximately 1; Root Mean Squared Error (RMSE) checks the goodness of prediction, and the approach with smaller value should be preferred because this means that predictions are closer to the observed values.

To make the comparison between the univariate and multi-source data fusion approaches more complete, the coefficient of determination (R²) between measured and estimated values and the correlation coefficient between the standardized error and the estimate (ρ) were also calculated. The latter, if the spatial model is indeed free of systematic errors, should be statistically not different from zero.

3. Results

3.1. Hyperspectral Data

Table 2 reports the results obtained with PLS analysis regarding the applied pre-treatment methods with a R²≥0.5, for the prediction of granulometric contents. All the selected models for clay and sand show a good performance with RIPQ > 2. The raw data of clay reflectance showed a model with a R²<0.5 using 'whole spectrum from 350-2500 nm (full range). Differently, when the wavelengths were selected from 360-440, 1560-2500 nm ranges, the performance improved sharply with data pretreated with SG. The first range comprises the absorbance peak of goethite (at around 420 nm), an iron oxide [83]; the second range includes the absorption peak of montmorillonite (at around 1900 nm), a component of clay, and that of kaolinite (at around 2200 nm), which is found to be strongly associated with clay [84].

The best model for sand reflectance data was obtained by considering the full range after pretreatment with SG (Table 2). A comparable R² (0.56) occurred when using ranges between 450-1350, 1430-1550 nm and 2060-2400 nm, where reflectance tended to increase. However, in the latter case the residual distribution can be assumed normal. The above results should thus suggest that the prediction of clay is mostly affected by those wavelengths corresponding to the maximum of absorbance and then to the minimum of reflectance. Differently for sand whose prediction seems to be more affected by those wavelength intervals in which reflectance tends to increase. However, it should be emphasized that in the interval (1430-1550 nm), not included in the clay model, the absorption peak of illite, a component of clay, occurs. This result might probably be due to the chemical affinity of this clay component, consisting mainly of silicates, with sand [85]. As regards the range between 2060-2400, which is shared by the clay and sand models, an explanation for its impact on clay prediction may probably be the occurrence of kaolinite peak absorption [84], while for sand the issue merits further investigation.

For silt no models reached R² >0.5 and the best model was obtained with the same ranges of the sand model. This probably depends on the chemical and mineralogical composition of the silt. The above results underline the importance of knowing the mineralogy of the soil before estimating the granulometric components with spectroscopy.

3.2. Geostatistical Processing

3.2.1. Deconvolution of Planet and Sentinel2 Data

The iterative approach described in Materials and Methods was implemented by fitting a nested spherical variogram model to the original areal satellite data and using a Gaussian function as a weight function according to the PSF information of the sensors on board the Planet and Sentinel 2 satellites. Table 3 shows the results of the structural analysis (variogram model determination) for each band of Planet and for bands B11 and B12 of Sentinel 2.

As can be seen from Table 3, a support ratio of 0.33 for Planet and 0.05 for Sentinel 2‘s B11 and B12 was used, roughly corresponding to the ratio of the reference pixel size (1 m) to that of Planet (3.17 m) and B11 and B12 (20 m).

This ratio is very crucial in deriving the variogram model of the deconvolved variable, as it determines the nugget ratio of the two (raw and deconvolved) variograms. Lacking direct knowledge of the deconvolved variable, this ratio is subject to great uncertainty and the estimation of the variogram model is not uniquely determined. The quality of the deconvoluted variogram model can only be estimated from a careful analysis of the deconvolved image and the possible existence of bias or artefacts with respect to the raw image. With regard to the variogram model for the Planet bands, it consists of 3 structures: a nugget, a spherical model with a short range varying between 8.48 m and 20 m, and a spherical model with a long range varying between 33.97 m and 52 m. For the B11 and B12 bands, given their close correlation, a multivariate approach was used with an estimated LMC including only one structure, that is, an exponential model with a range of 57 m. The greater simplicity of the variogram model of the Sentinel2 bands compared to those of Planet is probably due to the greater pixel coarseness. However, it should be noted that the longer-range structure is quite consistent between the two types of satellites. This leads us to think that spatial association structures of an average size of approximately 50 m do indeed exist in the study area, probably related to intrinsic variability of soil as surveyed also with different sensors [49,86].

Figure 2a shows on the left the maps of the Planet reflectance band after the filtering of the soil data because of the screen operated by the olive tree canopies; on the right side (Figure 2b) instead the corresponding whole soil images are shown at 1m resolution after applying deconvolution by kriging. Whereas for the raw B11 and B12 bands the soil maps were not retrieved due to pixel coarseness (spatial resolution 20 m).

Figure 3 instead represents the corresponding whole soil images of B11-B12 Sentinel2 bands at 1m resolution after applying deconvolution by kriging. In general, it can be said that the deconvoluted images show a high degree of coherence with each other, and no clear biases are evident, e.g. due to tree proximity, or artefacts resulting from the iterative fitting process and/or interpolation. Instead, all deconvoluted images display a sufficient degree of continuity derived from the raw images. Clearly the maps show different patterns of spatial association and distinct levels of reflectance depending on the wavelength.

Moreover, all deconvoluted maps exhibit spatial association structures, albeit with a different degree of smoothing depending on the pixel size of the raw data. Planet maps are characterized by great variability; however, it is possible to recognize structures characterized by higher reflectance values in the north-east and central area of the west side for bands 1, 2 and 3 (Figure 2a). Band 4 (Figure 2a) in the NIR range is both characterized by higher reflectance values than the previous maps and shows a different structure. The highest values are located to the NW and SW, and there is also a wide median strip parallel to the longest side of the field with lower values.

The B11 and B12 bands (Figure 3) in the SWIR range clearly differ from the previous bands both in their higher absolute values, because at these wavelengths the reflectance reaches the highest values in the entire spectrum, and in their spatial structure as well as in the higher degree of smoothing. This last, as mentioned above, depends on the greater difference between the size of the starting pixel (20 m) and that of the deconvoluted map (1 m). The spatial structure for both bands, which show a high degree of spatial correlation, is characterized by a broad continuous median strip of high values approximately corresponding to that strip of lower values observed for the NIR map of Planet.

Variations in soil reflectance can be due to 3 main factors: moisture, organic carbon content and clay. An increase in these factors should lead to a decrease in reflectance [87,88] . However, the response of soil is not as simple and linear as previously described, due to the multiple, multifaceted and non-linear interactions acting simultaneously. One effect here can be seen in the NIR map. If the large central area of high reflectivity in the SWIR can be interpreted as being due to a higher content of coarse textural component, as verified with the use of a Georadar [86], it is not clear how in the NIR map it corresponds to the lowest values of reflectance.

Before making predictions about the behavior of soils, one should first analyze its full spectrum to deduce its chemical composition. This type of information should then be supplemented by an investigation, e.g. geophysical surveys concerning the physical-structural characteristics of the soil.

3.2.2. Convolution of UAV Data

To estimate a LMC for the densely sampled UAV data (0.07-m resolution), it was necessary to reduce the number of samples since the calculation of the variogram is a computer time-consuming operation. For this purpose, a regular 1-m mesh grid was superimposed on the filtered soil data and only one point was chosen randomly within each cell [49]. This reduced the number of UAV pixels to 15.142, which enabled the calculation and modelling of direct and cross variograms. A point isotropic LMC was fitted to the data thus selected, comprising the following structures: a nugget effect, a spherical model with 6 m-range and a spherical model with-100 m range (Figure 4).

Like the Planet data, the UAV data detected the presence of two spatial structures: a short-range structure, linked to micro-variabilities such as mechanical variations in the soil and the presence of weeds, and another at longer range. The major difference among the three types of RS data consists in the longer-range structure: while Planet and B11 and B12 detect the presence of a structure at a range of about 50 m, probably to be assigned to the setting of the olive grove, UAV data disclose the existence of structures probably extending beyond the study field and linked to the soil characteristics of the surrounding area. This should make us reflect on the importance of using multi-spectral and different pixel satellite data to have a more complete characterization of the soil structure.

After the discretization of each 1m x 1m block into a 10 x 10 node grid, and the regularization of the variograms [70], the application of block cokriging produced the soil images of the UAV multispectral data shown in Figure 5.

The maps in the visible green and red (Figure 5) appear very similar to the corresponding ones in Planet with the macrostructures of high values recognizable to the NE and NW. The red edge map (Figure 5) is not present in Planet and is very similar to the NIR map (Figure 5). The latter, while evidencing the macrostructure of high values at the E, shows structural differences from its counterpart in Planet. It appears somewhat less structured, lacking the clear median structure of low values present in Planet, and with a greater component of random variability also found in the other UAV bands. This is probably due to the higher noise level of the UAV signal, which should be appropriately error-filtered.

3.2.3. Coregionalization Data Set

After migration of the raster data to the sample locations, a multivariate data set was obtained, the basic statistics of which are reported in Table 4.

From the analysis of Table 4, the soil is silty but with a good sandy component, which gives it good hydraulic and structural characteristics [86]. Furthermore, notice how the average hyperspectral estimates of the textural components are very close to the measured averages but show a considerable reduction in terms of variance.

In general, all variables are sufficiently centered and with no appreciable deviations from the normal distribution apart from the red bands in UAV and Planet, which are however little represented in the reflectance spectrum of this soil. It was therefore decided to apply the multivariate geostatistical analysis directly on the raw values of the variables without subjecting them to gaussian transformation, in order not to add another error component.

The linear correlation of the two variables of interest with the auxiliary variables (Table 5) is obviously higher with their hyperspectral estimates while it is generally quite low with the reflectance data from UAV, Planet and Sentinel2 data.

For clay, the highest correlations occur with red edge and NIR from UAV, and with B11. For sand, instead, the highest correlations are with red and green from UAV, while the relationship with B11 and B12 is not significant. This leads one to think that different soil texture components have different spectral behavior and therefore the optimal selection of the wavelengths used for survey cannot be separated from prior knowledge of the nature of the soils.

3.2.4. LMC of the Full Data Set

An isotropic LMC was fitted to the direct and cross experimental variograms of the multivariate data set, since the scarcity of samples prevented the detection of any anisotropies, consisting of three basic spatial structures: a nugget effect, a spherical model with a range of 9.91 m and a spherical model with a range of 39.55 m (Figure 6). As can be seen, there is sufficient consistency with the spatial structures detected by the individual types of radiation (UAV, Planet, Sentinel2) except for the longer-range component (100 m) detected only by the UAV data. However, of the total spatial variance only 14.47% and 16.53% are accounted for by the structured short- and long-range component, respectively, while the majority (69%) is nugget effect or unstructured random component. This is due to the paucity of soil samples (61) with an average sampling distance of 20 m, which for the variability of the study site is too coarse.

As far as the direct variograms (Figure 6) are concerned, a high dispersion of the experimental points with respect to the model is noted for practically all RS variables, due to the already observed scarcity of soil samples. This dispersion is reduced for the direct variograms of the soil variables and their hyperspectral estimates; however, the nugget effect is very high for the former and somewhat reduced for the latter.

The examination of the cross variograms related to clay and sand with the other variables (Figure 7) allows us to study the spatial relationships of the interest variables with the auxiliary ones.

The distance between the model (continuous line) and the dashed line of the intrinsic correlation (maximum spatial correlation between the two variables) [75] is taken as a quantitative measure of the degree of spatial correlation.

While a certain expected degree of spatial correlation can be seen between the soil variables and their hyperspectral estimates (40%-50%), unfortunately this correlation drops dramatically with the other variables. For clay, a somewhat decreasing relationship is observed with B11 and B12, which is positive up to distances of about 20 m, beyond which it becomes negative and then stabilizes on a negative value at about 40 m. Regarding sand, the spatial correlations with the RS variables are not significant even with B11 and B12. Maybe the direct correlations with red from UAVs and with NIR from Planet deserve some significance.

These results show the complexity of the relationship between soil and electromagnetic radiation, which besides depending on wavelength seems to be a function of scale. The latter can probably be interpreted as being due to the heterogeneity of the soil, i.e. the lack of stationarity associated with the presence of soil macrostructures of a different nature.

3.2.5. Maps of Clay and Sand from Data Fusion

Figure 8 shows the clay and sand maps obtained from the fusion (a) of all the variables under study. As can be seen, they are characterized by considerable short-range variability, despite the sparsity of the soil samples, transmitted by the auxiliary radiometric variables. However, macrostructures characterized by tendentially higher values in the north-east and in the median range for clay are noticeable. This highlights a certain positive relationship with B11 and B12 although the latter show a much higher degree of smoothing due to the low initial resolution.

In contrast, the sand (Figure 8) shows a certain similarity with the NIR band of Planet and the red band of UAV, with the median band showing lower values.

Figure 9 shows the maps of kriging standard deviations of the clay and sand cokriging estimates, respectively, from which the fiducial intervals [89] can be reconstructed, having initially supposed that the data distributions were sufficiently centered and symmetrical. These assumptions are therefore also reflected in the geometric distribution of the spatial association structures of error, which essentially reproduce the sampling pattern although the influence of auxiliary variables.

The maps above together with the assumptions mentioned allow for the determination of the fiducial ranges of the estimate and thus of the uncertainty associated with the prediction at each point. This evaluation is of fundamental importance, especially in risk analysis or decision-making [90].

3.2.6. Impact of Auxiliary Variables on Prediction of Soil Properties

To assess the benefit of the use of the auxiliary and RS variables on the prediction of the variables of interest, the previous estimates of clay and sand were compared with those obtained with the few samples of the variable alone (61) through cross-validation. Due to the sample scarcity, the direct variograms of clay and sand are pure nugget (supplementary material), therefore they are theoretically precluded from fitting an increasing monotonic mathematical model.

However, in geostatistics, variation at scales smaller than the sampling scale can be described by a variogram with nugget null and range equal to or less than the average sampling distance (20 m). Spherical models with a range of approximately 17 m were adapted in both cases (Figure 10).

Figure 8b shows the kriging maps of clay and sand. While recognizing a certain similarity with the homologous maps obtained from the data fusion (Figure 8a), a pronounced degree of smoothing and thus a noticeable reduction in spatial variance is clear.

The maps of kriging standard deviation maps (Figure 11) are very similar to the corresponding maps obtained from the fusion of the data.

The geometric distribution of the association structures of error did not change as was to be expected, while differences are only noticeable in the absolute values on the local scale.

To make the comparison between the two approaches more quantitative, the most relevant cross-validation statistics are reported in Table 6.

Examination of the table illustrates the marked improvement in precision and accuracy of estimation obtained from the fusion of both hyperspectral and multispectral auxiliary variables characterized by different wavelengths and different spatial resolutions. It is worth noting that statistically univariate kriging estimates could be globally accepted as unbiased and with an estimation error of an order of magnitude less than the standard deviation of the variable of interest. However, due to the sparse sampling and consequently the pure nugget character of the variogram, there is no association between measurement and prediction. This relationship is instead restored by the intervention of the auxiliary variables, albeit each with a different impact.

4. Discussion

In the present work, the problem of mixed pixels was addressed, i.e. when the RS images used had too coarse a resolution, while applications required a finer resolution, for example in Precision Agriculture. In the specific case study, it was necessary to select soil data from the remotely sensed images of an olive grove.

An important result that emerged from this work is that the best results are not always obtained using all available spectral information. In the case of hyperspectral data, it was in fact verified that the best estimation accuracy for clay was obtained not by using the entire spectrum but that specific ranges of wavelengths that also a previous work had shown to have the highest correlation with clay data. The same was true for sand even though the wavelength range changed, as was to be expected since the mineralogical composition of the textural constituents also changed.

An efficient way to build accurate spectroscopic models for predicting soil components may be to proceed in two steps: in the first step analyse the full spectrum by identifying the main absorption peaks; in the second step proceed to predict the variable of interest using different combinations of band ranges centred on the main absorption peaks.

Such a way of proceeding could be particularly advantageous in the choice of the filters to be applied to a hyperspectral radiometer on board a UAV to reduce its payload.

This selective choice is now possible using hyperspectral sensors, so the suggestion here is to take the objectives of the investigation into account. If these are very specific, e.g. searching for a certain soil or plant constituent, the suggestion is to use that narrow range of wavelengths that previous experimental results have shown to have a greater correlation with the variable of interest.

Finally, the work has emphasized, as was one of its main objectives, the importance of statistics in the joint analysis (fusion) of multi-source data, most of which are commonly remotely sensed data. The widespread availability of satellite data and the popularity of Machine Learning techniques have led to an increasing number of covariates being used, without any attention to their different measurement support and their different statistical properties. This can lead to incorrect conclusions, as spatial relationships have been shown to depend on the sampling support and sampling scale.

The use of auxiliary variables with a finer spatial resolution has made it possible to transfer soil variables, sampled at too coarse a scale, shorter-range information not otherwise obtainable even with interpolation techniques such as kriging. However, a legitimate question arises: is this added local variation related to the intrinsic variability of the variable of interest, or is it rather an artefact resulting from the use of auxiliary information? Indeed, this risk exists, but it will be the less the stronger the spatial correlation between the variable of interest and the auxiliary information. Unfortunately, in the case study, the correlation with the auxiliary radiometric variables was low. Therefore, although the results of the cross validation were encouraging, only an actual validation test, using an extended data set independent of the calibration data set, could prove the accuracy of the prediction.

The suggestion that can be made, as mentioned earlier, is not to use all the available auxiliary information, which today can also be considerable and cheap, but to make a prior selection based on statistical and physical considerations.

5. Conclusions

This paper addresses the problem of spatial statistical modelling applied to remote sensing data. An approach, based essentially on multivariate geostatistics, is presented that can fuse ground-based sample data and remote sensing grid data from multispectral satellite and UAV imagery with different spatial resolution. The problems of change of support and uncertainty, which are critical when jointly analyzing multi-source data with different statistical characteristics, are explicitly discussed. The geostatistical way of dealing with change of support and downscaling was to use a stochastic model of random function to represent spatial continua applied to remote sensing images.

This is just one way of addressing the crucial problem of spatial statistical modelling. It serves to illustrate the ongoing trend towards increasingly complex models that take advantage of both the computing power of computers and the development of statistical prediction models.

The contribution of this work was not only to present an efficient and flexible approach to the fusion of multi-source data, most of which is drawn from remote sensing. It was also intended to stimulate the development and application of spatial statistics to remote sensing.

Indeed, new statistical models of spatial data have recently been introduced and developed that would help not only in the understanding of reality but also in the interpretation of remote sensing data themselves. It is therefore desirable that remote sensing users be trained in these new techniques so that they know how to apply the most suitable spatial statistical model to extract the maximum information from their remote sensing data.