Inﬂuence of DEM Elaboration Methods on the USLE Model Topographical Factor Parameter on Steep Slopes

: Runo ﬀ erosion is an important theme in hydrological investigations. Models assessing soil erosion are based on various algorithms that determine the relief coe ﬃ cient using rasterized digital elevation models (DEMs). For evaluation of soil loss, the most-used model worldwide is the USLE (Universal Soil Loss Equation), where the most essential part is the LS parameter, which is, in turn, generated from two parameters: L (slope length coe ﬃ cient) and S (slope inclination). The most signiﬁcant limitation of LS is the di ﬃ culty in obtaining the data needed to generate detailed DEMs. We investigated three popular data generation methods: aerial photographs (AP), aerial laser scanning (ALS), and terrestrial laser scanning (TLS) by assessing the quality and e ﬀ ect of DEMs generated from each method over an area of 40 m × 200 m in Silesia, Poland. Additionally, the relationship between particular LS USLE parameter components was carried out based on its ﬁnal distribution. Our results show that resolution strongly inﬂuences DEMs and the LS USLE parameters. We found a strong relationship between the degree of height data resolution and the accuracy level of the calculated parameters. Based on our investigations we conﬁrmed the highest inﬂuence on the LS USLE came from the S parameter. Additionally, we concluded that in examinations over large areas, terrestrial laser scanners are not ideal; the beneﬁts of their additional accuracy are outweighed by the additional time and labor consumption; in addition, terrestrial-based scans are sometimes not possible due to ground obstacles the limited scope of most lasers. Aerial photographs or point clouds generated by aerial laser scanners are su ﬃ cient for most purposes connected with surface ﬂow, and further developments can be based on the use of these techniques for obtaining ground information for modeling erosion processes. and and writing—original P.K. and writing—review P.K. published manuscript.

