Mapping of Rill Erosion of the Middle Volga (Russia) Region Using Deep Neural Network

Soil erosion worldwide is an intense, poorly controlled process. In many respects, this is a consequence of the lack of up-to-date high-resolution erosion maps. All over the world, the problem of insufficient information is solved in different ways, mainly on a point-by-point basis, within local areas. Extrapolation of the results obtained locally to a more extensive territory produces inevitable uncertainties and errors. For the anthropogenic-developed part of Russia, this problem is especially urgent because the assessment of the intensity of erosion processes, even with the use of erosion models, does not reach the necessary scale due to the lack of all the required global large-scale remote sensing data and the complexity of considering regional features of erosion processes over such vast areas. This study aims to propose a new methodology for large-scale automated mapping of rill erosion networks based on Sentinel-2 data. A LinkNet deep neural network with a DenseNet encoder was used to solve the problem of automated rill erosion mapping. The recognition results for the study area of more than 345,000 sq. km were summarized to a grid of 3037 basins and analyzed to assess the relationship with the main natural-anthropogenic factors. Generalized additive models (GAM) were used to model the dependency of rill erosion density to explore complex relationships. A complex nonlinear relationship between erosion processes and topographic, meteorological, geomorphological, and anthropogenic factors was shown.


Introduction
Soil erosion became the main factor in degrading a fertile layer of agricultural lands. Intensive land use combined with natural factors creates conditions for increased development of soil erosion, including rills, erosion furrows, and ephemeral gullies. Often, these erosion forms turn into permanent gullies, completely removing the area from agricultural turnover. Therefore, it is not surprising that researchers worldwide pay special attention to the study of soil erosion and, in particular, rill erosion.
Both field and mathematical methods are the main methods for studying soil erosion. The field method is very accurate in estimating the volume of erosion changes on local areas, which somewhat complicates the spatial interpretation of the results at the region or landscape levels. Such methods include the classical methods of reference sites, landmarks, microprofiling, and measuring the length and width of washouts with a tape measure. Reference sites allow estimating the direct volume of soil washed away from the territory, providing an opportunity to artificially change the conditions of the "environment" [1]. Based on the data from such sites, the data for creating mathematical models of soil erosion assessment, which will be discussed below, were obtained. The method of landmarks and microprofiling makes it possible to accurately estimate the depth of washout in representative sites [2]. Studies using these techniques helped find a strong correlation between the width of washout and its depth, which, in turn, estimates the volumes of soil washed away from the territory [3,4]. At present, geodetic methods based on relief reconstruction using scanning data or photogrammetry, both ground and airborne, are being developed to assess erosion. The most accurate results can be achieved using ground-based laser scanning data not been carried out directly based on them before. Such solutions allow for achieving the best accuracy but are labor consuming and low productive. To solve this problem, the development of approaches based on the automation of rills washout extraction is necessary, which is the goal of this work. The source data used as a cartographic basis sets a list of possible approaches that can be applied to solve the purpose. The simplest of them is object recognition based on a thresholding approach, where the threshold defines the boundary of the reflectivity of the spectral data characteristic of the study object. Such methods are suitable for identifying different land-use types [43,44]. However, the result will be unpredictable when recognizing rills using a similar approach on different types of soils. Machine learning methods can be used to consider the spatial variability of environmental factors, such as the highly proven Random Forest or Support Vector Machine method [45,46]. However, such approaches only give a probability of soil erosion in a particular pixel, while not distinguishing the erosions themselves into a "tree-like pattern" in the landscape. In addition, the methods are very sensitive to the amount of input data-the more information used for analysis, the more consistent the results will be. On small watersheds, such approaches can be successfully applied; for large areas, their applicability is questionable.
Recently, there has been a rapid increase in the number of works related to deep neural networks (DNNs) for the semantic segmentation of remote sensing data. This has been facilitated by the increased quality of remote sensing data, user-friendly deep-learning framework availability, and the increased computing power available to researchers. At present, DNNs enable successful automated interpretation of anthropogenic objects [47][48][49], shorelines [50,51], land use [52][53][54], vegetation cover dynamics [55,56], and exogenous processes [57,58]. In all cases, the authors note the higher recognition accuracy of the objects of interest than other methods and emphasize scaling the trained models. Deep neural networks have not been applied to the recognition and mapping of rill erosion. Nevertheless, the importance of the problem under study and the promising use of artificial intelligence for solving this problem determines the necessity of developing an appropriate methodology.
Given the intensification of rill erosion, the necessity of mapping it, and the limitations of field survey and modeling, this study aims to develop and apply a deep neural networkbased methodology for automated rill erosion detection from remote sensing data.

