Satellite Multi-Sensor Data Fusion for Soil Clay Mapping Based on the Spectral Index and Spectral Bands Approaches

: Integrating satellite data at different resolutions (i.e., spatial, spectral, and temporal) can be a helpful technique for acquiring soil information from a synoptic point of view. This study aimed to evaluate the advantage of using satellite mono-and multi-sensor image fusion based on either spectral indices or entire spectra to predict the topsoil clay content. To this end, multispectral satellite images acquired by various sensors (i.e., Landsat-5 Thematic Mapper (TM), Landsat-8 Operational Land Imager (OLI), Advanced Spaceborne Thermal Emission and Reﬂection Radiometer (ASTER), and Sentinel2-MultiSpectral Instrument (S2-MSI)) have been used to assess their potential in identifying bare soil pixels over an area in northeastern Tunisia, the Lebna and Chiba catchments. A spectral index image and a spectral bands image are generated for each satellite sensor (i.e., TM, OLI, ASTER, and S2-MSI). Then, two multi-sensor satellite image fusions are generated, one from the spectral index images and the other from spectral bands. The resulting spectral index and spectral band images based on mono-and multi-sensor satellites are compared through their spectral patterns and ability to predict the topsoil clay content using the Multilayer Perceptron with backpropagation learning algorithm (MLP-BP) method. The results suggest that for clay content prediction: (i) the spectral bands’ images outperformed the spectral index images regardless of the used satellite sensor; (ii) the fused images derived from the spectral index or bands provided the best performances, with a 10% increase in the prediction accuracy; and (iii) the bare soil images obtained by the fusion of many multispectral sensor satellite images can be more beneﬁcial than using mono-sensor images. Soil maps elaborated via satellite multi-sensor data fusion might become a valuable tool for soil survey, land planning, management, and precision agriculture.


Introduction
Soil property prediction and mapping are necessary to determine a soil's production capabilities and, as a consequence, are excellent resources for many agricultural and environmental issues [1].In this context, digital soil mapping (DSM) [1,2] is widely used as a tool for soil information mapping around the world.DSM is the process of developing and populating spatial soil information systems using numerical models that infer the geographical and temporal changes of soil types and attributes from soil measurements, knowing, and associated environmental factors [1].The development of DSM as a valid option to meet the growing global need for geospatial soil data depends on its capacity to raise spatial resolutions, extend extents, and supply pertinent soil information [1].Remote sensing satellite images have been influential among the available input data for the DSM since they may give a chronological view of the scene under investigation [2,3], and enable assessing soil parameters over time [4].The extensive use of remote sensing satellite images may considerably increase the accuracy of DSM outputs.Indeed, remote sensing satellite images directly estimate some fundamental soil surface properties, such as soil texture (sand, clay, and silt content), carbon content, and soil salinity, e.g., [5][6][7][8].
In the literature, two common approaches have been widely applied to relate soil spectra extracted from remote sensing imagery to soil properties, including spectral indices, e.g., [9][10][11][12] and visible, near-infrared, and short wave infrared (VNIR/SWIR, 400-2500 nm) spectral bands (i.e., entire spectra) approach, e.g., [6][7][8]13].For instance, spectral datasets based on multi-and hyper-spectral images have been used by researchers to quantify the minerals and/or physicochemical properties of the soil.Spectral indices are based on the physical analysis of spectral reflectances, such as the slope or absorption band depth value, and by combining spectral reflectance from two or more wavelengths that indicate the relative abundance of features of interest, which allows the quantifying of minerals or physicochemical properties of the soils [14].For example, the absorption band depth values centered at 2206 and 2341 nm can be used to estimate the clay [15] and calcium carbonate (CaCO 3 ) [16] contents, respectively, whereas chromophores and iron oxides cause wide, weakly expressed absorption bands at wavelengths less than 1000 nm [17].
Many indices have been developed, such as the sand moisture index (SMI) [10] dedicated to sand content identification; the ferrous minerals index [18,19] and iron oxide index [18,19] created for iron content determination; the normalized soil moisture index (NSMI) [20] dedicated to moisture identification; the soil salinity spectral index (SSI) [21] used for salinity content identification; and the hue index [22] developed for soil color content identification.In addition, several spectral indices have been developed for clay content estimation, such as the Simple Ratio Clay Index (SRCI) [19], the SWIR fine particles index (FI) [10], and more recently, the Normalized Difference Clay Index (NDCI) [23].Reference [24] used linear regression models to establish the relationship between the NDCI derived from Landsat-TM images and clay content over the Kairouan plain in central Tunisia.In addition, [9] calculated the clay index from Hyperion bands at 2193, 2203, and 2214 nm.A SWIRFI index based on the spectral bands at 2130, 2205 and 2224 nm of the AISA-DUAL hyperspectral airborne data was employed by [14] in predicting clay content.
Remote sensing image fusion integrates spatial, spectral, and temporal data from airborne and satellite sensors with various resolutions to provide a multi-sensor image fusion that contains more extensive information than each of the original data [25].Data from multiple remote sensing sensors can be fused to significantly improve soil property prediction and mapping.For example, [4] indicated that using indices produced from multitemporal Landsat-8 data fusion instead of static terrain indices for DSM of soil parameters can improve the soil modeling and mapping for a wide range of applications.When combining gamma-ray spectrometry data with aerial photography photos, [26] discovered that the estimates of topsoil clay content improved.The authors of [27] established an empirical surface soil moisture (SSM) model via the fusion of soil moisture networks, remote sensing data (SMAP and Sentinel-1 images), and high-resolution land surface parameters (e.g., soil texture, terrain).
A possible fusion of multispectral images from many satellites can help achieve more detailed and abundant information than each satellite can provide, which can be helpful for improving soil studies and needs to be explored.This research aims (i) to assess the potential of multi-sensor bare soil images obtained by the fusion of various multispectral satellite images (i.e., Landsat-TM, Landsat-OLI, ASTER, and Sentinel-2A MSI) to predict clay content using spectral indices vs. spectral bands approaches and (ii) to compare their performances with those obtained from single-sensor images.Our hypothesis is that fusion satellite images provide high-quality data and more continuous, detailed, and abundant soil information than single-satellite images, based on both spectral indices and spectral band approaches, which can improve the digital mapping of soil properties.