equation [28] for S factor. Many researchers worldwide use algorithms proposed by Desmet and Govers based on raster resolution and unit area [31]. This solution allows the model to be adapted to varied reliefs and includes the influence of flow concentration [31,32]. The formula introduced by Nearing [33] can be used to estimate Winchell et al. [34]. Rodriguez and Suarez [35] have modified the equation by introducing and comparing various variants of the GIS technique approach. Panagos et al. [27] proposed a new approach for topographic parameter determination on areas approaching the size of European countries using the UCA method. Panagos et al. [27] carried out analysis of influence of DEM resolution on distribution of the LS factor. In the work they used two resolutions: 25-m and 100-m. The advantage of the UCA method, emphasized by the authors of this paper, is the improved ease of calculating LS. Hickey [36] and van Remortel et al. [23,37] describe new models for identifying gaps in LS that include changes in turning point of inclination and channels based on the definition of LS [17]. These methods can overcome some disadvantages of the UCA method [38], which does not, for example, take into account channels. These methods still have limitations; they use single flow direction (SFD) algorithms [39] for calculating the length of the slope, thus generating only partial results in DEMs [40]. Moore and Burch [29,41] noticed that higher indicators of erosion and deposition occurred in basin convergence, as was postulated by the USLE/RUSLE and the Chinese Soil Loss Equation (CSLE). Multiple Flow Direction (MFD) algorithms can cover both convergent and divergent flows and work better than the SFD algorithms for real areas [21,42]. Spatial exactness of slope length and the LS are necessary to estimate water erosion [43].
Most of the accessible LS algorithms are already built into GIS software programs, such as IDRISI, SAGA GIS, GRASS, and ArcGIS, and the raster calculator in ArcGIS ESRI. DEMs have become popular and are very convenient for generating representations of the constantly changing topographic surface of the Earth. In the literature, more and more papers present results of investigations carried out using models integrated with GIS techniques. Due to the spatial distribution of individual factors, both erosion processes and sediment transport are differentiated by the basin area, and GIS tools enable division of the total area into small, approximately homogenous areas with similar properties and relatively uniform precipitation distribution in the shape of network cells [44].
Some papers using USLE models with GIS include Jain et al. [45] who investigated erosion quantity in the sub-basin of the Sitlarao River, a 52 km 2 area located in India, and used the McCool equation [24] to determine L, taking the dimension of the network raster cell into consideration.
Lee [46] used the USLE model to evaluate erosion intensity in the Boun region of South Korea, an area of 68.43 km 2 . He used a topographic map with a scale of 1:5000, a soil map with the scale 1:25,000, and the program Landsat Thematic Mapper to display pictures with 30 m resolution to determine land cover. USLE model parameters were calculated and collected from space databases and converted to a network with 5 m resolution using ARC Info 8.1 software. LS was determined according to the method described by Moore and Burch [29] and Moore and Wilson [28].
Yoshino and Ishioka [47] used the USLE model to calculate the erosion threat in the Cidanau River basin in Indonesia, an area of 223.17 km 2 in area, and LS was calculated based on a DEM of 30 m pixels.
Bhattarai and Dutta [48] used a DEM resolution of 90 m for calculating erosion intensity in the Mun River basin in Thailand with an area of 69,000 km 2 . L and S were determined in the following cells. Lee and Choi [49] used the model for calculation of erosion intensity in the 273 km 2 area of the Bosung basin in South Korea. For L and S assessment, they based their calculations the following raster cells, using a 20 m resolution DEM generated from a topographic map in the scale 1:5000. the parameters: S based on the Nearing [33] and L for the following raster cells.
Zhang [50] used the USLE model for comparison of erosion intensity as calculated according to the modified ER-USLE method, in the Shangshe River basin in China, with an area of 7400 km 2 . Inclination and slope length were determined from a DEM with a 1 m resolution network. Perovic et al. [51] carried out investigations in the Nisava River basin, located in Serbia, over a 2.85 km 2 area. They used GIS techniques to assess spatial distribution and identify areas particularly threatened by water erosion based on 30 m resolution and DEM created out of 1:50,000 topographic maps.
Chen et al. [52] used the RUSLE model to present the distribution of water erosion threats in the Miyun River basin, a 47.5 km 2 area located in southern part of China. They used artificial neural networks (ANN) and a 30 m resolution DEM.
Mhangara et al. [53] used the RUSLE model to evaluate erosion risk in the 745 km 2 Keiskamma River basin located in the Republic of South Africa. To determine the LS parameter, they used the SA-TEEC GIS module in the ArcView program. Prasannakumar et al. [54] used this same model to illustrate the spatial distribution of water erosion risk in the Siruvani River basin, a 205.54 km 2 area located in India. The LS coefficient was determined using the ArcInfo ArcGIS program and a DEM. Griffin et al. [30] L = (m + 1)· As Desmet and Govers [31] L x m ·D m+2 ·(22.13) m D is the raster resolution; A(i,j) is the unit area at the entrance to cell (i,j); m is the exponent of slope length indicator; and x is the coefficient correcting the path length of flow through raster cells depending on flow direction and calculated based on exposure Nearing [33] S = −1.5 +

17
(1+e (2.3−6.1sinq) ) S as above; D is the grid cell size in meters; x (i,j) = sin a (i,j) + cos a (i,j) ; a i is the aspect direction of the grid cell (i,j); and m is related to the F ratio of the rill to inter-rill erosion Yoshino and Ishioka [47] LS = l 22.13 m · 0.065 + 0.045·s + 0.0065·s 2 l is the slope length (m); s is slope %; m is dependent on slope %, where: m = 0.5 for slope s ≤ 5%; m = 0.4 (3.5 < s < 4.5%); m = 0.3 (1% < s < 3%); and m = 0.2 (1.8% < s). The slope length l was defined as the length of a slope with the greatest incline in a given pixel L, m, S, β as above θ is slope inclination in degrees ( • ) Lee and Choi [49] L i = L, m,β as above x is raster cell length; and θ is slope inclination in degrees ( • ) Perovic et al. [51] LS(r) = (m + 1)·[A(ra 0 )] m ·(sinb(rb 0 )) n A is the supply area outflowing to a cell [m]; b is slope inclination in degrees ( • ); m and n are parameters (m = 0.4; n = 1.4); 0 is determined from USLE with a length of 22.1 m; and b 0 is the parameter arising from the USLE model and equal to 0.09. χ is the accumulation of surface flow calculated from a DEM using the delinearization module of the basin in the Arc View 9.2 program; λ is the dimension of the raster cell; and θ is the inclination in degrees ( • ) Kumar and Kushwaha [55] used the RUSLE-3D model when they investigated the spatial distribution of soil loss in the 44 km 2 subbasin of the Pathri Rao River located in India. The modified 3-D component in the RUSLE model replaces slope length with the coefficient of a 20 m resolution network cell divided by the up-slope contributing area. Saygm et al. [56] used the combined RUSLE/SDR method for calculating the volume of material supplied from the 2.62 km 2 basin the Saraykoy II Reservoir in Turkey.
Taking into account the development of new measurement techniques and the almost unlimited accessibility of DEM data with various spatial resolutions (cell dimensions), verification is needed. The primary purpose of examining these works was a comparative analysis of the LS coefficient, as determined by using DEMs originating from three data sources of various resolutions. Development of new techniques results in greater integration of the model with geospatial data obtained by various measurement techniques. New tools make it possible to assess particular parameters with more precision.
In our work, we asked the following questions: Can new solutions for obtaining and generating data for spatial analysis (e.g., terrestrial laser scans) be used as a measurement technique for generating a greater number of points? Is it necessary to attain an equilibrium between the number of points and analyses assignments and can we obtain similar results using faster, cheaper, but less precise techniques? To answer these questions, we decided to analyse three methods for obtaining spatial information of land: Point cloud-generated models based on aerial photographs, spatial data obtained by use of a terrestrial laser scanner, and data obtained using an aerial laser scanner. An additional aim was to use an ANN to evaluate the influence of individual components of the LS parameter on its final distribution. We focused on influence of resolution of DEM, generated based on geospatial data from TLS, ALS and AP, on the LS parameter components. We generated the particular components of the mode (S, L, slope, flow direction, flow accumulation, m and β) and we carried out analysis of absolute influence of its parts on the LS value.

