Assessment of Soil Degradation by Erosion Based on Analysis of Soil Properties Using Aerial Hyperspectral Images and Ancillary Data, Czech Republic

The assessment of the soil redistribution and real long-term soil degradation due to erosion on agriculture land is still insufficient in spite of being essential for soil conservation policy. Imaging spectroscopy has been recognized as a suitable tool for soil erosion assessment in recent years. In our study, we bring an approach for assessment of soil degradation by erosion by means of determining soil erosion classes representing soils differently influenced by erosion impact. The adopted methods include extensive field sampling, laboratory analysis, predictive modelling of selected soil surface properties using aerial hyperspectral data and the digital elevation model and fuzzy classification. Different multivariate regression techniques (Partial Least Square, Support Vector Machine, Random forest and Artificial neural network) were applied in the predictive modelling of soil properties. The properties with satisfying performance (R 2 > 0.5) were used as input data in erosion classes determination by fuzzy C-means classification method. The study was performed at four study sites about 1 km 2 large representing the most extensive soil units of the agricultural land in the Czech Republic (Chernozems and Luvisols on loess and Cambisols and Stagnosols on crystalline rocks). The influence of site-specific conditions on prediction of soil properties and classification of erosion classes was assessed. The prediction accuracy (R 2) of the best performing models predicting the soil properties varies in range 0.8–0.91 for soil organic carbon content, 0.21–0.67 for sand content, 0.4–0.92 for silt content, 0.38–0.89 for clay content, 0.73–089 for Fe ox , 0.59–0.78 for F ed and 0.82 for CaCO 3. The performance and suitability of different properties for erosion classes' classification are highly variable at the study sites. Soil organic carbon was the most frequently used as the erosion classes' predictor, while the textural classes showed lower applicability. The presented approach was successfully applied in Chernozem and Luvisol loess regions where the erosion classes were assessed with a good overall accuracy (82% and 67%, respectively). The model performance in two Cambisol/Stagnosol regions was rather poor (51%–52%). The results showed that the presented method can be directly and with a good performance applied in pedologically and geologically homogeneous areas. The sites with heterogeneous structure of the soil cover and parent material will require more precise local-fitted models and use of further auxiliary information such as terrain or geological data. The future application of presented approach at a regional scale promises to produce valuable data on actual soil …