Study Area
The study area is the Middle Volga region of Russia. The territory is located in the central part of the East European Plain [59]. The Middle Volga region's agro-industrial complex (AIC) is of all-Russian importance. The region holds one of the top places in the country to produce grain, including the valuable grain crop-wheat. Agriculture is characterized by high efficiency due to very favorable natural conditions. At the same time, the agricultural sector of the territory has a substantial impact on natural-territorial complexes, primarily on the soil cover. Given the intensity of plowing, cattle grazing, and natural conditions conducive to the development of exogenous processes, this area is highly affected by soil erosion. This situation is not new; it is not accidental that the territory of the Middle Volga region is historically called an erosion belt. At the same time, anti-erosion measures in recent decades are of a point and haphazard nature, mainly due to the poor study of the rate and development direction of the process, primarily rill erosion.
The study area is 345,477 sq. km. The macroregion includes Mari El, Tatarstan, Chuvashia, Saratov, Samara, Ulyanovsk, and Penza regions. The Middle Volga region is located in the forest landscape zone in the north, the forest-steppe landscape zone in the center, and the steppe landscape zone in the south ( Figure 1).
The relief of the Middle Volga region has a spotted asymmetry of slopes: the average height-146 m, maximum elevation-316.8 m, and minimum-13.2 m [60]. The region's climate is moderately continental and continental in the south, varying wildly from region to region. Average annual precipitation ranges from 471 to 549 mm in the north to 264-424 mm in the south.
The steepness of the slopes of the territory, in general, is favorable for agriculture; a quarter of the region is steep to 1 degree, and a significant part of the territory is characterized by steepness in the range of 1-2.5 degrees, rising to 5 degrees in the area of Bugulminsko-Belebeyevskaya upland [60]. Erosion risks are not uniform and vary from region to region [61]. Soil-forming rocks are not homogeneous and are mainly clayey and heavy loam. There are sandy soils and medium-loam in the north of the Republic of Mari El region and light loam in the forest-steppe zone (Penza Oblast). Soils vary; chernozems of various subtypes predominate, mostly leached (Haplic Chernozems), as well as ordinary (Calcic Chernozems), typical (also Haplic Chernozems), and southern (Haplic Chernozems too)(steppe soils). The remaining types have a clearly expressed zonal pattern, where the gray forests (Haplic Greyzems) are typical for the forest territories and north of the forest-steppe zone [60].

Materials and Methods
The development and implementation of a methodology for automatic rill erosion recognition involve several interrelated points:

1.
Selection and preparation of remote sensing data from space; 2.
Preparation of training and test sets; 3.
Training the neural network and evaluating the quality of rill recognition; 4.
Implementation of the neural network for the entire study area and vectorization of recognition results; 5.
Calculation of rill erosion length as a measure of rill erosion density index in the basins; 6.
Analysis of the obtained results.

Remote Sensing Data
Data from the Sentinel-2 satellite constellation was used for recognition. Unlike other sources of multispectral remote sensing data, such as Landsat, MODIS, SPOT, RapidEye and so on, Sentinel-2 imagery provides the best resolution in the target bands (10 m) at no charge. The cloudless composite data of the Sentinel-2 mission for the spring period, April-June, were used as input data to develop the rill erosion recognition technique. Near-infrared images with baseline resolution (10 m) were used for recognition. The composite was created using the Google Earth Engine (GEE) [62], product "Sentinel-2 MSI: MultiSpectral Instrument, Level-2A" (COPERNICUS/S2_SR). GEE allows the processing of big remote sensing data and facilitates some routine operations. A total of 2988 scenes from both S2A and S2B satellites were used to create the composite. To create the composite, the 2018-2019 images were filtered by date, images were cleaned from clouds and shadows, median pixel brightness values were calculated, and imagery was clipped to the boundaries of the study area and reprojected into the WGS 84/Pseudo-Mercator (EPSG:3857) projection coordinate system. A cloud probability raster was used to clear the imagery from the clouds. The S2 cloud probability is created using the sentinel2-cloud-detector library (LightGBM [63]). All bands are upsampled using bilinear interpolation to 10 m resolution before the gradient boost base algorithm is applied. The resulting 0 to 1 floating point probability is scaled to 0-100 and stored as a UINT8. Areas missing any or all of the bands are masked out. Higher values are more likely to be clouds or highly reflective surfaces (e.g., roof tops or snow). A 15% probability of having clouds on the scene was used as the cloud threshold.

Preparation of the Training Dataset
In the resulting spring composite, an area with visually the densest rill network and the presence of major land-use classes of water, anthropogenic, forest, and agricultural features was selected and clipped. For this area of more than 2500 square kilometers, continuous manual rill erosion recognition was performed to create a set of reference data samples (Ground truth, GT) ( Figure 2). Manual recognition was made by QGIS tools using Sentinel-2 data and ultra-high resolution (UHR) images (Bing, Google, and so on, available as a substrate at no charge). The main selection was made on the S2 images. In controversial cases, the UHR images were used. They were also used for the non-selection of permanent gullies. Due to the difficulty of classifying erosion forms on such a vast territory with different natural conditions, the differentiation of erosion forms was not carried out. Ephemeral gullies" and "rills" were defined as all erosion forms that can be interpreted on satellite images but do not belong to permanent gullies. Considering the resolution of the satellite image, rills are the most minor erosion form, which transforms into an ephemeral gully. The resulting sample was rasterized and reduced to the resolution of the satellite image fragment, and then both images were cut into 256×256 pixels patches. In total, 10,933 satellite image-binary mask pairs have been obtained this way ( Figure 3). The resulting rasters were further randomly transformed to augment the number of rasters by three times artificially. The resulting dataset was divided into training and test sets in the ratio of 1/5.   The resulting sample was rasterized and reduced to the resolution of the satellite image fragment, and then both images were cut into 256 × 256 pixels patches. In total, 10,933 satellite image-binary mask pairs have been obtained this way ( Figure 3). The resulting rasters were further randomly transformed to augment the number of rasters by three times artificially. The resulting dataset was divided into training and test sets in the ratio of 1/5. The resulting sample was rasterized and reduced to the resolution of the satellite image fragment, and then both images were cut into 256×256 pixels patches. In total, 10,933 satellite image-binary mask pairs have been obtained this way ( Figure 3). The resulting rasters were further randomly transformed to augment the number of rasters by three times artificially. The resulting dataset was divided into training and test sets in the ratio of 1/5.