Materials and Methods
The rapid development and accessibility of DEM data promotes the use of techniques that include picture transformation and software [57] in developing new approaches to assessing and improving the precision of the LS parameter. For modelling and assessment of the LS parameter, we used the ArcGIS program, version 10.5, with a set of modules included for data analyses. In the analyses, we choose a 40 m × 200 m field of ground under continuous tillage. The investigations were carried out based on three independent methods for the extraction and processing of geospatial data: TLS (Terrestrial Laser Scanning), ALS (Aerial Laser Scanning), and AP (Aerial Photographs).
The study was performed on a farm in the municipality of Szonowice  (Figure 1). The areas chosen eliminated artefacts from ground uncertainties connected with the use of agricultural practices-e.g., technical roads, balks, field roads, grass edges-to avoid an unaccounted-for influence on our investigations. Such influences have been shown to be significant, such as in [27], where local features of the landscape caused a proportional reduction in the LS coefficient.

Source of DEMs
As noted above, three methods were used to obtain geospatial data: Point clouds originating from various sources have various degrees of accuracy, but they have a common georeference, which makes it possible to carry out comparative analyses of the efficacy of their use. As reference data, the point cloud generated from Terrestrial Laser Scanning (TLS) had the highest precision and the greatest point cloud density.

Source of DEMs
As noted above, three methods were used to obtain geospatial data: Point clouds originating from various sources have various degrees of accuracy, but they have a common georeference, which makes it possible to carry out comparative analyses of the efficacy of their use. As reference data, the point cloud generated from Terrestrial Laser Scanning (TLS) had the highest precision and the greatest point cloud density.

Source of DEMs
As noted above, three methods were used to obtain geospatial data: Point clouds originating from various sources have various degrees of accuracy, but they have a common georeference, which makes it possible to carry out comparative analyses of the efficacy of their use. As reference data, the point cloud generated from Terrestrial Laser Scanning (TLS) had the highest precision and the greatest point cloud density.

Aerial Photographs (APs)
The point cloud from aerial photographs was generated using the SfM (Structure from Motion) method. Stereoscopic photographs of a given area were obtained from the Store of Institute of Geodesy and Cartography, Warsaw, Poland [58]. To orient the project and fit it into the National System of Geodetic Coordinates, PUWG 1992 (EPSG-2180), coordinates of photo points located in the area during photogrammetric flight were used. The photographs were assembled through the aerotriangulation process available in the Agisoft PhotoScan Professional program with 0.15 m accuracy. The margin of error, the size of field pixels, and remaining errors resulting from the fact that the point cloud was obtained using derivative information provided by photographs, should be added to the results as well.

