1.

INTRODUCTION

Grasslands should be managed in a manner to promote simultaneously both soil and rangeland quality, especially in countries that base their agricultural production on livestock farming. Therefore, there is a need for improved management of grazing resources which can result in productive bio-diverse grasslands, also mitigating soil erosion and enhancing carbon sequestration [1]. In this context, recent studies highlight the need for efficient monitoring of key soil descriptors to efficiently estimate soil loss by water erosion based on the Revised Universal Soil Loss Equation (RUSLE) model with higher accuracy [2] while state-of-the-art Artificial Intelligence (AI) and data mining techniques are considered critical to achieve this goal [3]. More specific, the soil layers then introduced into the RUSLE’s soil erodibility factor (K-factor), producing a more reliable soil erosion layer and with improved spatial resolution.

The existing soil explicit indicators provided information in moderate performance and coarse resolution, mainly relying on environmental covariates fed into machine learning models [4]. Several techniques, such as Random Forest and eXtreme Gradient Boosting (XGBoost), have resulted promising results, while deep learning techniques have also been used with some of them proposing a synergistic framework. On the other hand, the Sentinel-2 satellite has been extensively employed to map the soil texture from multispectral data. However, simple merging techniques may not fully exploit the complementary nature of the data, potentially resulting to information loss or misinterpretation.

Therefore, novel approaches are needed to tackle the challenge posed by the synergistic framework of Earth Observation (EO) data analytics, which require effective fusion of multispectral data with environmental and topographical covariates. In this study, we explore the potential of employing a deep learning architecture to obtain a new data representation from spaceborne Sentinel-2 information for the regression task. Concurrently, we employ a XGBoost regressor, using features extracted by a convolutional neural network (CNN).

2.

MATERIALS AND METHODS

2.3

Methodological Approach

2.3.1

AI Approach

The proposed approach comprises two steps starting from data collection and then the regression analysis using a hybrid CNN-XGBoost algorithm (Figure 1).

Figure 1.

Overall approach from data collection to regression

CNN, as feature generator, is initiated with an input layer designed to accept one-dimensional features with a length of 12 that corresponds to Sentinel-2 bands. Then, two layers with kernel sizes of 3x1 followed by Leaky ReLU activation functions have been utilized, having 48, and 24 filters, respectively. Subsequently, the feature maps are flattened and passed through two fully connected layers with 128, and 48 neurons, respectively, employing Leaky ReLU activation functions. Then, we utilized the spectral generated features along with the environmental covariates to feed an XGBoost algorithm which builds a series of decision trees sequentially. It should be mentioned that XGBoost is considered as a powerful ensemble model, since each subsequent tree adjusts the errors made by the previous one at each step. After a grid search the following hyper-parameters have been selected for the XGBoost: the ‘number of estimators’ was set to 60, the ‘minimum samples per split’ to 4, the ‘minimum samples per leaf’ to 2, the ‘maximum depth’ to 4, and the ‘learning rate’ to 0.05. The same hyperparameters were used to calibrate an XGBoost model that got as input features together the Sentinel-2 data and environmental covariates in order to make a comparison with our approach. We trained the model based on 50 random splits. The assessment of the regression performances is done considering the Root Mean Squared Error(RMSE), the concordance correlation coefficients (CCC) and the Ratio of Performance to Inter Quartile distance (RPIQ).

2.3.2

RUSLE Approach

The annual soil loss was estimated following the RUSLE empirical equation by using improved AI geospatial layers and open access EO datasets [9]:

where A is the average annual soil loss (ton/ha/yr) and following factors explained below:

• • R-factor—Rainfall erosivity

  For the R-factor (MJ·mm/ha/h/yr), the ERA-5 dataset was used covering a timeseries from 2007 to 2023 with 1km of spatial resolution. while the equation of Wichmeier and Smith (1978) [10] was used:

  

  where Pm is the monthly precipitation in mm and Pa is the annual precipitation in mm.

• • K-factor—Soil erodibility

  The K-factor [(t·ha·h)/(ha·MJ·mm)] was calculated by using the enhanced AI layers (SOC and soil texture) generated in this work, followed the methodology proposed by [2] and based on the equation described in [11]:

  

  where Sand, Silt, Clay and OC are the percentage contents of sand, silt, clay and organic carbon, while SN equates to 1 – Sand/100.

• • C-factor—Crop cover and management

  In this study, the C-factor is determined by taking the median value of the NDVI index, which is calculated using bands B4 and B8 from multi-temporal Sentinel-2 imagery with a spatial resolution of 10m, using the following epxpression:

  

  where α and b are unitless parameters that determine the shape of the curve relating to the values of the NDVI and C factors.

• • LS-factor—Slope length and steepness

  The RUSLE topographic factor (LS) defined by the combination of the slope length factor (L) and slope steepness factor (S), calculated based on the Copernicus DEM (30m) while the flow accumulation, downloaded from https://github.com/davidbrochart/flow_acc_3s, which is a 3 s flow accumulation derived form HydroSHEDS. At the end, the equation proposed by [12] was used:

  
  
  
  

  where L is the slope length factor; S is the slope steepness factor; λ is the horizontal projected slope length (based on flow accumulation); m is a variable length-slope exponent; β is a factor that varies with slope gradient; and θ is the slope angle.