Study Area
The study was conducted over the Lebna and Chiba catchments within the Cap-Bon region in northeastern Tunisia (Figure 1a,b).The Lebna and Chiba catchments cover an area of 210 and 220 km 2 , respectively [28].The catchments' upstream components, mainly in the northwestern region, have mountainous relief, while the downstream components, corresponding to the southeastern region, have a flat area (Figure 1b).The Jebel Abderrahmane is the western limit of the upstream section of the Lebna and Chiba watersheds, covered by scrubland and forest, with a maximum altitude of 650 m.The study area is an agricultural region with intense plant production, including fruit culture, field crops (cereals, legumes), fodder crops, grazing pastures, condiments, and olive cultivation.The area is distinguished by a broad range of parent materials, notably limestone, marls, sand, sandstones, dunes, and detrital deposits developed during the Cenozoic era.This lithological complexity has resulted in the evolution of several soil types, including Regosols, Calcic Cambisols, Vertisols, Cambisols, Rendzina, Ferralitic, and Iso-humic soils [29].Both catchments have a typical Mediterranean climate, with hot, dry summers and pleasant, chilly winters.Climate influence varies from a rainy regime at the top of the Jebel in the northwest, with more than 700 mm/year, to a semi-arid climate in the southeast, with less than 500 mm/year.The rainy season begins at the end of September and continues through February with moderate rain.The evaporative demand is roughly 1200 mm.In the summer, the highest temperature can reach 45 • C.

Soil Dataset and Sampling Collection
Soil samples were taken at a 0-5 cm depth during the soil studies conducted between October 2008 and November 2010.Five sub-samples were taken at random at every soil sample site inside a 10 × 10 m 2 square plot.The central geographical position was registered using a Garmin GPS navigation system and georeferenced.The clay content (granulometric fraction inferior to 2 µm) was measured using the pipette method (method NF 31-107, particle size distribution by sedimentation) after drying, homogenizing, and sieving to 2 mm of all the samples [30].

Multispectral Satellite Data and Bare Soil Images 2.3.1. Multispectral Satellite Data
Multispectral satellite images acquired by Landsat-5 TM [31], Landsat-8 OLI [31], ASTER [32], and Sentinel-2A MSI [33] were used for the purpose of this study (Table 1).The Landsat-TM and OLI images were level-1 precision terrain (L1TP) products.The ASTER images were level-1 precision terrain, corrected and registered, at-sensor radiance (AST_L1T) products.The Sentinel-2A MSI image was a level-1C top-of-atmosphere (TOA) reflectance data product.One scene of the Landsat-TM, Landsat-OLI, and Sentinel-2A data and two scenes of the ASTER data were selected to cover the entire study area and were registered on 14 August 2008, 27 July 2013, 30 August 2015, and 03 July 2004, respectively, (Table 1).The multispectral satellite images were chosen during the summer season to avert the soil moisture caused by rain, cloud cover, and the impact of green vegetation.All images were radiometrically calibrated, atmospherically corrected, and converted to surface reflectance using the QUick Atmospheric Correction (QUAC) method [34].The blue, green, red, and NIR1 bands of the Sentinel-2A image (i.e., bands 2, 3, 4, and 8A, respectively) were resampled to 20 m using the resample toolbox and the nearest neighbor technique [35] in ENVI software to guarantee that all Sentinel-2A spectral bands had an equal spatial resolution.The same treatment was carried out for the ASTER spectral bands with 15 and 30 m pixel sizes (band B1 to B3N and band B4 to B9), which were resampled to 30 m.The normalized difference vegetation index (NDVI) [36,37] was generated to mask the green vegetation in each image.Following a field visual inspection of different plots, a threshold value of NDVI equal to 0.22 was determined.Pixels with a reflectance value of less than 0.25 at 2100 nm, corresponding to bands 7 for the Landsat-TM and OLI, and bands 5 and 12 for the ASTER and Sentinel-2A, respectively, were designated dry vegetation coverage pixels and also masked.Pixels with a reflectance value of less than 0.28 at 1600 nm, corresponding to bands 5, 6, 4, and 11 of the Landsat-TM, OLI, ASTER, and Sentinel-2A, respectively, were considered as pixels covered by water and were also masked.The urban areas were also masked following a visual inspection.For further analysis, the green and dry vegetation, water, and urban areas were removed based on these diverse masks to maintain bare soil pixels in the satellite images, which we will refer to as the "bare soil images".