Terrestrial Laser Scans (TLSs)
Measurements were carried out with a terrestrial laser scanner (TLS) Leica P40 ScanStation. Accuracy of the device is stated as: range accuracy 1.2 mm + 10 ppm over its full range, with an angular accuracy of 8" horizontal, 8" vertical, which enables it to obtain results for 6 mm to 100 m heights. The device calibration process is carried out under laboratory conditions. In the field, measurements were carried out at a scanning resolution of 2 mm/10 m-i.e., at a distance of 10 m from the scanner, points are recorded horizontally and vertically every 2 mm. For georeferenced registration, uniformly distributed target points on the investigated field were used, and then the GNSS-RTK measurement was carried out. The RTK measurement was carried out with an accuracy of 0.03 m. The coordinate system used was the National Coordinate System, PUWG 1992 (EPSG-2180). For location of the point cloud within the global coordinate system, the Leica Cyclone program was used. As the accuracy of the fitting was below 0.01 m, the accuracy of every point in the point cloud set was 0.03 m.

DEMs Generation
To compare the three methods used to gather data on the investigated sites, a regular 1 × 1 m network was generated. The final comparison was carried out based on 6943 points. In the analyses, aberrant points attributed to resolution differences between generated clouds were rejected ( Figure 3).

Aerial Photographs (APs)
The point cloud from aerial photographs was generated using the SfM (Structure from Motion) method. Stereoscopic photographs of a given area were obtained from the Store of Institute of Geodesy and Cartography, Warsaw, Poland [58]. To orient the project and fit it into the National System of Geodetic Coordinates, PUWG 1992 (EPSG-2180), coordinates of photo points located in the area during photogrammetric flight were used. The photographs were assembled through the aerotriangulation process available in the Agisoft PhotoScan Professional program with 0.15 m accuracy. The margin of error, the size of field pixels, and remaining errors resulting from the fact that the point cloud was obtained using derivative information provided by photographs, should be added to the results as well.

Aerial Laser Scans (ALSs)
A point cloud was generated from aerial laser scans obtained from the Store of Institute of Geodesy and Cartography, Warsaw, Poland

Terrestrial Laser Scans (TLSs)
Measurements were carried out with a terrestrial laser scanner (TLS) Leica P40 ScanStation. Accuracy of the device is stated as: range accuracy 1.2 mm + 10 ppm over its full range, with an angular accuracy of 8" horizontal, 8" vertical, which enables it to obtain results for 6 mm to 100 m heights. The device calibration process is carried out under laboratory conditions. In the field, measurements were carried out at a scanning resolution of 2 mm/10 m-i.e., at a distance of 10 m from the scanner, points are recorded horizontally and vertically every 2 mm. For georeferenced registration, uniformly distributed target points on the investigated field were used, and then the GNSS-RTK measurement was carried out. The RTK measurement was carried out with an accuracy of 0.03 m. The coordinate system used was the National Coordinate System, PUWG 1992 (EPSG-2180). For location of the point cloud within the global coordinate system, the Leica Cyclone program was used. As the accuracy of the fitting was below 0.01 m, the accuracy of every point in the point cloud set was 0.03 m.

DEMs Generation
To compare the three methods used to gather data on the investigated sites, a regular 1 × 1 m network was generated. The final comparison was carried out based on 6943 points. In the analyses, aberrant points attributed to resolution differences between generated clouds were rejected ( Figure 3).

LSUSLE Calculations
The L and S layers of the original equation, and LS layer later identified as the "coefficient of topography," were calculated using a surface formulation [59] composed of unequivocally-defined coefficients produced by the selected algorithms, resulting in the digital generation of morphological features determined by hydrological slope parameters. To evaluate L, we used the ArcToolbox-