Training a Deep Neural Network
A convolutional neural network is a class of artificial neural network that uses convolutional layers to filter inputs for useful information. The convolution operation involves combining input data (feature map) with a convolution kernel (filter) to form a transformed feature map. The filters in the convolutional layers (conv layers) are modified based on learned parameters to extract the most useful information for a specific task.
Convolutional networks adjust automatically to find the best feature based on the task. Applications of convolutional neural networks include various image (image recognition, image classification, video labeling, text analysis) and speech (speech recognition, natural language processing, and text classification) processing systems, along with state-of-the-art AI systems, such as robots, virtual assistants, and self-driving cars.
LinkNet (Figure 4), a fully convolutional neural network for semantic image segmentation, was chosen as a neural network architecture [64]. In contrast to the well-known and well-proven U-NET architecture [57], which has been used for semantic segmentation of satellite images in geomorphological studies [65], LinkNet uses residual blocks instead of convolutional blocks in its encoder and decoder networks. The LinkNet deep neural network architecture efficiently exchanges the information received by the encoder with the decoder after each downsampling block. This turns out to be better than using pooling indices in the decoder, or just using fully convolutional networks in the decoder. This feature transfer technique not only gives good accuracy, but also allows a small number of parameters to be used in the decoder.

Training a Deep Neural Network
A convolutional neural network is a class of artificial neural network that uses convolutional layers to filter inputs for useful information. The convolution operation involves combining input data (feature map) with a convolution kernel (filter) to form a transformed feature map. The filters in the convolutional layers (conv layers) are modified based on learned parameters to extract the most useful information for a specific task. Convolutional networks adjust automatically to find the best feature based on the task. Applications of convolutional neural networks include various image (image recognition, image classification, video labeling, text analysis) and speech (speech recognition, natural language processing, and text classification) processing systems, along with state-of-theart AI systems, such as robots, virtual assistants, and self-driving cars.
LinkNet (Figure 4), a fully convolutional neural network for semantic image segmentation, was chosen as a neural network architecture [64]. In contrast to the well-known and well-proven U-NET architecture [57], which has been used for semantic segmentation of satellite images in geomorphological studies [65], LinkNet uses residual blocks instead of convolutional blocks in its encoder and decoder networks. The LinkNet deep neural network architecture efficiently exchanges the information received by the encoder with the decoder after each downsampling block. This turns out to be better than using pooling indices in the decoder, or just using fully convolutional networks in the decoder. This feature transfer technique not only gives good accuracy, but also allows a small number of parameters to be used in the decoder. The initial block contains a convolution layer with a kernel size of 7 × 7 and a stride of 2, then a max-pool layer with a window size of 2 × 2 and a stride of 2. Similarly, the last block performs a full convolution, taking feature maps from 64 to 32, then a 2D convolution. Finally, a full convolution is used as a classifier with a kernel size of 2 × 2. As practice has shown, this allows for better semantic aggregation, including dendro-like constructions (tree-like patterns on the landscape) and rill networks. By trial-and-error method, it became clear that the best results in decoding rills can be achieved using transfer learning, a deep neural network learning method that allows using the knowledge obtained about one deep learning problem and applying it to another, but with a similar problem. The encoders used for DenseNet image classification [66] were used in our case. In contrast to The initial block contains a convolution layer with a kernel size of 7 × 7 and a stride of 2, then a max-pool layer with a window size of 2 × 2 and a stride of 2. Similarly, the last block performs a full convolution, taking feature maps from 64 to 32, then a 2D convolution. Finally, a full convolution is used as a classifier with a kernel size of 2 × 2. As practice has shown, this allows for better semantic aggregation, including dendro-like constructions (tree-like patterns on the landscape) and rill networks. By trial-and-error method, it became clear that the best results in decoding rills can be achieved using transfer learning, a deep neural network learning method that allows using the knowledge obtained about one deep learning problem and applying it to another, but with a similar problem. The encoders used for DenseNet image classification [66] were used in our case. In contrast to similar models, the features are not summed up but concatenated (channelwise concatenation) into a single tensor before passing them to the next layer. The model learning and implementation algorithm was performed in the Python 3.7 programming environment using the Tensorflow library. The input of the neural network was a stack of image-mask pairs prepared in the previous step. EarlyStopping monitoring was used to prevent overtraining of the model, and the IOU Score metric, Jaccard's coefficient, was used to check the model's learning capability [67].