Multi-Sensor Image Fusion
The Sentinel-2A MSI bare soil image has a spatial resolution of 20 m, and to group with the Landsat-5 TM, Landsat-8 OLI, and ASTER bare soil images, it was resampled to 30 m using the nearest neighbor method [35] in ENVI software.The Landsat-5 TM, Landsat-8 OLI, ASTER, and Sentinel-2A bare soil images at 30 m pixel sizes were named TM_S, OLI_S, AST_S, and S2A_S, respectively.A fifth image was obtained by fusing bare soil pixels onto only the common areas between these images (i.e., TM_S, OLI_S, AST_S, and S2A_S images) and was denoted as Fusion_SB.The resulting Fusion_SB image with 31 spectral bands at a 30 m spatial resolution was performed using the ENVI layer stacking tool.

Spectral Index Images
Spectral index images of the clay contents, named TM_SI, OLI_SI, and S2A_SI, were generated using the spectral index proposed by [19] from the TM_S, OLI_S, and S2A_S bare soil images, respectively, as follows: where SWIR1 corresponded to the spectral range from 1550 to 1750 nm, 1570 to 1650 nm, and 1565 to 1655 nm, for the TM_S, OLI_S, and S2A_S images, respectively (Table 1).SWIR2 corresponded to the spectral range from 2080 to 2350 nm, 2110 to 2290 nm, and 2100 to 2280 nm, for the TM_S, OLI_S, and S2A_S images, respectively (Table 1).
A spectral index image of the clay contents, named AST_SI, was generated from the AST_S bare soil image, using the spectral index proposed by [38] as follows: where SWIR2, SWIR3, and SWIR4, corresponded to the spectral ranges from 2145 to 2185 nm, 2185 to 2225 nm, and 2235 to 2285 nm, respectively, and corresponding to the AST-band B5, AST-band B6, and AST-band B7, respectively.

Multi-Sensor Satellite Data Fusion of Spectral Index Images
All satellite clay spectral index images (i.e., TM_SI, OLI_SI, AST_SI, and S2A_SI) were fused on the only common area to form a single image with 30 m of spatial resolution, which reunited the four clay spectral index images, named Fusion_SI.This data fusion integrated the clay spectral index images derived from several multispectral remote sensing datasets at different spatial, spectral, and temporal resolutions to collect more abundant information.The spectral index processing and data fusion were implemented using the ENVI software.

Statistical Analysis 2.7.1. Soil Line Analysis
The soil line analysis approach [39] was used to determine if the applied masks may successfully recover bare soil pixels from each multispectral image (i.e., TM_S, OLI_S, AST_S, and S2A_S images).The reflectance values of the NIR and Red bands corresponding to the TM-bands B4 and B3, OLI-bands B5 and B4, AST-bands B3 and B2, and S2A-bands B8 and B4 (Table 1), were plotted on an XY graph, and regression analyses were conducted between them.The better the coefficient of determination (R 2 ) between the reflectance values of the NIR and Red bands, the nearer the scatter points were to the 1:1 line.Therefore, the applied masks were able to isolate the bare soil patches correctly.

Canonical Correlation Analysis
Canonical correlation analysis (CCA) [40,41] was carried out between the four multispectral images and only in the common bare soil pixels (i.e., on TM_S, OLI_S, AST_S, and S2A_S) to evaluate the spectral similarities between them.Each image's spectral bands were reduced to discrete canonical variables using a multivariate linear transformation, equivalent to substituting a group of variables with uncorrelated components [41].The higher the canonical correlations, the greater the spectral signatures were comparable amongst the multispectral images.This analysis was performed using SPSS software.