LS USLE Calculations
The L and S layers of the original equation, and LS layer later identified as the "coefficient of topography," were calculated using a surface formulation [59] composed of unequivocally-defined coefficients produced by the selected algorithms, resulting in the digital generation of morphological features determined by hydrological slope parameters. To evaluate L, we used the ArcToolbox-Spatial Analyst Tools-Hydrology modules [60] and the algorithm of Desment and Govers [31]. The raster network cell dimension is very important for any evaluation of S [61] because slope increases as the dimension of the raster cell decreases. After determining L and S, LS at various resolutions of network cells was calculated as well.
To calculate L, we used equation proposal by Desmet and Govers [31]: where: D is the raster resolution; A(i,j) is the unit area of feeding at the entrance to cell (i,j); m is the exponent of slope length indicator, where: β is the ratio of rill to inter-rill erosion for conditions when the soil is moderately susceptible to both rill and inter-rill erosion and θ is the slope angle; x is the coefficient correcting the length of flow path through raster cell, depending on flow direction and calculated based on exposure; L is the slope length coefficient (L) factor; 22.13 is the USLE unit plot length; and λ is the cumulative slope length in meters [26]. The first step in the generation of L was the evaluation of slope linear length. To this end, based on the DEMs, the layers of flow direction and flow accumulation were determined. Initially, we generated a network of direction using the one-way punctual model of flow, D8 (Figure 4). Calculations were carried out according to the equation proposed by Wilson and Gallant [62]: where: Z is the number of adjacent cells; h is the resolution of the GRID model; h∅ (i) is distance between cell centers-1 for those located in cardinal directions (N, E, S, W)-and square root of 2 for the remaining directions.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 21 Spatial Analyst Tools-Hydrology modules [60] and the algorithm of Desment and Govers [31]. The raster network cell dimension is very important for any evaluation of S [61] because slope increases as the dimension of the raster cell decreases. After determining L and S, LS at various resolutions of network cells was calculated as well.
To calculate L, we used equation proposal by Desmet and Govers [31]: where: D is the raster resolution; A(i,j) is the unit area of feeding at the entrance to cell (i,j); m is the exponent of slope length indicator, where: β is the ratio of rill to inter-rill erosion for conditions when the soil is moderately susceptible to both rill and inter-rill erosion = (sin /0.0896) 3.0(sin ) . + 0.56 (3) and θ is the slope angle; x is the coefficient correcting the length of flow path through raster cell, depending on flow direction and calculated based on exposure; L is the slope length coefficient (L) factor; 22.13 is the USLE unit plot length; and λ is the cumulative slope length in meters [26]. The first step in the generation of L was the evaluation of slope linear length. To this end, based on the DEMs, the layers of flow direction and flow accumulation were determined. Initially, we generated a network of direction using the one-way punctual model of flow, D8 (Figure 4).
where: Z is the number of adjacent cells; h is the resolution of the GRID model; h (i) is distance between cell centers-1 for those located in cardinal directions (N, E, S, W)-and square root of 2 for the remaining directions.  The network of flow directions enabled us to determine how water flows within the area and creates a network of flow accumulation in which every following cell on the flow path accumulates water from a higher part of the area. As a consequence, a map of flow accumulation could be produced for each method ( Figure 5).
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 21 The network of flow directions enabled us to determine how water flows within the area and creates a network of flow accumulation in which every following cell on the flow path accumulates water from a higher part of the area. As a consequence, a map of flow accumulation could be produced for each method ( Figure 5). For the evaluation of the S, the McCool et al. [24,25] method was used. They found that soil loss occurs faster in slopes that were steeper than 9%. Renard et al. [18] adopted this algorithm in RUSLE for the S-factor estimation based on the slope gradient: S = 10.8 × sin θ + 0.03, where: slope gradient < 0.09 S = 16.8 × sin θ − 0.5, where: slope gradient ≥ 0.09 where: θ is the gradient of slope in degrees.
Finally, we calculated the LS parameter as the product of the parameters L and S [63].

Statistical Analysis
To evaluate whether measurement data (as demonstrated by the LS parameter) had a normal distribution, we used the Kołmogorov-Smirnov test, and Statistica software, release 12. Homogeneity of measurement sequences for investigated features was checked using the non-parametric Kruskal-Wallis test [64]: where: and Ri is the total rank of the variables in a combined sample; ni is a random sample size; n is the sequence size; and k is the number of samples. Comparison of L and S evaluated the techniques for taking images to generate DEMs (Methods I-III) was carried out using the linear correlation and error measures [65] listed below: Mean error of prediction (MEP): Root mean square error (RMSE): Mean percentage error (MPE): For the evaluation of the S, the McCool et al. [24,25] method was used. They found that soil loss occurs faster in slopes that were steeper than 9%. Renard et al. [18] adopted this algorithm in RUSLE for the S-factor estimation based on the slope gradient: S = 10.8 × sin θ + 0.03, where: slope gradient < 0.09 S = 16.8 × sin θ − 0.5, where: slope gradient ≥ 0.09 where: θ is the gradient of slope in degrees. Finally, we calculated the LS parameter as the product of the parameters L and S [63].