• • P-factor—Support practices

  In this study, the estimated P-factor data developed by [13] was used which has a spatial resolution of 1km (https://esdac.jrc.ec.europa.eu/themes/support-practices-factor).

3.

RESUTLS

3.1

AI Results

Table 2 present the results obtained by XGBoost and CNN-XGBoost approaches on the Elassona region. We can notice that CNN-XGBoost approach outperforms the XGBoost for the prediction for SOC, Clay and Sand, with the only exception of Silt, where is is gained a lower RMSE. High values of SOC cannot accurately predicted, however our region is not characterized by SOC content values bigger than 2 g/kg.

Table 2.

Evaluation metrics considering CNN-XGBoost approach and XGBoost competing method

Parameters	XGBoost	CNN-XGBoost
	CCC	RMSE	RPIQ	CCC	RMSE	RPIQ
Clay	0.63	75.29	1.86	0.64	74.63	1.88
Sand	0.58	137.96	1.96	0.61	133.75	2.02
Silt	0.55	87.75	2.17	0.51	88.21	2.13
SOC	0.11	1.98	1.12	0.13	1.97	1.13

The best models were employed to produce spatial representations for soil texture and SOC content. Figure 2 illustrates the maps alongside the distribution of estimated values within the region of interest. The outcomes align with the distribution patterns observed in the Greek soil data archive for the specified soil variables.

Figure 2.

Soil spatial descriptors as estimated by the CNN-XGBoost algorithm

3.2

Soil loss estimations

First, each of the RUSLE’s factor was calculated (see sect. 2.3.2) and the multiplication of all the factors resulted in the final soil erosion map generation with 10 m of spatial resolution (Figure 3). Our study area characterized mainly by low to medium rainfall erosivity values while the LS-factor has mainly high values due to the steep slopes, that prevail in the area, and in combination with the intense stream network. The area has an average soil loss value of 4.6 ton/ha/yr with min and max values of 0 to 51 respectively. The 71% of the total area has a soil loss less than 5 ton/ha/yr while the 4% suffers of soil loss more than 20 ton/ha/yr.

Figure 3.

Spatial distribution of the produced RUSLE factors: a) C-factor, b) R-factor, c) LS-factor, d) K-factor and e) soil erosion

Considering that the spatial resolution will be the primary distinction between our soil erosion map and current available products, we opted to re-evaluate the readily available products. In that regard, we performed an additional simulation using SOC and soil texture layers from the SoilGrids platform [https://soilgrids.org/] and keeping the same datasets for the generation of the rest RUSLE factors (C, R, P and LS) producing a final soil erosion map with 250 m resolution. Although a similar pattern exists in the soil erosion products (Figure 4), critical differences are existed that can lead to erroneous estimations due to variations in the spatial distribution of soil layers and map resolution. For instance, certain areas within the sub-region categorized with low soil loss values <1 may exhibit significantly higher soil loss in coarser-resolution products.

Figure 4.

Soil erosion maps produced with different spatial resolution, using the current soil spatial products from SoilGrids platform (left) and the products as generated by the proposed CNN-XGBoost approach

4.

DISCUSSION

The accuracy of the estimations for the soil texture were deemed acceptable, with SOC content exhibiting the lowest predictive performance with an RMSE 1.98 g/kg (Table 2). Having our results in comparison with other studies in the literature, we can notice that the performance for soil texture is similar to where the research performed at regional scale. A significant percentage of recent studies have relied on data sourced directly from the specific region rather than from national archives, as we did. Therefore, our results can be attributed to this difference in data sources. Therefore, the absence of ground data should be noted as a limitation of the current study. Incorporating both field observations with EO data would enhance our ability to quantify and calibrate soil erosion AI-models more effectively.

Based on our results, Elassona region is generally characterized by low to moderate erosion levels (Figure 4). This is a significant result compared to the current estimates that result to significant uncertainties and higher soil loss estimations since they are long-term averages performed with empirical models. Through the lens of emphasizing soil ecosystem protection via dedicated monitoring as advocated by a set of policies, it is suggested that improved estimations, such as those proposed here, which integrate AI and high spatial resolution data, should be utilized. Enhanced soil products have demonstrated notable improvements when integrated into physical process models. This integration presents an avenue for further exploration, particularly in the realm of soil loss estimation [14]. The products illustrated in Figure 3 allow us to offer more timely and consistent estimations, facilitating the monitoring of soil loss on a scale able to propose best practices. Our model, which effectively fused spectral information with environmental covariates, has also facilitated the interpretation of results. Therefore, Shapley analysis [15] can be employed to further enhance our understanding of the contributions of different variables to the predictive outcomes. Moreover, further studies could explore techniques that integrate additional factors influencing soil erosion, such as management practices, through an interpretable approach [16] of post hoc analysis to derive recommendations on the most effective management practices (e.g., cover crops, buffer strips, etc.) for reducing soil loss [17].

5.

CONCLUSION

In our study, we proposed an approach to enhance the spatial representation of soil loss estimation by water erosion by leveraging cutting-edge deep learning techniques. Our approach involved integrating a CNN, to handle multispectral data from Sentinel-2, with a XGBoost regressor complemented by landscape features and bioclimatic variables. By training our CNN-XGBoost model on Greek soil samples from the LUCAS 2015 dataset, we achieved an improvement in soil input layers, resulting in a reduction of approximately 5% in RMSE. These enhanced spatial products were seamlessly integrated into the RUSLE framework, thereby enhancing the soil erodibility factor and yielding a soil erosion layer with unprecedented spatial resolution (10m). Our field mapping endeavors in Elassona, Greece, provided compelling evidence of the efficacy of our approach compared to existing soil products, that overestimate the current situation resulting also significant uncertainties. The current approach demonstrates a transformative potential in soil erosion monitoring able to evaluate the impact of various management practices and restoration policies.