Introduction
Soil erosion is one of the most significant and widespread forms of soil degradation in Europe [1].In the Czech Republic, soil erosion of different types (water, wind, and tillage) is identified as the most severe type of soil degradation [2], which significantly affects soil functions, production of agricultural crops and the quality of water resources.The degradation risk by erosion is currently relatively well known at both local and regional scale due to evaluation using erosion models [3].According to these calculations, nearly 50% of arable land in the Czech Republic is endangered by water erosion and about 10% by wind erosion [4].Nevertheless, the widely applied erosion models (mostly based on USLE/RUSLE models) only provide the estimation of potential erosion.The calculation of the actual soil losses due to soil erosion and estimation of real long-term soil degradation are still insufficiently developed in spite of being essential for policy and management purposes [5,6].
Currently, the assessment of soil degradation by erosion is solved mainly at a local scale using different methods, e.g., direct field monitoring, observations and measurements obtained from distributed point datasets (identification of erosion features and soil profile truncation, assessment of Caesium-137, risk elements and other relevant soil properties) [7][8][9][10][11].At a regional or global scale, information on the level of soil degradation by erosion is often based on expert knowledge approach consisting in the extrapolation of acquired local data using ancillary data.This is due to the lack of reliable methods, since detailed soil mapping is not financially feasible and it is extremely difficult particularly at this scale [12].Nevertheless, estimates based on the physical evidence are much required [5,13].
Application of remote sensing represents a possible solution in the assessment of soil degradation by erosion at a regional and global scale.The main advantage of remote sensing methods consists in the possibility of rapid data acquirement from large areas in a detailed spatial resolution.However, its potential has not been fulfilled and there are still many theoretical and technical gaps and unanswered questions [13,14].Continual progress in measurement techniques, image processing algorithms, development of computation models, new sensors and new satellite missions have brought new data and methods.Use of hyperspectral data in the erosion assessment represents a promising method, particularly in the context of the expectations forthcoming spaceborne hyperspectral sensors with high signal-to-noise ratio (SNR) and pixel size from one to several tens of metre, such as German Environmental Mapping and Analysis Program (EnMAP), Italian Hyperspectral Precursor of the Application Mission (PRISMA), NASA's Hyperspectral Infra-Red Imager (HyspIRI), Japanese Hyperspectral Imager Suite (HISUI), French HypXIM, israel-italian Spaceborne Hyperspectral Applicative Land and Ocean Mission (SHALOM) [15] or Chinese TianGong-1 [16].
Land degradation by soil erosion can be assessed directly using different remote sensing methods [14,17,18]: (i) spatial delineation of degraded land and erosion features (rills, gullies, and sediment depositions) using visual interpretation of aerial images [19], computer pre-processing [20][21][22] or automatic classification methods [23][24][25]; (ii) soil loss measurement methods based on digital elevation models (DEM) comparison (laser altimetry or radar interferometry); and (iii) assessment of soil degradation stages and monitoring of soils affected by soil erosion and their properties.These methods are based on classification techniques or use different mathematical and statistical procedures to assess the correlation between erosion signs and their spectral reflectance [15,26,27].Soil erosion and accumulation affect chemical and physical properties of upper layer of soils; thus, spectral reflectance is changed and eroded soils have different spectral response from non-eroded "healthy" soils [28][29][30][31][32][33][34].
Soil properties with direct relation to the spectral signature and with the best prediction performance on one hand and a high variability due to erosion-accumulation processes on the other hand include soil organic carbon (SOC) concentration, particle size distribution (particularly clay content), CaCO 3 concentration, Fe and Al oxides ratio and water content [15,[35][36][37][38][39][40].Other soil properties used for determination of erosion-affected areas are soil colour, pH, structure or coarse fragments content [41][42][43][44][45].
Imaging spectroscopy or generally analysis of soil spectral characteristics has a great potential for soil erosion assessment [36,46], as the hyperspectral data, thanks to their high spectral resolution, allows for more efficient quantification of spectral features in comparison with, e.g., multispectral data.Thus far, only few studies have employed this method in soil erosion assessment [47][48][49], a majority of them in semi-arid conditions [30,33,45,[50][51][52].Several studies dealt with the monitoring of specific soil properties related to soil erosion, such as development of soil crust [50] or changes in soil moisture [47,48].Hyperspectral data have scarcely been used for distinguishing soils influenced by soil erosion.Hill et al. [51,52] used linear spectral mixture analysis and k-means clustering for classification of erosion stages based on skeleton content determination.Schmid et al. [33] distinguished erosion classes using support vector machine classification method and development of image-derived endmembers for each erosion class.Lin and Zhou [31,45] have analysed a spectral response of different eroded soils in subtropical China.They have compared a direct identification of erosion stages by analysing laboratory hyperspectral data and an indirect identification using prediction of soil properties; the latter performed better due to mixing and crossover effect among erosion groups.
However, there are still many limitations in the use of the image hyperspectral data for soil erosion assessment.On one hand, there are restraints connected to data collecting, such as atmospheric attenuation, signal-to-noise ratio, image resolution and, above all, the Bidirectional Reflectance Distribution Function (BRDF) effect [36].On the other hand, the limitations are often caused by surface characteristics making the direct soil assessment impossible or difficult (soil covered with vegetation, litter, dust or soil crust).Moreover, the spectral information is considerably influenced by soil properties (soil moisture, soil roughness, soil texture or size of soil aggregates) and heterogeneity of environmental settings.This is valid mainly in an erosion relief, where the combination of soil loss and redeposition may lead to similar soil surface properties in eroded and accumulated soils [52,53].A successful application of hyperspectral data then requires a very thorough planning of flight campaigns, precise image corrections and a sufficient amount of reference soil data.In temperate climate, a multi-temporal approach must be considered when assessing larger areas, as in each period some of the soils are covered with vegetation.
In spite of the progress in knowledge and algorithms (especially in terms of atmospheric correction), more studies are needed for refinement of this method and widening its application.
The aim of this study is to analyse the possibilities of imaging spectroscopy using aerial hyperspectral data for the assessment of soil degradation by erosion on arable land in temperate agriculture regions with different dominant soil unit.The specific objectives of this work were to: (1) quantify the prediction performance of selected soil surface properties derived from hyperspectral data in context of various soil environment, spectral pre-processing methods and regression methods; and (2) determine erosion classes that represent soils differently impacted by soil redistribution by applying the predicted soil surface properties as site-specific indicators and input data in classification models of erosion classes distinguishing.

Regional Settings
Four study sites representing the most extensive soil units of agricultural land in the Czech Republic were chosen.Soil regions are represented by Chernozems and Luvisols on loess and Cambisols and Stagnosols on crystalline and sedimentary rocks.Locations of the study sites are depicted in Figure 1.Sites are similar in terms of terrain characteristics (dissected relief including a set of following terrain units: side valley, toe-slope, plateau and back-slope), land management (long-term tillage, no conservation practices, plough depth 25 cm, 5-6 course rotation based on the Norfolk system) and climatic conditions (rain erosivity).All study sites are located within regions with high rates of both potential and actual degradation by erosion according to Maps of long-term average annual soil loss by water erosion [3] and database of Monitoring of soil erosion on agricultural land [8].The types of erosion influencing the soil redistribution include mainly sheet and rill erosion and tillage erosion.Wind erosion is presumed to have a minor impact at the Chernozem site (Šardice).Environmental settings and details of each locality are given in Table 1.
and rill erosion and tillage erosion.Wind erosion is presumed to have a minor impact at the Chernozem site (Šardice).Environmental settings and details of each locality are given in Table 1.

Flight Campaign (Imaging Spectroscopy) and Image Pre-Processing
Four image scenes were acquired during four flight campaigns from the VNIR (visible and near-infrared) sensor-Compact Airborne Spectrographic Imager (CASI1500; 370-1040 nm) and SWIR (short-wavelength infrared) sensor-Shortwave infrared Airborne Spectrographic Imager (SASI600; 960-2440 nm).Both sensors are developed by Itres Ltd. (Calgary, AB, Canada).The sensors are pushbroom sensors with a Field Of View of 40°.Sensors acquired data in 72 spectral bands in the VNIR with a full width at half maximum (FWHM) of 10 nm (CASI) and 100 bands in the SWIR with FWHM ≈ 15 nm (SASI).Sensors were mounted on Cessna 208B Grand Caravan aircraft.The flight campaigns took place during dry conditions (minimum of five days after last rain), and

Flight Campaign (Imaging Spectroscopy) and Image Pre-Processing
Four image scenes were acquired during four flight campaigns from the VNIR (visible and near-infrared) sensor-Compact Airborne Spectrographic Imager (CASI1500; 370-1040 nm) and SWIR (short-wavelength infrared) sensor-Shortwave infrared Airborne Spectrographic Imager (SASI600; 960-2440 nm).Both sensors are developed by Itres Ltd. (Calgary, AB, Canada).The sensors are pushbroom sensors with a Field Of View of 40 • .Sensors acquired data in 72 spectral bands in the VNIR with a full width at half maximum (FWHM) of 10 nm (CASI) and 100 bands in the SWIR with FWHM ≈ 15 nm (SASI).Sensors were mounted on Cessna 208B Grand Caravan aircraft.The flight campaigns took place during dry conditions (minimum of five days after last rain), and study fields were prepared for seeding-ploughed and harrowed without vegetation or litter.Details of individual flight campaigns are given in Table 2. Data acquisition and pre-processing (geometrical and atmospheric corrections) were realized by the Flying Laboratory of Imaging Systems (FLIS) [54,55] operated by Global Change Research Institute CAS (Brno, Czech Republic).The radiometric correction was performed using the RadCorr Ver.9.2.6.0 (Itres Ltd.) and by means of calibration parameters obtained in laboratory.Data output was given in radiometric units (µW•cm −2 •sr −1 •nm −1 ).Atmospheric correction was applied to remove the effect of atmospheric influences and convert radiance values into at-surface reflectance.MODTRAN radiative transfer model incorporated into the program ATCOR-4 ver.7.0 was used for this purpose.Algorithm BREFCOR was used for BRDF effect reduction.The geometric correction and geo-referencing of images was performed using the Geocor tool (Itres Ltd.), based on data recorded by the on-board GPS/IMU sensors and digital elevation model.Image data were resampled by nearest neighbour method and transformed into the UTM map projection.
Image data were masked by bare soil mask (by means of knowledge of study sites and NDVI computing) and resampled to 6-metre spatial resolution (due to random noise reduction in data and corresponding to the size of sampling plots) before processing.Spectral bands highly influenced by the absorption in the atmosphere and bands on the edge of spectra influenced by noise were removed from the dataset.The data finally include 49 bands for CASI (400-750 and 770-880 nm) and 53 bands for SASI (1000-1080, 1200-1300, 1490-1770, 2060-2370 nm).Spectral data extraction from data cubes, as well as other data manipulation, was done using the R software [56].

Data Collection and Soil Analysis
Soil sampling was carried out during four field campaigns following each flight campaign.Fifty samples were taken from each site in an optimized network of borings fashioned using cLHS (conditioned Latin hypercube sampling [57]) stratified random strategy.Terrain attributes and spectral data obtained within hyperspectral campaign were used as feature space variables.This approach ensures the cover of maximal variation of each variable.Exact position of sampling sites was measured by GPS Trimble GeoXM receiver, with a post-processing accuracy of approximately 1m.Soil unit, profile stratigraphy, soil depth, and thickness of horizons were determined by description of gouge auger core.The composite soil samples for analysis of SOC, texture classes, CaCO 3 and Fe oxides content were taken from each site at 0-10 cm depth.

Statistical Analysis of Soil Properties
Different multivariate regression methods were used to establish relationships between image spectral data and soil properties.Different pre-treatment and regression methods were tested.Selected soil properties data were used in regression models as response variables, and spectral reflectance data (transformed spectra) as predictors.Terrain derivatives were used for A horizon thickness establishing.The best resulting regression models were then applied on all the image data with the aim of predict spatial variability of particular soil surface properties and create maps.All spectroscopic and statistical analyses were implemented in R environment [56].

Pre-Processing
Before performing statistical analyses, mathematical pre-treatments were applied to adjust raw reflectance spectra.These transformations were performed to improve the prediction accuracy by noise reduction and mitigation of influence of environment and sensing properties.Continuum removal (CR), Savitzky-Golay filter (third order polynomial smoothing and 5 bands window widths) with first (SG 1st) and second derivatives (SG 2nd) [60] and standard normal variate transformation (SNV) were applied.Reflectance was transformed into absorbance (log (1/R) as well.Transformation of CASI and SASI spectra were carried out separately as spectra was not acquired with same bandwidth.Bands on the edge of spectra affected by pre-treatment due to moving window averaging were removed from final dataset.R package Prospectr [61] was used for spectra pretreatment.

Calibration and Validation
The dataset was divided into training set used for fitting models and validation set used for assessing the prediction accuracy of each model.Dividing in a 3:1 ratio was performed by random stratified sampling.The procedure of prediction modelling was carried out using caret package [62] of the R software.In the first step, the training set was used for fitting the model.Model performance was assessed through 5-fold cross-validation of the training set.Random search method generated by caret package, providing the automatic search of parameter values, was used for selecting the best parameters of models.Single model with yielding smallest root mean squared error of cross validation (RMSE CV ) value was selected for subsequent validation on the validation set.The final prediction accuracy was assessed with root mean squared error of prediction (RMSE P ) and coefficient of determination (R 2 ).The ratio of standard deviation to standard error of prediction (RPD) was also calculated alternatively to R 2 .Goodness of fit was visually inspected through a plot depicting the observed values against the predicted values.

Assessment of Soil Erosion Classes
Assessment of the different erosion classes on the study sites was performed in several steps: (1) definition of groups of site-specific erosion classes based on soil data from field campaigns; (2) assessment of soil properties for distinguishing erosion classes; (3) classification of spatial data into erosion classes; and (4) validation of results based on comparison with point soil data.Flowchart of the whole process is depicted in Figure 2.

Definition of Groups of Site-Specific Erosion Classes
Soil samples data were classified into groups based on description of soil profiles (soil profile stratigraphy, soil depth or A horizon thickness), field observation (evidence of erosion such as ploughed subsoil, increase in skeleton content) and expert knowledge of soil erosion evidences in different pedological conditions on study sites.Setting of defining parameters of each class is site-specific depending on soil type and parent material.The quantitative setting of each erosion stage at individual study sites was based on previously performed regional studies [71,72].Four groups (erosion classes of soils) were established at each study site (see Table 3): Non-eroded soils (NE)-unchanged autochthon soils with no or negligible evidence of material removal or accumulation; Eroded soils with various stage of degradation-moderately eroded (ME) and strongly eroded (SE) soils-soil profiles with evidence of material removal and soil profile truncation; and Accumulated soils (AC) formed by material re-deposition in concave parts of relief-soil profiles with evidences of accumulation of new material manifested by increased thickness of A horizon or burial of former surface horizons.

Definition of Groups of Site-Specific Erosion Classes
Soil samples data were classified into groups based on description of soil profiles (soil profile stratigraphy, soil depth or A horizon thickness), field observation (evidence of erosion such as ploughed subsoil, increase in skeleton content) and expert knowledge of soil erosion evidences in different pedological conditions on study sites.Setting of defining parameters of each class is site-specific depending on soil type and parent material.The quantitative setting of each erosion stage at individual study sites was based on previously performed regional studies [71,72].Four groups (erosion classes of soils) were established at each study site (see Table 3): Non-eroded soils (NE)-unchanged autochthon soils with no or negligible evidence of material removal or accumulation; Eroded soils with various stage of degradation-moderately eroded (ME) and strongly eroded (SE) soils-soil profiles with evidence of material removal and soil profile truncation; and Accumulated soils (AC) formed by material re-deposition in concave parts of relief-soil profiles with evidences of accumulation of new material manifested by increased thickness of A horizon or burial of former surface horizons.

Assessment of Soil Properties for Erosion Classes Distinguishing
Each of erosion classes was consequently characterized by a range of values of distinctive soil surface properties (soil erosion indicators, [30]) having specific spectral reflectance characteristics that can be obtained using imaging spectroscopy (SOC, texture classes, CaCO 3 and iron oxides).We presume that in ME and SE, the removal of topsoil and its mixture with subsoil leads to qualitative and quantitative changes in surface soil properties manifested in spectral reflectance characteristics [30,33,41,42,73].However, the soil properties of deposited material can vary significantly (both topsoil and subsoil material can be accumulated), strongly depending on site-specific conditions (topography, parent material, soil properties of parental soil, type, intensity and frequency of erosion events).This fact can significantly influence the information value of surface soil properties and the corresponding reflectance data and their applicability in defining the AC class [71,73].For this purpose, thickness of A horizon was used as an additional indicator for the definition of this class.The thickness was spatially predicted using the data from borings and derivatives of DEM (altitude, slope, curvatures, catchment area, and topographic wetness index).The 5 × 5 m 2 resolution DTM (DMR 4G ® ) distributed by Czech Office for Surveying, Mapping and Cadaster-ČÚZK, obtained using laser scanning with a total mean height error of 0.3 m was used.
Soil properties applicability to distinguish individual erosion classes was analysed using multiple range test (parametric and non-parametric, depending on statistics of normality and variability in the set).Only the properties that showed significant differences in at least two erosion classes (by means of unique confidence interval) were used for distinguishing the erosion classes.Soil properties which achieved low prediction accuracy using hyperspectral image processing (R 2 < 0.5) were not used as defining parameters.

Classification of Spatial Data into Erosion Classes
Classification of erosion stages was carried out using fuzzy C-means (FCM) clustering algorithm [74].C-means function from R package e1071 [75] was used for clustering execution.The method creates partitioning of dataset in multivariate space by maximizing similarity between samples within the same cluster and dissimilarity among different clusters.The probability of cluster membership of samples is computed.Final classification into classes is based on the highest probability.
Spatial data of previously selected soil properties derived from hyperspectral images and DEM were used as input data layers.Initial cluster centres were determined as means of soil properties values for defined erosion classes.

Validation of Results
Summarizing the validation of results was performed by confusion matrix.Overall and individual class accuracy was obtained for all study sites.Since the independent samples were not available, point data classes determined under the field campaign were used for the validation.

Descriptive Statistics of Soil Samples
Summary statistics and correlation matrix for the soil samples from the study sites are shown in Tables 4 and 5.The SOC content analysis in topsoil has revealed differences among the study sites.Low contents were observed at Jičín site (Luvisols), where values ranged from 0.7% to 1.41% (mean 1.03%), the highest variability and low to high contents were observed at Šardice site (range 0.84%-2.62%,mean 1.44%).On the other sites (Cambisols) Nová Ves (range 0.56%-1.44%,mean 1.07%) and Přestavlky (range 0.61%-1.88%,mean = 1.19%), the content of SOC is low to medium.The soil texture analysis has revealed silt loam (partially silty clay loam) at Jičín site, prevalence of loamy soils at Šardice site (locally more sandy or clayey), sandy loam and loam at Nová Ves and silt loam and loam at Přestavlky.CaCO 3 content was analysed only on sites with loess parent material.At Jičín site, presence of CaCO 3 was identified only in 3 topsoil samples (~0.1%).At Šardice site, CaCO 3 content varies from 0% to 10%.Fe oxides in soils extracted by two reagents (Fe ox and Fe dith ) were analysed only on Cambisols sites.The content of both was identified higher at Přestavlky site (Fe ox mean 0.77%; Fe dith mean 3.59%) than at Nová Ves site (Fe ox mean 0.21%; Fe dith mean 1.24%).

Prediction of Soil Properties by Imaging Spectroscopy
Different multivariate techniques (ANN, PLS, RF, SVM) were used to predict soil properties spatial distribution using raw spectroscopic data and data adjusted by different pre-treatment methods (log, SG 1st, SG 2nd, CR and SNV).The results of multivariate modelling indicate different prediction accuracy according to soil properties, number of soil samples used for calibration, study sites, used multivariate technique and pre-treatment methods.Table 6 summarizes the results of calibration and validation by prediction models with the highest accuracy according to lowest RMSE P values.Values of R 2 and RMSE in training and validation sets for all models and different types of pre-treatment are shown in Figure S1 (in supplementary material).Plots of measured vs. predicted soil parameters for the validation datasets are given in Figure 3.
In summary, the best prediction accuracy was achieved for SOC content, where R 2 value in the validation set ranged for the best models within sites between 0.8 and 0.91.Consequently, these models yielded RMSE P between 0.07% and 0.12%.SOC is the most frequently predicted soil property using hyperspectral imaging data.Reported values of R 2 and RMSE for SOC prediction from previous works are highly variable [37] and depend on many local conditions and used techniques (type of sensors, corrections, site conditions, size of study area, number of samples etc.).Nevertheless, the prediction of SOC distribution is generally successful thanks to good spectral response of SOC and the values reported in similar studies (R 2 p ranged between 0.65 and 0.96) are comparable to our results [76][77][78][79].
Prediction accuracy of textural classes is highly variable among study sites and more variable than in SOC prediction.R2 p and RMSE P values vary between 0.21-0.67R 2 p and 2.49%-9.04%RMSE P for sand content, 0.4-0.92R 2 p and 2.75%-7.06%RMSE P for silt content and 0.38-0.89R 2 p and 1.47%-2.88%RMSE P for clay content.In general, lower value of R 2 p for certain textural class was achieved at the sites where its content is low.This is, for example, the case of Jičín site where a very low prediction accuracy of sand (R 2 p 0.21) was assessed.Alternatively, the prediction accuracy of sand and silt which are not spectrally active in VNIR-SWIR region can be influenced by a correlation between these properties and spectrally active properties (SOC, clay, CaCO 3 and Fe forms).This is probable at Jičín site where a good model performance for silt (R 2 p 0.92) can be given by its strong correlation with clay (R 2 p 0.89) or at Šardice site where sand (R 2 p 0.67) is highly negative correlated with clay (R 2 p 0.89) (Table 5).Another factor influencing the variability of the prediction accuracy can be a different level of BRDF effect caused by diverse hyperspectral sensing conditions.This is valid for Přestavlky and Nová Ves sites where the BREFCOR corrections have not led to full elimination of BRDF effect.A high variability of model performances was reported in several studies dealing with texture prediction by RS [80][81][82][83][84][85].[84] or even better [86,87] in comparison to other studies.In conclusion, the analysis showed that the distribution of soil properties in eroded landscapes can be successfully predicted using the spectroscopic data.The prediction accuracy is adequate in majority of sites and predicted properties; however, a significant worsening in the model performance was observed in case of properties with low content [36].Moreover, the presented models are locally-trained (for few field blocks in maximum) which increase their prediction capability in comparison with similar regional models [76].In conclusion, the analysis showed that the distribution of soil properties in eroded landscapes can be successfully predicted using the spectroscopic data.The prediction accuracy is adequate in majority of sites and predicted properties; however, a significant worsening in the model performance was observed in case of properties with low content [36].Moreover, the presented models are locally-trained (for few field blocks in maximum) which increase their prediction capability in comparison with similar regional models [76].Methods of machine learning SVM (71% of cases) and ANN (19%) showed the best prediction accuracy.At Šardice site, the PLS model (10%) showed the best accuracy for two textural classes (similar to machine learning models).PLS technique has been traditionally used in many studies dealing with soil spectroscopy and its ability to prediction of soil properties from spectral data was proven, but with varying accuracy [36].In this study, PLS models reveal slightly lower accuracy than other methods, however not significantly worse.More recently, machine learning techniques have been used in imaging spectroscopy [30,76,[88][89][90][91]. Whereas ANN method is not so frequently used due to the need of a large training dataset, training demands and a tendency to over-fitting, the SVM has become the most frequent method overcoming the difficulties of the ANN technique [27].Contrarily to the mentioned ANN shortcomings, the tendency to over-fitting by ANN was not fully proven in our study.The problem with over-fitting was observed in RF in all cases.This finding is in contrast to Viscarra Rossel and Behrens [92] and it is most probably influenced by insufficient optimization of model parameters by random search method for automatic search of parameter values.General applicability of SVM as a robust model with low sensitivity to noise in data [76,89] and its ability to predict soil properties across hyperspectral images with small amount of samples was proven in our study as the SVM was the method of the first choice in the majority of predicted properties at the study sites.
Regarding the spectra pre-treatment methods, significant differences of prediction accuracy were not observed.The majority of the best performing models (33%) used the techniques of absorbance transformations.Raw data (reflectance), SG (1st), and SG (2nd) were used in 19%, CR and SNV in 5% of the finally selected models.The methods of pre-treatment performed better (according to R 2 ) than reflectance and absorbance data in case of SOC and Fe ox prediction at Přestavlky and Nová Ves sites, where the conditions of hyperspectral data acquisition were less favourable, which led to stronger influence of BRDF effect.This effect can be highlighted at these sites by cementing effect of iron oxides correlated to soil aggregation [36,93].
In case of other properties, the relationship between prediction ability and the use of a pre-treatment method was not observed.Spatial distribution of all predicted soil properties for all sites is shown in Figure 4.

Assessment of Soil Erosion Classes
Four erosion classes were distinguished on the study sites.The highest number of profiles was identified as non-eroded (11 at Šardice, 29 at Přestavlky, 30 at Jičín and 40 at Nová Ves) except for the Chernozem site Šardice where 23 profiles were classified as strongly eroded soils.Moderately eroded soils were identified in 5-8 cases at each site.Accumulated soil class was identified only at Šardice (9) and Jičín (12) sites.

Assessment of Soil Properties for Erosion Classes Distinguishing
Results of selection and assigning soil properties to each erosion class are shown in the Table 7.The table shows selected soil properties and their centres (means) used for erosion classes classification.The performance and suitability of different properties for erosion classes' classification are highly variable at the study sites.In summary, the SOC content is applicable at every site except for Jičín, where the difference in the values among erosion classes was not significant.At three study sites (except for Jičín), SOC can be used for distinguishing of NE and SE class with significantly the highest and lowest SOC content, respectively.Potential of SOC for ME classification is low.In contrast to Schmid et al. [30] a Hill and Schütt [42] who assessed the erosion stages in semi-arid climate, our study showed difficulties in distinguishing the AC class using SOC content in topsoil.At Šardice site, the SOC content in topsoil decreases in AC class.This is due to an advanced stage of erosion degradation at the study site when loess is exposed in eroded parts of slopes and redeposited in the accumulation positions [53,72].Zádorová et al. [94] described this process as a retrograde soil development typical for the most vulnerable loess regions with dissected relief.Moreover, the combination of preferential (interill erosion) and non-preferential (rill and tillage erosion) removal of soil material influences the resulting variability of SOC distribution in topsoil [95-97].In case of texture, the classification potential of at least one textural class was proven at all the study sites.However, the accuracy of prediction model for these properties did not reach required limit of R 2 ; therefore, only sand at Přestavlky, silt at Nová Ves and Jičín and clay at Jičín were used for classification.At Přestavlky site, SE was significantly different from NE and ME in sand content.This result can be expected at sites with high sand content where the fine particles are selectively removed from the eroded soils.At Nová Ves, high silt content can be used for distinguishing NE class.AC class was significantly distinguishable from NE at Šardice site by clay content.Iron oxides applicability was observed only at Nová Ves site.
CaCO 3 content was analysed only at Šardice site (only negligible amount was observed at Jičín site), where it showed a good potential to distinguish AC class from NE and ME class.The AC and SE classes were not significantly different, as the CaCO 3 content is similarly high due to truncation of soil profile and admixture of loess in the plough layer in case of SE class and accumulation of CaCO 3 -rich material in case of AC class.The study showed that the CaCO 3 content can be used as an erosion classes indicator only in regions with CaCO 3 -rich subsoil that is exposed at soil surface due to A horizon removal or tillage, similarly to [30,33].
Due to poor results of tested soil properties in AC class distinguishing, we used A horizon thickness as an auxiliary indicator at Šardice and Jičín sites, where the AC class was identified.The A horizon thickness prediction using DEM modelling was previously successfully used in Zádorová et al. [71,72] for identification of deep colluvial soils.In our study, the A horizon thickness showed very good potential in AC distinguishing at both sites.It was also significant for NE class in Šardice; at Jičín site, it did not differ significantly in any other class except for AC.

Classification of Spatial Data into Erosion Classes
Results of classification based on FCM and previously determined class centres are shown in Figure 5.In total, 1.24 km 2 of strongly eroded soil was classified.This class is arranged in descending order; 32.9% at Nová Ves, 29.1% at Šardice, 23.3% at Jičín and 16.9% at Přestavlky.At Přestavlky site, the NE class (46.2%) followed by ME class (36.8%) cover the rest of the area.At Nová Ves site, the area of moderately eroded soil (35.4%) exceeded the area of non-eroded soils (31.7%).At both localities, SE class can be found only on the most exposed terrain position (in context of both water and tillage erosion).The pattern of different erosion classes is more complex at Šardice and Jičín sites.The AC class is present mainly in the terrain concavities (side valleys) closely neighbouring to strongly and moderately eroded soils at the back slopes (20.7% at Šardice; 13.6% at Jičín).Non-eroded soils were classified in flat and nearly flat relief (plateaus) (22.0% at Šardice; 40.5% at Jičín).These findings are in accordance with previous studies performed in regions with similar conditions.Zádorová et al. [94], Schmitt and Rodzik [98] or Terhorst [99] observed a high heterogeneity of soil cover, depth and SOC stocks and a significant area covered by strongly eroded soils in Chernozem and Luvisol loess regions.The Cambisol regions on crystalline rocks were typical in terms of less pronounced erosion evidences with lower incidence of both accumulated soils and soils with severe erosion removal [71].
strongly and moderately eroded soils at the back slopes (20.7% at Šardice; 13.6% at Jičín).Non-eroded soils were classified in flat and nearly flat relief (plateaus) (22.0% at Šardice; 40.5% at Jičín).These findings are in accordance with previous studies performed in regions with similar conditions.Zádorová et al. [94], Schmitt and Rodzik [98] or Terhorst [99] observed a high heterogeneity of soil cover, depth and SOC stocks and a significant area covered by strongly eroded soils in Chernozem and Luvisol loess regions.The Cambisol regions on crystalline rocks were typical in terms of less pronounced erosion evidences with lower incidence of both accumulated soils and soils with severe erosion removal [71].

Validation of Results
The overall accuracy of erosion classes' classification varies across the sites (Table 8).High accuracy was obtained in Šardice site (82%), moderate accuracy was obtained in Jičín site (67.3%).Přestavlky and Nová Ves achieved 51.1% and 52.6%, respectively.This is regarded as a low accuracy.Significant deficiencies of the spectral modelling were expressed in case of classification of

Validation of Results
The overall accuracy of erosion classes' classification varies across the sites (Table 8).High accuracy was obtained in Šardice site (82%), moderate accuracy was obtained in Jičín site (67.3%).Přestavlky and Nová Ves achieved 51.1% and 52.6%, respectively.This is regarded as a low accuracy.Significant deficiencies of the spectral modelling were expressed in case of classification of two classes of eroded soils (SE and ME).The most frequent classification error is in all sites linked with misclassification of moderately eroded class to non-eroded or to strongly eroded class or vice versa.In general, the moderately eroded class achieved lowest values of both producer and user accuracy.This problem was also reported in Schmid et al. [30] who identified the majority of misclassifications in the moderately eroded class.According to the results of classes' separability and model validation, distinguishing two classes of eroded soils shows to be too detailed and classification of one class of eroded soil would be sufficient as the two classes overlap in a number of soil properties.Distinguishing transitional classes can be performed using the fuzzy classification and comparison of membership in non-eroded and eroded classes.However, distinguishing of different classes of eroded soils can be reasonable and bring good results (77% in overall accuracy) in certain conditions as was proven by Schmid et al. [30] in a large, severely degraded semi-arid region on CaCO 3 -rich parent material.Zádorová et al. [94] reported similar degree of misclassification in case of slightly and strongly accumulated soils at a Chernozem study site.The slightly accumulated soils showed the lowest overall accuracy due to aggregating the properties of other classes.At Přestavlky site, two SE samples were classified as NE.The profiles are situated on the slope with high amounts of rock fragments on the surface (up to 20%) and evidence of erosion features; the misclassification to NE class is evident.At Nová Ves site, 4 NE soil profiles were classified as SE.The pattern of eroded and non-eroded soils is very complex at this site, strongly depending on the tillage erosion.The site is distinctive by abrupt change of soil depth within short distances.Misinterpretation of erosion classes at this site can be therefore associated to an insufficient spatial accuracy of input data (GNSS, hyperspectral, DEM) for the exact delineation of the soil variability at short distances.
Few misclassifications between SE and AC classes were observed at Šardice (1) and Jičín (5) sites.The misclassified points were situated near the inflex points of the slope where removal, transport and sedimentation can act in dependence to the intensity and the type of erosion.The resulting surface soil properties correspond to such a complex process and can overlap in the studied classes.This fact is proven by the fuzzy membership that is evenly distributed to more erosion classes.Two samples in the north part of Jičín site identified as AC were predicted as SE.This error can be again explained by the mixed effect; the points are situated in the narrow transitional strip between Luvisol and gleyic Chernozem nearby the stream.
The validation of results showed high to moderate potential of used methods for erosion classes' classification.In our study, the high applicability of spectral data is restricted to Chernozem and Luvisol loess study sites (Šardice and Jičín).This can be attributed to two facts.The first is the homogeneity of parent material and soil properties of former dominant soil cover at these sites, now changed by erosion.The observed differences in the soil properties are then closely linked to erosion processes and thus well performing in the prediction models.In contrast, the variability of soil properties at the sites on heterogeneous crystalline and sedimentary rocks with more heterogeneous soil cover (Přestavlky and Nová Ves) is a result of more co-acting processes and local conditions and is not necessarily linked to the soil erosion.The second factor is the intensity of the soil cover change due to erosion processes.Erosion leads to extreme (Šardice) and significant (Jičín) redistribution of soil material at the loess sites and results in development of severely degraded soils with truncated soil profile on one hand and several metres deep colluvial soils on the other hand [53].The particular erosion classes are then very distinctive and form a specific soil mosaic.The intensity of soil erosion at Cambisol study sites is less pronounced; the accumulation class has not been identified and SE class covers marginal areas.Thus, the soil properties variability given solely by soil erosion is presumably less significant and the soil pattern driven by erosion material transport is less evident.Similar results were reported by [71] who compared the soil units and SOC stock prediction DEM based models in Chernozem, Luvisol and Cambisol study sites with similar area.They reported a very good model performance in Chernozem and Luvisol sites and a poor performance at the Cambisol site.The studies from semiarid regions dealing with soil erosion stages prediction using spectral data [30,33] showed high accuracy of the model, similar to the results at Šardice site.However, the comparison is not fully relevant according to differences in site area and conditions.

Conclusions
The study demonstrates the potential of different soil properties predicted using hyperspectral data for the assessment of soil degradation by erosion at four pedologically different study sites influenced by soil redistribution due to accelerated erosion.The erosion impact at the study sites was evaluated by distinguishing four erosion classes representing different stages of soil redistribution by erosion.
The study showed that: (i) soil properties prediction can be successfully performed using spectral data and adequate prediction method; and (ii) selected soil properties are applicable for the assessment of soil degradation by erosion.The selected predictive properties, best performing predictive methods and classification models accuracy differed in the study sites.The accuracy of classification models was influenced by variability of soil units and parent material.The presented approach was successfully applied in Chernozem and Luvisol loess regions where the erosion classes were assessed with good overall accuracy (82% and 67%, respectively).The model performance in two Cambisol regions was rather poor (51%-52%).At study sites with less pronounced soil degradation, the restriction of the erosion classes from four to three can bring better results.However, at severely degraded sites, the limited number of classes may lead to a significant loss of information and decrease in models' applicability in conservation management.The sites with heterogeneous soil properties and parent material will require more precise local-fitted models and use of further auxiliary information such as terrain or geological data.
The key requirement for a successful use of hyperspectral data in soil predictive modelling is the application of high quality input spectral data with precisely performed corrections and limited influence of physical factors.The approach presented in the study can be applied exclusively on bare soil.Thus, the direct observation is in case of temperate agricultural regions limited to periods with minimum vegetation cover.More images from different periods are needed to cover wider areas.The study sites selected for the study and the periods of viewing were chosen with aim to minimize the roughness, wetness and non-photosynthetic vegetation cover.Further research is needed to find effective methods to filter out these factors and facilitate the use of hyperspectral imaging at regional scale, especially in context of increased use of no tillage farming.Implementation of satellite hyperspectral sensors as well as further research on application of multi-and super-spectral sensors represent other effective tools in application of presented approach.
The applicability of the presented approach at the global scale has considerable limitations.To develop and apply a global model of a reasonable accuracy, utilizable for management purposes, similar conditions (atmospheric, surface roughness etc.) within the flight campaigns or a very precise image pre-processing (namely the elimination of the BRDF effect) are needed.Thus, regional models performed at a watershed scale or at a level of a geologic/pedologic/agricultural region represent, from the practical point of view, a more applicable option as a reasonable accuracy of the models can be provided more easily.
The presented approach promise, mainly at the local a regional scale, to produce valuable data on actual soil degradation, present structure of soil cover and redistribution of soil properties due to soil erosion, that will be usable for soil conservation policy purposes.The acquired data can be directly used by farmers for adjustment of management practices according to the level of soil degradation by erosion.At present, such information on real erosion impact on different soils is not available in the required detail and quality.

Supplementary Materials:
The following are available online at www.mdpi.com/2072-4292/9/1/28/s1, Figure S1: R2 and RMSE of prediction models for training and validation datasets by study site, multivariate technique and pre-treatment method.

Figure 1 .
Figure 1.Location of the study sites (yellow border) with position of hyperspectral images.

Figure 1 .
Figure 1.Location of the study sites (yellow border) with position of hyperspectral images.

Figure 2 .
Figure 2. Scheme of the method for soil erosion assessment.

Figure 4 .
Figure 4. Spatial distribution of soil properties at the study sites based on multivariate prediction using the best performing model.

Figure 4 .
Figure 4. Spatial distribution of soil properties at the study sites based on multivariate prediction using the best performing model.

Figure 5 .
Figure 5. Maps of erosion classes at the study sites as derived by means of fuzzy C-means method using the mean values of selected soil properties as class centres: (a) Přestavlky; (b) Šardice; (c) Nová Ves; and (d) Jičín.

Figure 5 .
Figure 5. Maps of erosion classes at the study sites as derived by means of fuzzy C-means method using the mean values of selected soil properties as class centres: (a) Přestavlky; (b) Šardice; (c) Nová Ves; and (d) Jičín.

.
Study sites and their characteristics.

Table 1 .
Study sites and their characteristics.

Table 2 .
Details of individual aerial hyperspectral campaigns.

Table 3 .
Erosion classes distinguished at each site according to soil profile stratigraphy and thickness of A horizon.
Figure 2. Scheme of the method for soil erosion assessment.

Table 3 .
Erosion classes distinguished at each site according to soil profile stratigraphy and thickness of A horizon.

Table 4 .
Descriptive statistics of soil properties measured in collected soil samples.

Table 5 .
Correlation matrix of soil properties.
* correlation is significant at the 0.05 level; ** correlation is significant at the 0.01 level; *** correlation is significant at the 0.001 level; non-signed values are non-significant.

Table 6 .
Accuracy of the best prediction models (the best performing model according to lowest RMSE P values).Models dealing with iron oxides and CaCO 3 content performed with R 2 p 0.73-089 for Fe ox , 0.59-0.78for Fe d and 0.82 for CaCO 3 (only at site Šardice).Fe ox , Fe d (iron oxides in general) and CaCO 3 are only rarely studied soil properties by imaging spectroscopy.The high to moderate prediction accuracy for Fe ox , Fe d and CaCO 3 is similar

Table 7 .
Class centres used for determination of erosion classes by fuzzy C-means method.

Table 8 .
Confusion matrix of the erosion stages classification.