Statistical Analysis
To evaluate whether measurement data (as demonstrated by the LS parameter) had a normal distribution, we used the Kołmogorov-Smirnov test, and Statistica software, release 12. Homogeneity of measurement sequences for investigated features was checked using the non-parametric Kruskal-Wallis test [64]: where: and R i is the total rank of the variables in a combined sample; n i is a random sample size; n is the sequence size; and k is the number of samples. Comparison of L and S evaluated the techniques for taking images to generate DEMs (Methods I-III) was carried out using the linear correlation and error measures [65] listed below: Mean error of prediction (MEP): Root mean square error (RMSE): Mean percentage error (MPE): Model efficiency (ME) [66,67]: where: C TLS i is the values obtained by the TLS; C ALS/AP i are the ones obtained by ALS or AP; and n is the number of data points; C is the mean value. The relative and absolute influence of the individual elements generating the LS coefficient was evaluated by use of the ANNs of the multi-layer perceptron (MLP) type using the Statistica software, release 13.1. A total of 70% of all variables were used for the learning process; 15% were used for validation; and 15% were used to test the model. A quasi-Newton algorithm with a Broyden-Fletcher-Goldfarb-Shanno (BFGS) modification was selected for the learning of neural network. The sum of squares (SOS) was treated as the error function. Relationships between the LS parameter and model parts as independent variables were applied in all data sets. The MLP consists of three layers of neurons: input, hidden and output (it is written as m-n-p where m, n, p are numbers of perceptrons in the following layers). Each neuron had a number of inputs-from outside the network or the previous layer-and a number of outputs-leading to the subsequent layer or out of the network. The computed output response of a is based on the weighted sum of all its inputs according to an activation function [68].