Pearson's Correlation Analysis
Pearson's correlation coefficient (r) was performed between the clay values and both the spectral band and spectral index values to evaluate the relevance of the derived remote sensing data for the clay content prediction.The p-value indicates whether the clay values were statistically significantly correlated with the spectral bands and spectral index in the Pearson's correlation analysis.

Spatial Regression Models
The predictions of the soil clay content using TM, OLI, AST, S2A, and fusion images were performed following two approaches: (a) using the spectral index of each multispectral image (i.e., TM_SI, OLI_SI, AST_SI, S2A_SI, and Fusion_SI) and (b) the VNIR-SWIR spectral bands of each multispectral image (i.e., TM_S, OLI_S, AST_S, S2A_S, and Fusion_SB), as predictors, for a total of ten predictor images (Figure 2).The soil sample datasets for each image were split into training (75%) and validation (25%) using the Unscrambler software package [42].The response variables (i.e., Y clay , clay content) were ordered in increasing rank.The lower clay content sample was first attributed to a validation set, followed by the next three samples in the training set.The process was then carried out by successively placing the next sample in the validation set and the subsequent three samples in the training set.This division guaranteed that the clay content in the training and validation sets was distributed similarly.In this study, we relied on AutoML to select the best regression model and corresponding hyper-parameters on the given data.We chose to run the Auto-WEKA [43], which is an AutoML system based on the popular machine learning and data mining platform, WEKA software [44].Auto-WEKA was the first AutoML system that considered the problem of simultaneously selecting a machine learning algorithm and optimizing its hyperparameters.
Among the popular machine learning algorithms, the Multi-Layer Perceptron (MLP) [45] network with a backpropagation (BP) algorithm has been widely used in practice for classification and regression problems.The MLP is a neural network that has a number of hidden layers of neurons.The connections between the neurons have a weight, which is obtained from a backpropagation algorithm.In the experiments, only one hidden layer was used in the MLP.The number of neurons in this layer, the learning rate, and momentum parameters were optimized through internal cross-validation in order to maximize the precision and accuracy of the results.The WEKA Multisearch package was used for this purpose.Feature selection was used to select the most pertinent and discriminative bands from the Fusion_SB image.Thus, it could identify and reduce as much irrelevant and redundant information as possible, thereby leading to the good predictive performance of the soil clay content.A mean decrease accuracy (MDA) [46] based on multilayer perceptron with the backpropagation learning algorithm (MLP-BP) evaluated the available band's characteristics and ranked all bands based on the variable importance criterion for each band.MDA utilized the capability of the MLP-BP to evaluate the spectral bands of the Fusion_SB image during the training stage.After the MDA rankings, the spectral bands that had zero or negative contributions to the predicting clay content were removed from the Fusion_SB image.The remaining bands were then examined in descending order of MLP-BP models, with the smallest and most accurate model being kept.In more detail, if the introduction of a spectral band provided no increase in the model accuracy, then the said band was removed from the Fusion_SB image.At the final step of the procedure, we obtained a fusion image with the most discriminant features and this was denoted as Fusion_S.

Model Performances Assessment
The performance in the validation sets of the MLP-BP models was evaluated by the coefficient of determination of validation (R 2 val ), the root mean square error of prediction (RMSEP), the mean absolute error of prediction (MAEP), the bias of prediction on validation sets, the ratio of the performance to the deviation on the validation set (RPD, [47]) and the ratio of performance to interquartile on the validation set (RPIQ, [48]).The predicted clay values and the produced maps were realized using the prediction tools of the WEKA software [44] and the ArcGis software package [49].

Spatial Structure Analysis
The spatial structure of the predicted clay contents was evaluated by examining their variogram characteristics, which quantified their spatial dependence.The semi-variance for the soil clay values to be displayed on a variogram was calculated using the formula [50]: where h is the lag distance, γ(h) is the semi-variance of the soil clay content at a lag distance h, N(h) is the total population of sample pairs at lag distance h, z is the value of the clay content, x is the position of point i on the landscape, and z(x i + h) is the clay value at the position x i + h.The variograms were calculated for the observed values dataset and maps estimated by spectral index (i.e., TM_SI, OLI_SI, AST_SI, and S2A_SI), spectral bands (i.e., TM_S, OLI_S, AST_S, and S2A_S), and image data fusion (i.e., Fusion_SI and Fusion_S) using Surfer-15 software (Golden Software, LLC).