Statistical Analysis of the Obtained Results
Further research was conducted to analyze the influence of various natural-anthropogenic factors on the intensity of erosion processes in the study area. Classical correlation, regression analysis, and analysis of variance (AOV) (one-way [68] and multivariate [69]) were used to analyze the relationship between natural and anthropogenic factors on the density of the rill network. The AOV was performed using the Statsmodel package for Python 3. To analyze complex relationships, an analysis based on generalized additive models was used. The generalized additive model (GAM) is a generalized linear model in which a linear predictor depends linearly on unknown smooth functions of some predictor variables, and interest is focused on inference about these smooth functions. GAM uses spline functions, functions that can be combined to approximate arbitrary functions. GAM introduces penalties for the weights so that they remain close to zero, which effectively reduces the flexibility of the splines and reduces the possibility of overfitting. The smoothness parameter, which is normally used to control curve flexibility, is then adjusted using cross-validation. The natural-anthropogenic factors presented in Table 1 [60] were used for the analysis. Table 1. Natural-anthropogenic factors of rill erosion selected for analysis [60].

Results and Discussion
The trained RECNN (rill erosion convolutional neural network) was tested on an independent test dataset. The recognition accuracy was 0.62, F1-measure was 0.76, and loss-function was 0.27. Qualitatively analyzing the obtained evaluation results ( Figure 5), we can note a relatively high level of recognition of the rill erosion network. It is crucial that not a single case of detection of gullies or ground roads, which are abundant in the study area, instead of rills, was recorded.