Results and Discussion
Calculations of the LS USLE parameter were carried out by using the original equation proposed by Desmet and Govers [31] and implemented using an automatic geological analysis (generated by ArcGIS) system, which contributed to the precision evaluation of flow accumulation. The data set of the LS coefficient was calculated using DEMs obtained from the three methods. This connected approach combining GIS software tools with DEMs of high-resolution data was used with success on regional and European scales by Panagos et al. [27].
Spatial distribution of the LS parameter and its components was determined based on DEMs generated by the ARCGIS program as presented in Figure 6.
The parameter S was between 0.06 and 4.17; the parameter β was between 0.06 and 1.89; parameter m was between 0.06 and 0.65; the parameter L was between 1.01 and 3.84; flow direction was between 1 and 128; flow accumulation was between 0 and 2868, slope was between 0.19 • and 16.13 • ; LS was between 0.06 and 6.32 ( Table 2). The lowest and at the same time the highest values were generated by means of Method I, Aerial photographs. The highest values of L were seen in Method II, Aerial laser scans. Panagos et al. [27] obtained parameters for Poland that include a mean of 0.52 and a standard deviation of 0.86; comparable values to those generated within this study, despite the significantly different DEM resolutions. The parameter S was between 0.06 and 4.17; the parameter β was between 0.06 and 1.89; parameter m was between 0.06 and 0.65; the parameter L was between 1.01 and 3.84; flow direction was between 1 and 128; flow accumulation was between 0 and 2868, slope was between 0.19° and 16.13°; LS was between 0.06 and 6.32 ( Table 2). The lowest and at the same time the highest values were generated by means of Method I, Aerial photographs. The highest values of L were seen in Method II, Aerial laser scans. Panagos et al. [27] obtained parameters for Poland that include a mean of 0.52 and a standard deviation of 0.86; comparable values to those generated within this study, despite the significantly different DEM resolutions. Comparison of the LS USLE models generated based on the geospatial data from the APs, ALSs, and TLSs are presented in Figure 7 and Table 3. The greatest accordance both within the DEMs and the LS parameter appeared between the ALS and TLS methods, determination coefficient R 2 was 0.9992 ( Figure 7C) and 0.8185 ( Figure 7F) respectively. The least accordance characterized the methods AP and TLS in a case of the parameter LS USLE R 2 = 0.501 ( Figure 7D).  ANNs were designed to validate the relative and absolute influence of the chosen parameters when evaluating the LS parameter. To this end, the MLP network type was chosen. An analysis of the quality of MLP networks should be uniformly obtained from the data sets of the LS parameter. The quasi-Newton algorithm can easily be optimized to forecast the relationship between the chosen geospatial parameters and the LSUSLE parameter. For the AP method, the activation functions were the logistic and tanh for hidden and output perceptrons, respectively; for the ALS method, the activation functions were logistic and linear, respectively; and for the TLS methods, the activation functions were tanh and exponential, respectively. Quality of learning for all three methods was extra satisfactory at 1 or near 1. Similarly, quality of testing and validating were near 1 or 1. The ANN's   ANNs were designed to validate the relative and absolute influence of the chosen parameters when evaluating the LS parameter. To this end, the MLP network type was chosen. An analysis of the quality of MLP networks should be uniformly obtained from the data sets of the LS parameter. The quasi-Newton algorithm can easily be optimized to forecast the relationship between the chosen geospatial parameters and the LS USLE parameter. For the AP method, the activation functions were the logistic and tanh for hidden and output perceptrons, respectively; for the ALS method, the activation functions were logistic and linear, respectively; and for the TLS methods, the activation functions were tanh and exponential, respectively. Quality of learning for all three methods was extra satisfactory at 1 or near 1. Similarly, quality of testing and validating were near 1 or 1. The ANN's learning sample error for the examined methods were 1.8 × 10 −5 , 5 × 10 −6 , and 0, respectively; for the testing sample, it was 2.4 × 10 −5 , 4 × 10 −6 , and 0, respectively; and for the validation sample, it was 1.8 × 10 −5 , 6 × 10 −6 , and 0, respectively (Table 4a,b). Analysis of the sensitivity of ANN modelling by MLP showed that geospatial parameters S and Slope had the greatest absolute influence (50.7% and 28.4%, respectively) for the AP method; parameters S and L were the greatest influence (67.2% and 11.3%, respectively) for the ALS method; and parameters S and Slope were the greatest influence (45.6% and 42%, respectively) for the TLS method (Table 5).  The LS parameter was calculated with great accuracy using high resolution DEMs. As the visual analysis of evaluated parameters ( Figure 5) shows, DEM resolution influences the spatial distribution of both parameter LS and its components. High resolution DEMs reveal geomorphological changes in the ground with greater precision, and as a consequence, the evaluated values have greater accuracy. While Panagos et al. [27] carried out their analyses of erosion and particular parameters on the model of Europe, using DEMs with a 25-m resolution (earlier, 100-m resolution), the investigations we carried out show this is all possible on a smaller scale.