Bare Soil Images and Spectral Pattern Descriptions
The TM image had the highest bare soil pixels percentage (63.14%),followed by the S2A (57.37%),AST (52.50%), and OLI (41.62%) images (Figure 3).The bare soil's surface variation was related to the different acquisition years of the multispectral scenes.According to the soil line analysis, the OLI image best described the bare soil patterns in the study region, with the highest coefficient of determination (R 2 = 0.90).In contrast, the AST image had the lowest representation (R 2 = 0.79) of all the examined images (Figure 3).The four multispectral remote sensing datasets had 34.6% of the common bare soil pixels.The canonical correlation analysis between these common bare soil pixels showed that the TM and OLI, as well as the OLI and S2A, had the highest coefficient (0.90) (Table 2).The percentage of bare soil pixels over both fusion images (i.e., Fusion_SI and Fusion_SB) was equal to 34.6%, which corresponds to the common bare soil area between all derivedmultispectral images.Based on the spectral characteristics of multispectral satellite images (i.e., TM, OLI, AST, and S2A images) (Figure 4), the reflectance values rose as the clay content dropped [8].The TM and OLI spectra appeared very similar (Figure 4a,b).The reflectance values of the sandy soils were greater than those of the clay soils; however, the differences were noticeable in the NIR-SWIR wavebands regardless of the sensor (Figure 4).

Descriptive Statistics of Clay Soil Samples and Correlation with Bare Soil Images
The entire soil sample dataset contained 262 soil measurements.Due to the different number of bare soil pixels in each image, the dataset was reduced to 221, 174, 196, 218, and 149 soil data for the TM, OLI, AST, S2A, and fusion images, respectively (Table 3).The descriptive statistics results of the sample datasets used for each satellite image revealed a significant degree of similarity between the clay content values of the soil sample datasets.Based on the analysis of variance (ANOVA) results, no significant variations in the analyzed clay content values were observed at the 0.01 alpha level.The coefficient of variation (CV, reported in percent) was between 41 and 45% in all datasets.The TM soil samples had the highest CV, indicating that the TM soil sample set was the most varied in these datasets.In all datasets, the interquartile and standard deviation ranged from 283 to 302 g/kg and 179 to 182 g/kg, respectively.Finally, a soil sample with a clay content of 772 g/kg was reported in all datasets (Table 3).The Pearson's correlation coefficients between the spectral bands of the TM, OLI, AST, and S2A bare soil images and the measured clay values (Table 4) showed the significantly correlated bands of multispectral images with the measured clay values (p-value less than 0.05).All the spectral bands of the satellites' multispectral images revealed a negative correlation with the soil clay content values.All datasets' clay values demonstrated high correlation coefficients with the Red, RedEdge, NIR, and SWIR bands of the multispectral images (Table 4).The clay content was most correlated with the spectral bands TM-B7, OLI-B7, AST-B8, and S2A-B12, with r values equal to −0.72, −0.71, −0.67, and −0.72, respectively (Table 4).This most correlation is due to the clay minerals' characteristic hydroxyl absorption around 2200 nm, caused by a mixture of vibrations related to the OH and OH-Al-OH bonds [15,51].All clay value datasets were significantly correlated with the clay spectral index of each multispectral image, although the S2A_SI had the highest correlation (r = −0.71,Table 5).

Selected Features from the Spectral Bands of the Fusion Image
The importance of each spectral band of the Fusion_SB image was derived from the MDA based on the MLP-BP algorithm and plotted on Figure 5 in order to investigate the feature rankings.The spectral bands S2A-band B12 (2100-2280 nm) and TM-band B7 (2080-2350 nm) were the most important for predicting the soil clay content.This spectral range had a central wavelength of around 2200 nm, which corresponds to an absorption characteristic caused by the combination of OH-Al bending and OH stretching modes [15].Other variables that exerted important impacts were the S2A-band B11, followed by the OLI-band B7 and AST-band B8 (Figure 5).The AST-band B7, AST-band B6, TM-band B5, AST-band B5, S2A-band B8, S2A-band B8A, AST-band B4, and OLI-band B6 were identified as the factors group with a moderate importance to predict soil clay content.The other variables (i.e., S2A-band B7, AST-band B9, S2A-band B6, S2A-band B5, and S2A-band B4) emerged as having weak impacts.Although these variables might not be particularly predictive by themselves alone, some can increase the model accuracy through their interactions with other variables [52].Additionally, thirteen bands had no impact on the prediction performance (Figure 5).Finally, from the Fusion_SB image, the variable selection using the MDA based on MLP-BP selected a feature subset of eight attributes for the clay prediction.This feature subset included the spectral bands S2A-band B12, TM-band B7, OLI-band B7, TM-band B5, S2A-band B8, OLI-band B6, S2A-band B7, and AST-band B9, which were selected to generate the Fusion_S image.