Results and Discussion
The trained RECNN (rill erosion convolutional neural network) was tested on an independent test dataset. The recognition accuracy was 0.62, F1-measure was 0.76, and lossfunction was 0.27. Qualitatively analyzing the obtained evaluation results ( Figure 5), we can note a relatively high level of recognition of the rill erosion network. It is crucial that not a single case of detection of gullies or ground roads, which are abundant in the study area, instead of rills, was recorded. The trained and tested model was decided to apply to the entire study area ( Figure 6). For the obtained geometry, the length was calculated, which was aggregated to a grid of basins [60] (Figure 7). Analyzing the obtained maps, the presence of specific clusters becomes evident. The first cluster-the highest density of rill erosion-is located in the left-bank part of the Volga River in the Saratov region. The second cluster-the area of medium erosion-characterizes the right bank of the Volga in the Saratov region and Zakamye of the Republic of The trained and tested model was decided to apply to the entire study area ( Figure 6). For the obtained geometry, the length was calculated, which was aggregated to a grid of basins [60] (Figure 7). Analyzing the obtained maps, the presence of specific clusters becomes evident. The first cluster-the highest density of rill erosion-is located in the left-bank part of the Volga River in the Saratov region. The second cluster-the area of medium erosion-characterizes the right bank of the Volga in the Saratov region and Zakamye of the Republic of Tatarstan. The Republics of Mari El and Chuvashia territories are characterized mainly by low rill erosion, forming the third cluster. Tatarstan. The Republics of Mari El and Chuvashia territories are characterized mainly by low rill erosion, forming the third cluster.   Tatarstan. The Republics of Mari El and Chuvashia territories are characterized mainly by low rill erosion, forming the third cluster.   In total, information in the study area is aggregated on 3037 basins with an average area of 117.83 sq. km. The average density of rill erosion over the entire study area is 0.23 km/sq. km, and the maximum density is 3.79 km/sq. km. In general, the distribution of the total length of rill erosion in the study area is lognormal; that is, most of the measurements contain, in general, relatively small values of the total length (Figure 8). In total, information in the study area is aggregated on 3037 basins with an average area of 117.83 sq. km. The average density of rill erosion over the entire study area is 0.23 km/sq. km, and the maximum density is 3.79 km/sq. km. In general, the distribution of the total length of rill erosion in the study area is lognormal; that is, most of the measurements contain, in general, relatively small values of the total length ( Figure 8). For most of the selected factors, no strong correlations were found with the total length of rill erosion. Only some factors revealed statistically significant (p < 0.001) weak relationships (r > |0.3|). For example, an inverse correlation describes the relationship between total rill erosion length and forest cover, hydrothermal coefficient, mean annual precipitation, mean May-August precipitation, and mean warm-season precipitation. The inverse relationship with forest cover is more than obvious and only states a well-known fact-the less forest cover of the territory, the greater the risk of soil and gully erosion. The inverse relationship with the other indicators is not so obvious. Selyaninov's hydrothermal coefficient of wetness (GTK) is a characteristic of the level of moisture availability of the territory. It is widely used in agronomy for general climate assessment and allocation of zones with different levels of moisture availability to determine the usefulness of the cultivation of various crops. The factor is closely related to the amount of precipitation in the territory. The inverse relationship here is explained not so much by the very nature of the dependence "the more precipitation, the less erosion", which cannot be, as by the influence of climatic conditions in general. This is partly confirmed by a direct relationship between the total length of the rill network and such factors, such as the sum of active temperatures and the average annual maximum temperature for the year. Apparently, with increasing temperature and climate change due to global warming, the character of precipitation changes-they may fall rarely, but with great intensity [70][71][72][73][74]. Unfortunately, it is impossible to check this hypothesis directly due to the lack of data on mean annual precipitation intensity changes. However, the same conclusions are drawn by other researchers who have worked in the study area [31,32,75,76]. Nevertheless, there is still no clear answer to this question.
Another question is how classical methods of relationship analysis can assess the effect of a single factor on the dependent variable. Linear dependencies are rare in nature, and polynomial dependencies of high degrees are even rarer. The analysis demonstrated poor applicability of classic correlation and regression analysis for revealing the relationships between natural-anthropogenic factors and the length of the erosion network. Some examples (Figures 9-11) demonstrate that the dependencies are close to the logistic or exponential type but more complex. Therefore, it was decided to use generalized additive models (GAM) for dependence analysis [77]. These models allow estimating the For most of the selected factors, no strong correlations were found with the total length of rill erosion. Only some factors revealed statistically significant (p < 0.001) weak relationships (r > |0.3|). For example, an inverse correlation describes the relationship between total rill erosion length and forest cover, hydrothermal coefficient, mean annual precipitation, mean May-August precipitation, and mean warm-season precipitation. The inverse relationship with forest cover is more than obvious and only states a well-known fact-the less forest cover of the territory, the greater the risk of soil and gully erosion. The inverse relationship with the other indicators is not so obvious. Selyaninov's hydrothermal coefficient of wetness (GTK) is a characteristic of the level of moisture availability of the territory. It is widely used in agronomy for general climate assessment and allocation of zones with different levels of moisture availability to determine the usefulness of the cultivation of various crops. The factor is closely related to the amount of precipitation in the territory. The inverse relationship here is explained not so much by the very nature of the dependence "the more precipitation, the less erosion", which cannot be, as by the influence of climatic conditions in general. This is partly confirmed by a direct relationship between the total length of the rill network and such factors, such as the sum of active temperatures and the average annual maximum temperature for the year. Apparently, with increasing temperature and climate change due to global warming, the character of precipitation changes-they may fall rarely, but with great intensity [70][71][72][73][74]. Unfortunately, it is impossible to check this hypothesis directly due to the lack of data on mean annual precipitation intensity changes. However, the same conclusions are drawn by other researchers who have worked in the study area [31,32,75,76]. Nevertheless, there is still no clear answer to this question.
Another question is how classical methods of relationship analysis can assess the effect of a single factor on the dependent variable. Linear dependencies are rare in nature, and polynomial dependencies of high degrees are even rarer. The analysis demonstrated poor applicability of classic correlation and regression analysis for revealing the relationships between natural-anthropogenic factors and the length of the erosion network. Some examples (Figures 9-11) demonstrate that the dependencies are close to the logistic or exponential type but more complex. Therefore, it was decided to use generalized additive models (GAM) for dependence analysis [77]. These models allow estimating the dependencies in the form of spline functions, which allows one to make better approximations and comprehensively study the dependence of a factor on a variable [24,78]. dependencies in the form of spline functions, which allows one to make better approximations and comprehensively study the dependence of a factor on a variable [24,78].   dependencies in the form of spline functions, which allows one to make better approximations and comprehensively study the dependence of a factor on a variable [24,78].   Modeling of generalized relationships was performed in a package for building generalized additive models in the Python language [79]. During modeling, statistically significant (p < 0.001), stable, and strong relationships were found between the density of rill erosion and such natural factors as mean basin elevation ( Figure 12A), mean basin slope ( Figure 12B), basin elevation range ( Figure 13A), basin forest cover ( Figure 13B), average air temperature in January ( Figure 14A), average annual long-term maximum temperature ( Figure 14B), average warmseason precipitation ( Figure 15A), average basin hydrothermal coefficient ( Figure 15B), longitude ( Figure 16A  Modeling of generalized relationships was performed in a package for building generalized additive models in the Python language [79]. During modeling, statistically significant (p < 0.001), stable, and strong relationships were found between the density of rill erosion and such natural factors as mean basin elevation ( Figure 12A), mean basin slope ( Figure 12B), basin elevation range ( Figure 13A), basin forest cover ( Figure 13B), average air temperature in January ( Figure 14A), average annual long-term maximum temperature ( Figure 14B  Modeling of generalized relationships was performed in a package for building generalized additive models in the Python language [79]. During modeling, statistically significant (p < 0.001), stable, and strong relationships were found between the density of rill erosion and such natural factors as mean basin elevation ( Figure 12A), mean basin slope ( Figure 12B), basin elevation range ( Figure 13A), basin forest cover ( Figure 13B), average air temperature in January ( Figure 14A), average annual long-term maximum temperature ( Figure 14B), average warmseason precipitation ( Figure 15A), average basin hydrothermal coefficient ( Figure 15B    Generalized models confirmed the known relationships between soil erosion and natural factors and revealed complex nonlinear relationships. It is observed that the density of the erosion network is stable with an increase in the average elevation of the basin, but up to 100 m, after which the growth of density is replaced by a decrease in the density and its stabilization ( Figure 12A). In our opinion, this dependence should be considered in combination with the effect of average slope steepness on the density of the erosion network. In this case, the relationship is more superficial and closer to a direct linear relationship-the higher the steepness, the greater the density ( Figure 12B). However, when analyzing steepness in basins with an average elevation of less than 100 m, a pattern becomes apparent-steepness does not exceed 3 degrees, and the slopes of these basins are  Generalized models confirmed the known relationships between soil erosion and natural factors and revealed complex nonlinear relationships. It is observed that the density of the erosion network is stable with an increase in the average elevation of the basin, but up to 100 m, after which the growth of density is replaced by a decrease in the density and its stabilization ( Figure 12A). In our opinion, this dependence should be considered in combination with the effect of average slope steepness on the density of the erosion network. In this case, the relationship is more superficial and closer to a direct linear relationship-the higher the steepness, the greater the density ( Figure 12B). However, when analyzing steepness in basins with an average elevation of less than 100 m, a pattern becomes apparent-steepness does not exceed 3 degrees, and the slopes of these basins are Generalized models confirmed the known relationships between soil erosion and natural factors and revealed complex nonlinear relationships. It is observed that the density of the erosion network is stable with an increase in the average elevation of the basin, but up to 100 m, after which the growth of density is replaced by a decrease in the density and its stabilization ( Figure 12A). In our opinion, this dependence should be considered in combination with the effect of average slope steepness on the density of the erosion network. In this case, the relationship is more superficial and closer to a direct linear relationship-the higher the steepness, the greater the density ( Figure 12B). However, when analyzing steepness in basins with an average elevation of less than 100 m, a pattern becomes apparent-steepness does not exceed 3 degrees, and the slopes of these basins are well suited for farming, particularly as arable land. Indeed, in these areas, the percentage of the plowed area exceeds 50%.
The effect of relief morphometry is also confirmed by the dependence of the density of rill erosion on the height dispersion index in the basin-the difference between the maximum and minimum heights. In general, the smaller the height difference, the lower the density, but fundamentally, the most significant decrease in rill erosion intensity is observed in basins with a small height difference (up to~70 m), after which the intensity of density reduction is not as pronounced and essentially remains constant ( Figure 13A). The influence of anthropogenic development on the intensity of erosion processes is also confirmed by the inverse relationship between forest cover and the density of the rill network ( Figure 13B).
Climatic factors have no less influence on the character of erosion processes in the study area, partly complementing the topographic and economic predisposition of the area to the development of rills. The relationship between the density of rill erosion and the average air temperature in January is interesting-the lower the mean annual temperature in the coldest month of the year, the more intense the process's pattern is ( Figure 14A). It seems to be related to the soil freezing capacity, which ensures the intensity of the erosion processes. However, it is impossible to verify this reliably due to the lack of adequate models of soil freezing in the study area. Temperature influences the character of erosion processes in summer as well-there is a strong direct correlation between the intensity of the process and the maximum air temperature starting from 36 degrees Celsius ( Figure 14B). Initially, the authors attributed this to the potentially high evaporation capacity at such temperatures and the possible high precipitation intensity. However, analysis of the dependence of erosion processes on the precipitation layer in the warm season and the zoning map of the study area by the maximum temperature rejected this hypothesis. The maximum air temperatures are associated with erosion processes indirectly through the prevailing soil types, as the conditions for the emergence of one or another type. This is confirmed because the areas with maximum temperatures over 36 degrees are located in the transition zone of semi-deserts and dry-steppes with prevailing chestnut and solonetz soil types. These types of soils are characterized by a large number of cracks reaching a depth of 5-6 cm [80]. Even with an aggregate not large annual precipitation layer, the passing rains create a continuous surface runoff on the dry soils, which erodes the soils along the available cracks, forming stable washouts.
In general, the relationship between the intensity of rill erosion and moisture is direct; there is a steady trend to an increase in the density of the erosion network with increasing precipitation in the warm period up to~270 mm, after which the role of precipitation decreases while maintaining the rule-"much precipitation-more intense erosion" ( Figure 15A). Stabilization of the trend, in this case, is explained by the fact that 270 mm is a kind of boundary between steppe and forest-steppe zones, in which there is a clear seasonality of precipitation. By the time of the most intensive precipitation (summer-autumn period), vegetation with good soil-protective properties develops on these territories [24], which becomes the dominant factor of soil erosion intensification. The influence of moisture factors most clearly shows the relationship between the density of rill erosion and the value of the hydrothermal coefficient of the territory ( Figure 15B). The density of the erosion network decreases with increasing GTK up to the value of 0.7-a clear border between arid area and wetted area. The intensity of soil-protecting vegetation increases with the increase of the SCC, but the precipitation becomes sufficient for the vegetation to perform its soil-protecting functions less effectively when moving into the wetted zone.
The analysis of the influence of the geographical location of the basin on the intensity of erosion processes also demonstrated statistically significant relationships. The density of rill erosion decreases linearly as one move to the east of the study area ( Figure 16A), which is explained by the fact that the western part of the study area has a higher right bank of the Volga River and, as a result, higher values of slope steepness, which, as shown earlier, directly affect the intensity of the process. The latitudinal pattern is not so linear and more complex ( Figure 16B)-with moving from the south to the north up to the zone of typical steppes, there is a systematic increase in the density of rill erosion, after which there is some decrease in the intensity of soil erosion and its stabilization, after which, starting from the northern part of the forest-steppe zone there is another gradual increase in the intensity of erosion activity. The reasons for this, as previously shown, lie in the zonal features of the study area, changes in the degree of moisture, as well as the types of plant communities, including agricultural ones, growing in a particular area.
It is worth separately considering the influence of soil-geological conditions of erosion processes. A one-factor analysis of variance was conducted to analyze the relationship between the prevailing soil type and total rill erosion length ( Figure 17). study area ( Figure 16A), which is explained by the fact that the western part of the study area has a higher right bank of the Volga River and, as a result, higher values of slope steepness, which, as shown earlier, directly affect the intensity of the process. The latitudinal pattern is not so linear and more complex ( Figure 16B)-with moving from the south to the north up to the zone of typical steppes, there is a systematic increase in the density of rill erosion, after which there is some decrease in the intensity of soil erosion and its stabilization, after which, starting from the northern part of the forest-steppe zone there is another gradual increase in the intensity of erosion activity. The reasons for this, as previously shown, lie in the zonal features of the study area, changes in the degree of moisture, as well as the types of plant communities, including agricultural ones, growing in a particular area.
It is worth separately considering the influence of soil-geological conditions of erosion processes. A one-factor analysis of variance was conducted to analyze the relationship between the prevailing soil type and total rill erosion length ( Figure 17). The basic soils of taiga and coniferous-broadleaved forests, as well as hydromorphic soils, are least susceptible to erosion activity. The smallest total length of rills is observed on sod-podzol soils (shallow-podzol and deep-podzol). These taiga soils, as a rule, are covered by mixed forests with a dense canopy and good protection against the formation of surface runoff. Moving southward, there is a slight increase in the total length of the rill network on gray forest soils of various subtypes. These soils are actively involved in agricultural management and have suitable hydrometeorological conditions. Weak susceptibility to erosion of territories with hydromorphic soils is, in our opinion, evident and connected, by and large, with land use on such soils-these are meadows with maximum soil-protective properties. The next group of soils is more affected by erosion processessteppe soils. These are all kinds of chernozems, starting from leached chernozems. Southern chernozems, which are peculiar "record-breakers" by erosion susceptibility among The basic soils of taiga and coniferous-broadleaved forests, as well as hydromorphic soils, are least susceptible to erosion activity. The smallest total length of rills is observed on sod-podzol soils (shallow-podzol and deep-podzol). These taiga soils, as a rule, are covered by mixed forests with a dense canopy and good protection against the formation of surface runoff. Moving southward, there is a slight increase in the total length of the rill network on gray forest soils of various subtypes. These soils are actively involved in agricultural management and have suitable hydrometeorological conditions. Weak susceptibility to erosion of territories with hydromorphic soils is, in our opinion, evident and connected, by and large, with land use on such soils-these are meadows with maximum soil-protective properties. The next group of soils is more affected by erosion processes-steppe soils. These are all kinds of chernozems, starting from leached chernozems. Southern chernozems, which are peculiar "record-breakers" by erosion susceptibility among typical steppes soils, are distinguished. Southern chernozems are in the cluster with the highest density of erosion network (Figure 7). Southern chernozems are spread in the southern part of the steppe zone. They are formed in the conditions of semiarid climate under soddy-cereal medium steppes. Grass cover is sparse, clearly expressed summer semi-rest period for most dominant cereals. Soil-forming rocks are represented mainly by loess and loess-like loams, often containing easily soluble salts, as well as eluvial-deluvial deposits. Water conditions of soils are solid. Agricultural development of southern chernozems is high: in the European part of Russia, it exceeds 50%, while moving to the east, plowing decreases and the number of pastures increases. The main crops grown are cereals (wheat and corn) and legumes; considerable areas are occupied by industrial crops (sugar beet and tobacco), vegetable, and melon crops. Plowed soils are prone to water and wind erosion, structural degradation, and slitization under irrigation. The last group is soils with a maximum total length of erosion network-soils of dry steppes and semi-deserts. These are different subtypes of chestnut soils and solonetz soils, a description of which was given earlier.
The one-factor analysis of variance was also applied to analyze the relationship between erosion activity and the predominant type of soil-forming rocks ( Figure 18). In this case, grouping types is problematic due to implicit clustering; however, certain patterns are still observed. The increase in total length of rill erosion directly depends on clay content in the soil-forming rock and-in its water-holding capacity, drainability. The minimum, near-zero erosion network length is observed on weakly erodible shales. Analyzing the rocks on which erosion is observed, the minimum length of the erosion network falls on the group of rocks conditionally united into the so-called "sands"-these are sandy rocks and sandstones. These rocks have the maximum infiltration capacity, not allowing the formation of washout surface runoff and coming to the soils with minimal humus horizons. Slightly better eroded areas are composed of clayey and loamy rocks, underlain by sandy and sandy loam rocks, as well as light loam rocks themselves. Medium loamy rocks and rocks of different textures with a predominance of loam and clay are better eroded. On average, easily erodible chalky rocks-limestone and other carbonate rocks-are even better erodible. They usually form gray forest soils and leached chernozems, heavily involved in agricultural production. The most significant total length of the erosion network is expected to be in water-resistant clay and heavy loam rocks, on which surface runoff is formed, eroding fertile soil layers. Separately, sandy loam rocks, for which we could not find a reliable relationship-these rocks, draining well the surface runoff, can be resistant to erosion processes and, folding typical chernozems, due to active economic activity falls in the area of maximum rill erosion. The analysis of variance was conducted separately by landscape zones, and the result can be observed in Figure 19. The results confirm the earlier conclusions-the zones of dry and typical steppes and semi-deserts are most susceptible to erosion processes. The analysis of variance was conducted separately by landscape zones, and the result can be observed in Figure 19. The results confirm the earlier conclusions-the zones of dry and typical steppes and semi-deserts are most susceptible to erosion processes. The analysis of variance was conducted separately by landscape zones, and the result can be observed in Figure 19. The results confirm the earlier conclusions-the zones of dry and typical steppes and semi-deserts are most susceptible to erosion processes. Figure 19. Relationship between the rill erosion network length and landscape zone. Figure 19. Relationship between the rill erosion network length and landscape zone.
The analysis of the relationship between erosion processes and anthropogenic load also confirms the previous conclusions ( Figure 20). The score of the anthropogenic load is given based on the development of the territory, the presence of settlements, industry, and roads of various types, as well as population density. The analysis of the relationship between erosion processes and anthropogenic load also confirms the previous conclusions ( Figure 20). The score of the anthropogenic load is given based on the development of the territory, the presence of settlements, industry, and roads of various types, as well as population density. The intensity of erosion processes increases with increasing anthropogenic load, but not linearly. The maximum intensity of erosion is observed under moderate anthropogenic load. Decrease of the total length of erosion forms with further increase of anthropogenic load is related to decrease of the area not occupied by cities, manufactories, and other anthropogenic objects. The intensity of erosion processes increases with increasing anthropogenic load, but not linearly. The maximum intensity of erosion is observed under moderate anthropogenic load. Decrease of the total length of erosion forms with further increase of anthropogenic load is related to decrease of the area not occupied by cities, manufactories, and other anthropogenic objects.
The cumulative effect of anthropogenic load, forest cover, prevailing soil type, and soil-forming rocks has the most significant contribution to the development of erosion processes in the study area ( Figure 21). The intensity of erosion processes increases with increasing anthropogenic load, but not linearly. The maximum intensity of erosion is observed under moderate anthropogenic load. Decrease of the total length of erosion forms with further increase of anthropogenic load is related to decrease of the area not occupied by cities, manufactories, and other anthropogenic objects.
The cumulative effect of anthropogenic load, forest cover, prevailing soil type, and soil-forming rocks has the most significant contribution to the development of erosion processes in the study area ( Figure 21).

Conclusions
Deep neural networks can effectively solve geomorphological problems, allowing large-scale mapping of erosion processes. For the first time in the world, a detailed assessment of the density of rill erosion and its relation to the main natural-anthropogenic factors is given. In this case, preparation of the training sample is the most time-consuming stage, taking 2/3 of the whole time of the study. The trained model allows highly efficient and at an acceptable level of recognition of the erosion network on high-resolution images, in our case using Sentinel-2 images. We investigated which Sentinel-2 bands allow the most efficient visual recognition at the preliminary stage. As a result of numerous experiments, the best "readability" of the image was demonstrated by a separate near-infrared band, winning in comparison with RGB-synthesis, synthesis of artificial colors in various combinations, and other channels, separately. Another sensitive point of the research is the choice of neural network architecture for semantic segmentation. Experiments were carried out with the most common architectures for such tasks, including the well-known U-Net architecture, as well as FPN, PSPNet, and LinkNet in combination with the most popular encoder groups, such as VGG, ResNet, SE-ResNet, ResNeXt, SE-ResNeXt, SENet154, DenseNet, Inception, MobileNet, and EfficientNet. The best recognition accuracy on the test set showed the FPN + EfficientNet combination; however, when applying the trained model, the recognition results were poor. The best results in terms of real-world applications were achieved using LinkNet + DenseNet (DenseNet201). In the future, it is planned to use the trained model to recognize the erosion network for 2019-2020 to assess the dynamics of erosion processes in time.
The analysis made it possible to assess the key factors influencing the intensity of rill erosion processes at a high level of confidence. In some cases, existing knowledge about the influence of certain factors on soil erosion was confirmed. Particular interest is the possibility of mapping the territories with the highest risk of erosion processes due to natural-anthropogenic conditions, which will allow producing precise and effective measures to minimize the risks. Under the conditions of intensive agricultural development, no reduction of anthropogenic load can be discussed; however, proper crop rotations, regular sowing of perennial grasses, and proper plowing will reduce the existing risks many times.
The use of generalized additive models (GAM) to analyze dependencies made it possible to describe complex relationships between the density of rill erosion and climate, geomorphological, and other geographic factors. Unfortunately, not all parameters could be used in the analysis; in the future, detailed modeling of the dependence of erosion processes on factors of soil freezing, and the intensity of precipitation, especially changes in intensity in time series, should be conducted.