LS-Factors
S, which measures the effect of slope steepness, had the greatest influence, while L influences the length of the ground slope, and in our case, L turned out to be the least essential factor in modelling the LS parameter. This can be explained by the fact that the measuring field was regular.
Flow accumulation is the total amount of water that inflows into every ground cell [69]. The traditional method of wobbulation requires repetition for the purpose of obtaining flow accumulation, what was introduced by O'Callaghan and Mark [39], as well as Tarboton [70], and is still used because of its simplicity and ease of use. However, this method is time consuming, especially in a case of vast areas, requiring substantial transformation of data [71][72][73][74][75]. Yao et al. [76] discovered an efficient method of DEM creation by sweeping various places in various ways. Su et al. [77] used the drainage-basin tree algorithm for determining how flow accumulation efficiency could be improved. Bai et al. [78] used a sequence of sorted pixels (by elevation) to collect iterations of an area of an elevation from top to bottom. Collection of flow information based on DEMs depends on flow direction, while most of the other models solve flow accumulation without regard to flow direction.
DEMs often represent of ground surface topography as a Cartesian network, triangle irregular network (TIN), or contour flow networks [79]. Regular network DEMs are convenient for calculations and conversions; the simple structure allows it to be widely used for topographic analysis and visualizations of the ground in hydrological applications [80]. Because DEMs approximate ground characteristics, the accuracy of topography, which is connected with its hydrological applications, depends on such DEM features as slope analysis, the density of river networks, flow paths, and the topographic index. DEMs used in hydrological models must be first and foremost hydrologically correct. Most DEMs contain topographic depressions, defined as areas without outlet and often called sinks and concavities [81]. In regular network DEMs, relief lowerings are areas having one or more adjacent cells which are lower than all the surrounding cells [82][83][84].
There is a strong correlation between the generalization of height data and the level of accuracy of the calculated mass of eroded material [14]. Most DEMs currently produced contain a great number of depressions that represent incorrect, false topographic errors [81,85]. The L particularly depends on the accuracy of DEMs and the precision of various calculation methods [86]. The accuracy of DEMs is also quite susceptible to horizontal resolution, vertical precision, and the density of sampling points. With the decrease of horizontal resolution, slope gradient decreases [87], and the total length of flow and density of drainage also have a decreasing tendency [88] as the concrete area of the basin increases [71]. The density of DEM sampling points on an increasing slope causes a decrease in the basin area [89][90][91]. Horizontal resolution and the density of sample points influences the total length of flow, and the slope gradient and alimentation area influence the L parameter.
The precision of transformation of the L coefficient is very sensitive to the chosen calculation methods and basic attributes that originate from DEMs, as presented above. The raster digital DEMs representing continuity of ground height above the basic level are commonly used in automatic hydrological analyses and in the generation of various basin features.
Based on the results of the LS parameter (Models I, II, and III), the A parameter (long-term average annual soil loss of the USLE model) values were calculated The particular parameters of the model (R = rainfall erosivity factor, K = soil erodibility factor, C = cropping management factors, P = conservation practices factor) were determined based on the investigations carried out for Europe by Panagos et al.  Depending of exactness of input geospatial data, DEM from TLS and from ALS are comparable, what causes comparable results of A. This is result of character of cloud points from laser scanning, that keep geometry of objects. In a case when AP was used and point cloud was generated based on the pictures, the highest values of the parameter were obtained, what is result of use of cloud, being derivative of pictures and many treatments and processes of input data adaptation. Use of various types of input data can cause changes even up to 20%.
The obtained results were worked out using the values of R, K, C, and P and incorporating the existing conditions and variability of the LS parameter and the results of various exactness of source materials analyses for the geospatial data.
The obtained results can be extrapolated for slopes and total catchments, but the extensions are subject to limitations and problems, such as: technical and logistic restrictions on equipment, cost of equipment, and accessibility of devices and applications. Attention should be paid to the scale of the area being examined to choose the optimum measurement techniques. Exactness and precision of measurement are generally hampered as the area of investigation grows. Model precision is a derivative of the technology used and the nature of the geospatial data collected; exactness directly influences the quality of the DEM models generated.
Along with the source of data coming from the direct field measurements using AP, ALS, and TLS, additional data sources may be available like: SRTM and ASTER data and other satellite images (i.e., materials from Landstat, Modis, Sentinel 1, and Sentinel 2) for construction of digital elevation models. Taking into account the lower resolution, such data may useful in large-scale global applications, for example by regions or country or continent, but it is not recommended as a substitute for local scale (field data).

Conclusions
We presented DEMs and LS under the USLE model generated using three popular methods of data collection used in environmental investigations, AP, ALS, and TLS. The results of parameters generated by each method were compared. We confirmed that the highest influence shaping the LS parameter was the S parameter (index of ground slope). Additionally, we concluded that in investigations over vast areas, terrestrial laser scanners are not optimal due to the additional labor and time consumption they require. We further established that measurements sometimes cannot be carried out via TLS because ground obstacles the limit the practical extent of laser work.
In our opinion, aerial photographs or clouds of points from an aerial laser scanner are most suitable in such situations. The data generated from these methods turned out to be highly precise (comparable with the terrestrial laser scanner data). However, these methods were compared over a small area, without plant cover and without ground obstacles such as balks, reducing the uncertainty of measurement for the scope of this project. To assess the accuracy of the other two methods, the DEM and the LS parameter generated by means of the terrestrial laser scanner were taken as base or reference values because of their potentially higher accuracy.
We investigated the influence of various parameters on the LS coefficient having topological connections (connected with relief) and compared results from the terrestrial laser scanner with results from the two other methods. In spite of the fact that terrestrial laser scanner had the lowest error measure values of (that is, it was the most accurate method), we think this method has faults-such as costs and time of obtaining and transforming and the time needed for scanning of ground-when compared with the other methods. TLS may be most successful in micro applications evaluating various environmental parameters; for macro purposes, aerial laser scanners and aerial photographs are considerably quicker and cheaper, and at the same time have low error measures. Of these two methods, we believe the aerial laser scanner is most preferred.