Spectral Index vs. Spectral Bands Prediction Performances
Based on the Auto-WEKA results, the MLP-BP was selected as the best regression model for the datasets.The MLP-BP models generated from spectral index images (i.e., TM_SI, OLI_SI, AST_SI, and S2_SI) performed poorly regardless of the sensor (0.43 ≤ R 2 val ≤ 0.51, and 126.42 ≤ RMSEP ≤ 140.24 g/kg, Table 6); however, the MLP-BP model built from the spectral index image fusion (Fusion_SI) provided good performance (R 2 val = 0.61 and RMSEP = 114.18g/kg, Table 6), outperforming the MLP-BP models built from individual spectral index images (Table 6).
The MLP-BP models created from the multispectral images (TM_S, OLI_S, AST_S, and S2A_S) showed accurate results (0.60 ≤ R 2 val ≤ 0.71, and 96.45 ≤ RMSEP ≤ 115.64 g/kg, Table 6), whereas the clay prediction model built using the S2A_S spectra provided the higher performance, taking advantage of the sentinel-2 MSI sensor's ten spectral bands.The clay prediction model created using the AST_S spectra performed worse, although the ASTER sensor had nine spectral bands.The MLP-BP model built from the fusion image (Fusion_S) performed excellently (R 2 val = 0.78, and RMSEP = 87.84g/kg, Table 6), outperforming the MLP-BP models built from individual multispectral images (Table 6).

Predicted Soil Maps and Spatial Structure Analysis
The MLP-BP models were applied to their corresponding bare soil images to provide the predicted clay content maps (Figure 6).The predicted clay content maps derived from the spectral index approach showed a globally similar distribution of clay spatial patterns regardless of the percentage of bare soil areas, as was the case for the maps derived from the VNIR-SWIR spectral bands approach.On these predicted clay maps, the south-eastern part of the studied area presented low concentrations of soil clay, whereas the north-eastern part of the area presented high concentrations, regardless of the map (Figure 6).Moreover, the clay content map predicted using the Fusion_S showed clay content values slightly higher than the other maps, regardless of the maps region.
For each satellite sensor, both the spectral index and VNIR-SWIR spectral bands images were able to represent the variability of the soil clay content; however, with the spectral images, the interquartile and standard deviation values were higher than those obtained with the spectral index images.Both the Fusion_SI and Fusion_S clay content maps, which had the best performance (Table 6), showed high interquartile and standard deviation values in contrast to the other predicted maps.Depending on the images used to create the predicted clay maps, the IQ varied from 187 g/kg (obtained from the AST_SI image) to 279 g/kg (obtained from the Fusion_S image).Thus, the standard deviation varied from 120 g/kg (obtained from the AST_SI image) to 161 g/kg (obtained from the Fusion_S image) (Figure 6).
Figure 7 shows a zoom over a small area of the predicted clay maps, characterized by strong variations in the soil patterns on a small scale, with a close succession of clay-rich areas (brown color) and clay-poor areas (sky blue color) oriented northwest/southeast, corresponding to marl and sandstone outcrops, respectively.Following the VNIR-SWIR spectral bands approach, the predicted clay maps exhibited more clearly short-range variations in soil patterns, with a close succession of clay-rich areas and clay-poor areas, than the maps following the spectral index approach (Figure 7).Between the clay maps based on the fusion of the spectral index and VNIR-SWIR spectral bands images, the Fusion_S showed slightly higher clay values than the Fusion_SI values.It should be noted that the changes in the bare soil spatial coverage are small among the Landsat-5 TM, Landsat-8 OLI, ASTER, and Sentinel-2A images (Figure 7).Histograms analysis of the predicted maps over the Lebna and Chiba catchments showed a normal distribution of clay content, with the exception of the histogram of predicted from the TM_SI around 500 g/kg, which represents an extremely weak density (Figure 8).The histograms of predicted values from the TM, AST, and S2A images showed a histogram peak of around 300 g/kg, while the histograms of predicted values from the OLI image showed a histogram peak of around 500 g/kg (Figure 8).This may be due to the training database's clay values, which were slightly higher than the other databases (see Table 3 and Figure 8).Moreover, the distribution of clay contents predicted using the spectral index images (TM_SI, OLI_SI, AST_SI, and S2A_SI) shifted to slightly higher values compared to those obtained from the VNIR-SWIR spectral bands images.This means that the pedological variability of the study area is more represented in the predicted clay content maps obtained using the VNIR-SWIR spectral bands images (TM_S, OLI_S, AST_S, and S2A_S) than in the spectral index images (TM_SI, OLI_SI, AST_SI, and S2A_SI).Thus, the maps predicted using the fusion images exhibited a wider distribution of clay content values than the other images (Figure 8).Finally, the variograms of the measured clay content values (dots in Figure 9) had higher variance than the predicted values, regardless of the sample datasets used for each satellite image.Thus, the spatial structures of the estimated clay contents obtained by the VNIR-SWIR spectral bands images exhibited very high sills compared to the spectral index images (Figure 9).By contrast, the variance obtained by the spectral index image fusion exhibited a very high sill and strong variations compared to those obtained by a single sensor spectral index image.Finally, the variogram obtained by the VNIR-SWIR spectral bands image fusion was higher.

Clay Predictions Depending on the Multispectral Sensor
The models based on multispectral images (i.e., TM_S, OLI_S, AST_S, and S2A_S) provided acceptable accuracies with R 2 val values between 0.60 and 0.71, RMSEP values between 96.45 and 115.64 g/kg, MAEP values between 73.29 and 91.86 g/kg, RPD values between 1.60 and 1.90, and RPIQ values between 2.59 and 3.01 (Table 6).This performance variance could be explained by differences in satellite remote sensing datasets' acquisition dates (i.e., date ranges), as previously shown by [53].Moreover, though the bands of each sensor cover similar spectral domains, they were not exactly the same (Table 1).Additionally, different acquisition conditions and view angles between the sensors can result in differences between the recorded datasets [54].The model performance obtained from the S2A_S image was superior to that obtained by [6], who also used the ten bands of the Sentinel-2A (S2A) and a PLSR model.Similar results were obtained by [8] using Sentinel-2 MSI and Landsat-8 OLI images for clay content prediction, with R 2 values of 0.68 and 0.62, respectively.Moreover, [29] also found similar performance using the MLR model and nine bands of the ASTER image, while [55] obtained a high prediction accuracy of soil clay content using a mean spectral reflectance from bare soil pixels along a Landsat-TM time series.

Clay Predictions Obtained from Spectral Indexes vs. Entire Spectra
Whatever the sensor, the models built using the entire VNIR-SWIR spectra outperformed the models built using spectral indexes (Table 6).The models built from the four spectral index images (i.e., TM_SI, OLI_SI, AST_SI, and S2A_SI) performed poorly (0.43 ≤ R 2 val ≤ 0.51, 126.42 ≤ RMSEP ≤ 140.24 g/kg, 99.65 ≤ MAEP ≤ 112.29 g/kg, 1.34 ≤ RPD ≤ 1.45, and 2.14 ≤ RPIQ ≤ 2.30, respectively; Table 6), and this could be explained by the lack of spectral information.These prediction accuracy results were in agreement with those obtained by [23], who used a Simple Ratio Clay Index (SRCI) and a Normalized Difference Clay Index (NDCI) based on Landsat-OLI imagery for clay content mapping in Indonesia.Accuracy assessment of the regression models based on independent validation showed that both clay indexes provided low performances, with R 2 values of 0.42 and 0.40, respectively [23]; however, [24] used Simple Linear Regression (SLR) to establish the relationships between many NDCI-Landsat TM and clay content over the Kairouan plain in Tunisia and obtained results with R 2 values that varied between 0.48 and 0.73.

Added-Value of Fusion Images
The best-performer models were obtained using multi-sensor image fusion incorporating spectral bands or spectral index images (Table 6).Effectively, the model with the best performance for the spectral bands approach was founded on the Fusion_S image, with an R 2 val of 0.78, an RMSEP of 87.84 g/kg, an MAEP of 69.99 g/kg, an RPD of 2.12, and an RPIQ of 3.25 (Table 6).Similarly, for the spectral index approach, the model based on the Fusion_SI image was superior to those found on the TM_SI, OLI_SI, AST_SI, and S2A_SI (Table 6).
On the other hand, based on the same soil sample set (i.e., 149 samples), the models using multispectral images (i.e., TM_S, OLI_S, AST_S, and S2A_S) provided acceptable results with R 2 val values between 0.48 and 0.69, RMSEP values between 99.94 and 132.31 g/kg, MAEP values between 82.89 and 109.93 g/kg, RPD values between 1.40 and 1.86, and RPIQ values between 2.16 and 2.86 (results not shown).These performances, however, always remained lower than those obtained using the Fusion_S (Table 6).The models based on the four spectral index datasets (i.e., TM_SI, OLI_SI, AST_SI, and S2A_SI) and the same soil samples (i.e., 149 samples) performed poorly in estimating the clay content with R 2 val values between 0.21 and 0.53, RMSEP values between 125.03 and 162.89 g/kg, MAEP values between 99.36 and 130.59 g/kg, RPD values between 1.14 and 1.49, and RPIQ values between 1.75 and 2.28 (results not shown), and these performances always remained lower than those obtained using Fusion_SI (Table 6).
The results of this paper suggest that the models based on satellite multi-sensor data fusion were the most precise for prediction and mapping of soil clay content.Moreover, the improved accuracy of the results obtained using fusion images might be attributed to complementary spectral, spatial, and temporal information derived from the fusion of various satellite sensor datasets.For example, [4] demonstrated the potential of using fusion indices derived from multitemporal Landsat 8 imagery to predict static soil properties such as OC, sand, and CCE.Meanwhile [26] investigated multisource data fusion for topsoil clay mappings, such as proximally measured gamma (γ) ray spectrometry, apparent electrical conductivity (ECa), four terrain attributes, and the digital numbers (DNs) of an aerial photo.The results improved when the gamma-ray spectrometry data was integrated with the digital numbers (DNs) of the aerial photos [26].Ref. [27] established an empirical surface soil moisture (SSM) model via the fusion of soil moisture networks, remote sensing data, and high-resolution land surface parameters.The results indicated that the data fusion-based model could retrieve the SSM temporal dynamics [27].

Future Research
In this study, the best results were obtained via the fusion of four multispectral satellite images with feature selection using an MDA based on MLP-BP, which allowed the selecting of the most significant, pertinent, and discriminative spectral bands to predict clay content (i.e., eight bands), instead of just using six bands for the TM or OLI, for example.Thus, the modelling process produced an increase of about 10% in the prediction accuracy compared to that obtained from single-sensor images.Nonetheless, only the variable selection using the MDA based on the MLP-BP algorithm was used to select the key features impacting clay prediction; however, future research could consider using additional feature selection methods, including feature subset evaluation and feature ranking methods belonging to three groups: filter, wrapper, and embedded techniques.Finally, the incorporation of additional data sources to generate the fusion image, such as thermal infrared bands, synthetic aperture radar (SAR) and light detection and ranging (LiDAR) data, can improve the accuracy of the results for soil property prediction.

Conclusions
Spectral bands images restricted to bare soil pixels, derived from Landsat-TM, Landsat-OLI, ASTER, and the Sentinel-2A satellites, were more suitable for obtaining the soil clay map than those derived from the spectral index images of the same satellites.Moreover, the satellite's multi-sensor data fusion, which included multispectral bands or spectral index images, was revealed to be a promising tool for mapping soil attributes.Multi-sensor images generated by the fusion of multispectral satellite data can provide a voluminous amount of bare soil information for DSM.To reinforce these conclusions and the guidelines, this study could be enlarged to (i) evaluate additional feature selection algorithms and methods, (ii) integrate additional data sources into the fusion image, (iii) assess the influence of the spatial resolution of the fusion image on soil property prediction, and (iv) include additional soil properties and additional pedological settings.

Figure 1 .
Figure 1.Location of the study area over (a) Tunisia (red areas) and (b) Sentinel-2A MSI (Red, Green, Blue: bands 4, 3, 2, respectively) image acquired in 2015 over the Lebna (northeastern) and Chiba (southeastern) bordering watersheds (black polygons).The yellow dots represent the locations of the soil samples.

Figure 2 .
Figure 2. Flowchart methodology for the study.

Figure 3 .
Figure 3. Bare soil pixels of the TM image (band B1: 485 nm), OLI image (band B2: 482.6 nm), ASTER image (band B1: 556 nm) and S2A image (band B2: 496.6 nm) at 30 m spatial resolution.The percentage values correspond to areas with bare soil.Relationship between NIR and Red bands are presented to show the soil line for each image.The coefficient of determination (R 2 ) value is as an indication of bare soil representation.

Figure 4 .
Figure 4. Comparison of clay and sandy soils spectral patterns (blue-cyan and orange-yellow spectra, respectively) based on the same pixels of (a) Landsat-5 (TM), (b) Landsat-8 (OLI), (c) ASTER (AST), and (d) Sentinel-2A (S2A) images.Note that the TM and OLI images have six bands, while the AST and S2A images have nine and ten bands, respectively.

Figure 5 .
Figure 5. Variable importance using mean decrease accuracy (MDA) based on Multilayer Perceptron (MLP) with the backpropagation (BP) learning algorithm.

Figure 8 .
Figure 8. Histograms of the clay content (a) estimated by Landsat-5 (TM) and 165 measured sites of the training set (165M), (b) Landsat-8 (OLI) and 130 measured sites of the training set (130M), (c) ASTER (AST) and 147 measured sites of the training set (147M), (d) Sentinel-2A (S2A) and 163 measured sites of the training set (163M), and (e) fusion images (Fusion) and 111 measured sites of the training set (111M) and following the spectral index (SI) and VNIR-SWIR spectral bands (S) approaches.

Table 2 .
Canonical correlation analysis between the common bare soil area of multispectral images, with the first canonical function.

Table 3 .
Statistics of the sample datasets used for the satellite images.

Table 4 .
Pearson's correlation coefficient (r) between the spectral bands of TM, OLI, AST, and S2A bare soil images and measured clay value datasets.

Table 5 .
Pearson's correlation coefficient (r) between the spectral index of TM, OLI, AST, and S2A bare soil images and measured clay value datasets.

Table 6 .
Model performances of clay content (in g/kg) for validation datasets following spectral index and spectral bands approaches.Notes: SI: spectral index.S: spectral bands.NB: number of spectral indices or bands used in predictions.NS: number of samples.H: number of neurons in the hidden layer.L: learning rate.M: momentum.The highest performances are in bold.