Point and pixel inclusive machine learning models for exploring gully erosion susceptibility

Abstract Sample point based spatial model derived from Machine Learning (ML) algorithms often generalizes the spatial pattern of an event. The present study has tried to highlight how far it is acceptable and can it replace with pixel based modeling? The present study has presented a comparative view of pixel and sample point based modeling of gully erosion susceptibility of the upper Mayurakshi basin to assess the predictabilities. Random forest (RF), Support Vector Machine (SVM), and (ADB) have been applied for developing pixel and point based models in Python and WEKA software environments respectively. The models show that 14–25% area located mainly in the upper parts of the study unit is very highly susceptible to gully erosion. Based on the accuracy and performance level using area under curve (AUC) of Receiver operating curve, sensitivity, precision, F1 score, MCC, pixel based ensemble models are superior to point based modeling. The point-based models often have an inferior agreement with training and testing data. So, pixel based models could not be replaced with point based models. RF is found as the best representative model. The study, therefore, recommends using pixel based modeling for this or a similar purpose. Since the models have figured out the gully erosion susceptible areas, it would be a useful tool for related planning processes.


Introduction
Machine Learning (ML) algorithms based on spatial models coupled with and Geographical Information System (GIS) has been used to address different issues in a wide array of disciplines, including epidemiology (Auchincloss et al. 2012), climate change (Mansfield et al. 2020), natural resource (Frey 2020), environmental vulnerability (Pal and Debanshi 2021), environmental hazards . Gully erosion is one of the major issues that have become a global concern with the increasing coverage of degraded land worldwide (Secretariat of the UNCCD 2013). The global coverage of degraded land is about 30% (Nkonya et al. 2016) and about 44% in India (Mythili and Goedecke 2016). Soil erosion by water, especially in the gully, is the most common soil degradation process (Saha et al. 2020;Panagos et al. 2021). Gully erosion featuring deep channels eroded by running water causes that degradation locally by forming badlands over lateritic terrain (Tang et al. 2013).Gully erosion is capable of causing acute problems by causing the formation of badlands (C anovas et al. 2017;Wuepper et al. 2021), eroding fertile soil of agricultural land (Han et al. 2018), clogging sediment in a reservoir (Dutta 2016), depleting of groundwater table (Tilahun et al. 2016). Thus, sustainable practices for soil conservation should be adopted to secure future development (Vanmaercke et al. 2016;Arabameri et al. 2020). The present study area, the Mayurakshi river basin, suffers from the severe problem of gully erosion due to its geo-environmental location over the degraded eastern fringe of the Chottanagpur plateau (Bhattacharyya 2014;Debanshi and Pal 2020). Gully erosion susceptibility mapping is essential for initiating such conservation practice (Zabihi et al. 2018;Debanshi and Pal 2020). A gully erosion susceptibility mapping can be done by detecting the relationship between the existing gully occurrences and geo-environmental gully conditioning factors (Nicu 2021;Rahmati et al. 2017).
ML algorithms are becoming increasingly applied tools for assessing gully erosion and are considered robust over the traditional statistical tool for their accurate predictability (Azareh et al. 2019;Arabameri et al. 2020;Saha et al. 2020). Several software and programming languages like R studio (Arabameri et al. 2018;Hosseinalizadeh et al. 2019a), Waikato Environment for Knowledge Analysis (WEKA) (Hosseinalizadeh et al. 2019b), Python (Wang et al. 2020;Can et al. 2021) are used for implementing the ML algorithms. Several studies reported the robustness of ML models over the multi-criteria decision approach models, especially when feature space is complex and the target variable is involved with the input dataset in a non-linear relationship (Rodriguez-Galiano et al. 2012;Churpek et al. 2016;Xie et al. 2018;Jung and Lee 2019). ML algorithms have the potential to identify the complex non-linear relationship (Mosavi et al. 2018) and eventually produce the spatial model of gully erosion susceptibility based on associated predisposing factors of gully erosion (Saha et al. 2020;Borrelli et al. 2022). Integrating remote sensing and GIS with ML and ensemble algorithms has made mapping any environmental issue, including gully erosion, easier and more precise (Pal and Paul 2020;Saha et al. 2020;Rafiei Sardooi et al. 2021). For example, Arabameri et al. (2018) and Hosseinalizadeh et al. (2019a) used R programming language and GIS technique to assess gully erosion susceptibility and reported greater performance accuracy and preciseness than WEKA is a software that provides a set of ML and ensemble algorithms and leverages data mining and model building (Ahmadpour et al. 2021). Many studies have used the WEKA tool and GIS for developing spatial models (Oh et al. 2019;Hosseinalizadeh et al. 2019b;Pham et al. 2020). Since the software provides a set of already programmed algorithms, the parameters can be optimized to achieve the best performance Pham et al. 2020). Recently, the python programming language and associated libraries are also used for developing ML-based predictive spatial models (Wang et al. 2020;Arabameri et al. 2021;Can et al. 2021).
In this context, in recent times, the generalization of the spatial predictive model developed based on ML and Deep Learning (DL) has been in discussion (Maxwell et al. 2016;Neyshabur et al. 2017;Kawaguchi et al. 2017;Maxwell et al. 2020aMaxwell et al. , 2020b. This generalization was sometimes the different geographical regions from which the model has been trained and applied (Maxwell et al. 2020b) or sometimes the disparate data set . However, in most cases, the spatial models integrated with ML algorithms and GIS limit the type of analysis performed (Hogland and Anderson 2017). The process of such model building involves multiple sequential steps viz. (1) building dataset using GIS, (2) importing the dataset in ML algorithm performing software, (3) training of the model using a set of samples of that imported dataset where the target variable is known and define a relationship between target variable and explanatory variable, (4) applying the model on the rest of the dataset, and then (5) producing a spatially explicit surface in GIS environment using the output of the predictive model. Producing final outputs in this manner is often challenged by multiple issues like learning multiple software, importing and exporting a large dataset, long processing time, and the requirement for more extensive storage space Anderson 2014, 2015). While producing spatial models of a larger geographical area, these limitations intensify, mainly because of the issue of storage space. Since the remote sensing datasets are becoming more refined day by day, a large dataset is needed to be handled to take all the pixels of a large geographical region into account (Hogland and Anderson 2017). Model building using a previously programmed set of algorithms like the WEKA tool suffers from limitations due to its lesser storage capacity (Naudts et al. 2004). To overcome this issue, sample point based model building is practised. This method has a fast model building benefit but often gets subjected to generalizing or smearing the nuances of the place located further from the chosen sample sites or points. This problem reduces the precision of the mapping, and eventually, the initiation of the conservation practice is compromised. This is an acute problem and could be resolved by taking all the pixels of the study area into consideration while producing any model. Python libraries give such an opportunity to handle large datasets containing all the pixels of a larger geographical region. It is an open space, and the algorithms need to be programmed (Brownlee 2016), making it difficult for non-expert users. In addition, machine learning needs hyper-parameter optimization to increase the prediction accuracy, which requires additional computation (Kotthoff et al. 2019). Hyper-parameters are the parameters that are set before the training process, and optimization can be achieved by manual search or automatic search method . The auto-optimization function of the WEKA tool, to some extent, eases the hyper-parameter optimization process (Kotthoff et al. 2019), but it has limited scope, and it optimizes the hyper-parameters up to 70% (Brownlee 2019). Python library enables to overcome this limitation with a broader spectrum of hyper-parameters optimization (Brownlee 2016). Therefore, python libraries are considered superior for optimizing the hyper-parameters than WEKA (Mitrpanont et al. 2017). Here the research gap could be pointed out whether the time and expertise expensive laborious work of pixel inclusive modelling are essential to replace the representative sample point based modelling for spatial analysis. If yes, then to what extent the precision of the pixel inclusive models greater than the time effective and fast approach of representative sample point based modeling. To answer this query, the pros and cons of sample point based modeling using WEKA tool and pixel inclusive modeling using Python libraries have been decided. Considering this issue and addressing the aforementioned research gap, the present study has intended to examine the better credential of Sample Point Based Modeling (SPBM) and Pixel Based Modeling (PBM) to explore gully erosion susceptibility.

Study area
The Mayurakshi river basin is located over the Chottanagpur plateau fringe region of eastern India (Figure 1a) between 23 15 0 N to 24 34 0 15 00 N latitude and 86 58 0 E to 88 20 0 30 00 E longitude and covers more than 5403 km 2 (Figure 1b and 1c). The entire basin consists of the Rarh tract dominated by secondary lateritic formation. The rivers that originated from the Chottanagpur Plateau took part in such formation Kapat 2003, 2009;Pal 2016). Geologically (Figure 1b), Granitic gneiss is a dominant rock type that belongs to the Pleistocene age that covers most of the upper and middle catchment of the basin (Figure 1e and 1f). Granet-biotite gneiss is also found over a considerable area in the north-eastern part of the basin. Geomorphologically (2c), the upper and middle catchment of the basin consists of a denudational origin pediment-pedeplain complex, and in a few patches, moderately dissected hills and valleys originating from the denudational process are found. On the other hand, lower catchment mostly belongs to flood plain geomorphologically comprised of fluvial origin older alluvial plain, older flood plain, and active flood plain. Geologically also, this part of the basin belongs to newer and older alluvium. Eventually, this part of the basin was not subjected to severe soil erosion; instead, an alluvial deposition by the western tributaries took part behind the formation of this part (Bandyopadhyay et al. 2014). Apart from that, the confluence part of the eastern side of the basin is merged with the Ganges delta. Almost all the gully occurrences are found in the upper and middle part of the basin. In many cases these gullies appear as 1st and 2nd order streams. Araujo and Pejon (2015) and Joshi et al. (2016) reported the same in their study area also. In this part the numbers of 1st and 2nd order streams are respectively 4827 and 1091 (extracted from SOI topographical maps). Geologically this part belongs to Granitic gneiss formation and Granet-biotite gneiss formation and geomorphologically belongs to denudational origin pediment-pedeplain complex and denudational origin moderately dissected hill and valleys. Therefore, these geological and geomorphological divisions have been taken as an area of interest in this study, and the lower part of the basin has been excluded from the investigation (Figure 1b-d). This area of interest covers more than 3500 km 2 which is about 65% of the entire basin area. Figure  2 presents the channel network of some gully erosion prone sites.

Materials
For this study, 18 data layers have been prepared using various data sources like USGS Landsat images with 30 m spatial resolution (path/row: 139/43 and 139/44) to prepare land use/land cover (LULC), Ferrous mineral index, Normalized Difference Vegetation Index (NDVI), Bare Soil Index (BSI), Water Supplying Vegetation Index (WSVI) and monthly fluctuation of a temperature map of the basin area. The map of geological divisions has been taken from the Geological Survey of India, whereas BHUVAN (an Indian web-based utility) has been subscribed for extracting the geomorphological map of the study area. Mouza-wise soil texture map has been derived from National Informatics Centre (NIC) to prepare a soil erosivity map. The average monthly and annual rainfall map of the Department of Science and Technology (DST) has been employed in this study for rainfall information and its further derivations. The spatial data layers of the slope, aspect, and stream power index (SPI) have been prepared using SRTM DEMs with 30 m spatial resolution. Survey of India (SOI) toposheets of 1972 (scale 1:50,000) have been used for developing the base map of the study, along with preparing the maps of the drainage network and water covered area and extracting the location of gully head bundh ( Figure 3). Recent Google Earth imagery (22/10/2021) has also been subscribed for extracting the location of check dams. All Landsat imageries and SRTM DEMs have been preprocessed before final layer preparation.

Software taken for PBM and SPBM
Generally, SPBM is considered a cost-effective and quick approach for environmental studies (Chabuk et al. 2020). Representation of the surface dynamics through these sample points minimizes the variability between two subsequent points and generates specifiable continuity in the surface (Zhu et al. 2020). PBM is another approach where all the pixels are considered for the generation of the final output (Riyadi et al. 2020). In pixel based modelling, the model has been applied to a standard dataset (70% for training and 30% for validation), and the final trained model has been applied to the entire dataset or pixel (Zhu et al. 2020). There are different platforms available for sample point-based modelling (SPBM) and pixel based modelling (PBM), such as Google Cloud ML Engine, WEKA, Amazon Machine Learning (AML), Shogun, Apache Spark MLlib, Pytorch, Keras, TensorFlow, R-studio (Gonz alez et al. 2020). These platforms run various programming languages like Python, Java, Cþþ, and R to run different machine learning algorithms. Recently, Python and R have been used to model various environmental phenomena for efficiency and big data handling capability. Whereas software like WEKA and KNIME is used to run relatively small sample point datasets. For this study, we have used WEKA for sample point based modelling (SPBM), and Python has been used for pixel based modelling (PBM).

Preparing data layers for contributing factors of gully erosion susceptibility
Total eighteen factors under four-factor clusters like (1) erodibility (soil texture, LULC, geological aspect, BSI, ferrous mineral, geomorphology) (2) erosivity (soil erosivity, yearly average temperature, monthly fluctuation of temperature, WSVI) (3) resistance (NDVI, distance from GHB, water coverage, distance from check dams) and (4) topographical factors (slope, curvature, SPI, distance form 1st and 2nd order stream) have been taken for gully erosion susceptibility mapping. These factors have been selected based on their importance in determining the gully erosion of the study region and the taking experience from the existing literature like Debanshi and Pal (2020) and Saha et al. (2020). Spatial data layers have been prepared for each factor. After preprocessing the SRTM, DEM has been employed to derive the data layers of slope, curvature, and SPI. Several studies (Meliho et al. 2018;Shit et al. 2020) have used these parameters for measuring soil and gully erosion. Slope, SPI directly influences the erosion rate and also, to some extent, determines the distribution of gullies (Kert esz and Gergely 2011). Gullies expand rapidly over the steep slope compared to flat land (Marden et al. 2012). The SPI indicates the erosive power of water, considering the discharge proportional to the catchment area (Debanshi and Pal 2020). The surface analysis tool of ArcGIS software (v-10.2) has been used to facilitate slope and curvature derivation (Figures 3r and 2o), while the raster calculator tool has been used to calculate SPI (Figure 3p) of the basin area with the help of Eq. (1) (Conforti et al. 2011). Similarly, satellite imagery of the Landsat 8 series has been used to derive the data layers of LULC, BSI, ferrous mineral, monthly fluctuation of surface temperature, NDVI, and WSVI. LULC of the region can influence the gully development in many ways (Gelagay and Minale 2016). The bare ground is more prone to gully development than a dense forest (Feizizadeh et al. 2013). Supervised classification based on maximum likelihood classifier has been carried out in ERDAS Imagine software (v-9.2) for LULC identification and mapping ( Figure 3b). Since the bare ground positively contributes to gully development of its lower resistance rate (Jahantigh and Pessarakli 2011). Equation (2) has been implemented in the raster calculator for deriving BSI (Elfadaly et al. 2017) (Figure 3d). The inbuilt function of the ERDAS Imagine software has been used to produce the data layer of ferrous content ( Figure 3e). It has been recognized as a predictor of gully erosion because of the capability of ferrous minerals to initiate chemical weathering in the soil. Studies like Kapat 2003, 2009) recognized this factor for gully development in this region. The temperature fluctuation is the key factor for the mechanical disintegration of rock (Bandfield et al. 2011), and it could be an important predictor of gully development. Monthly Land Surface Temperature (LST) has been calculated (Eq. (3)) following the guideline provided by Landsat Project Science Office (2002), and its Coefficient of Variation (CV) has been calculated for the spatial measurement and preparation of monthly fluctuation of surface temperature data layer ( Figure 3i). Keeping the ability of denser forests to protect the soil surface from being eroded, the NDVI has been frequently incorporated (Arabameri et al. 2018;Roy et al. 2020;Zhou et al. 2021) in gully erosion studies. Calculation of NDVI (Eq. (4)) (Townshend and Justice 1986) and spatial mapping (Figure 3k) have been performed in the raster calculator of ArcGIS. The spatial data layers of LST and NDVI have been used to calculate the WSVI (Eq. (5)), which has been incorporated in this study as an indicator of drought severity (Figure 3j) (Alshaikh 2015). 1st and 2nd order streams often behave like gullies and accelerate soil and gully erosion (Araujo and Pejon 2015;Joshi et al. 2016). To arrest the erosion of soil and its deposition in the reservoir of the Massanjore dam which was constructed for hydro-electric generation, total of 242 smaller reservoirs have been installed at different parts of the basin. These reservoirs feature concrete barrier (locally called bundh) in front of it mostly constructed at the head of the gully and therefore, called Gully Head Bundh (GHB). These GHBs and check dams are capable of arresting soil erosion promoted by 1st and 2nd order streams (Pal and Debanshi 2018). The Euclidian distance tool of Arc GIS software has been employed to compute the spatial distance of the basin area from 1st and 2nd order streams, GHBs, check dams, and generate the spatial data layers at a 30 m Â 30 m pixel scale (Figure 3q, 3l, and 3n). 1st and 2nd order streams have been digitized; the locations of GHBs have been extracted from the SOI toposheets, while the locations of check dams have been extracted from the Google earth imagery. The frequency of existing gullies per unit area of each geological division (scale-1:200,000) has been calculated, and a raster layer with a 30 m Â 30 m cell size has been generated for each division, considering its gully frequency for facilitating the input of spatial model development. The same procedure has been repeated for geomorphological divisions also. For producing the spatial data layer of soil texture (Figure 3a), the proportion of coarser sand in the top soil has been considered since coarse-textured is very prone to gully formation and erosion (Pal 2016). The proportion of coarser sand in the top soil layer of the upper catchment has been measured by extracting soil samples from 43 sites and testing with the help of a digital sieve shaker. The village level soil texture maps (1:100,000) of NIC have also been subscribed to get the soil texture data. The percentage of coarser sand in the soil surface has been identified based on USDA (1999) soil classification system. The DST-provided rainfall map has been used for preparing the data layer of average yearly rainfall (Figure 3h). Soil erosivity has also been incorporated in this study to consider the impact of rainfall intensity on the soil. Since the concentrated rainfall in fewer months causes higher erosion than the uniform rainfall for the entire year (Pal and Debanshi 2018), the erosivity of the rainfall has been considered as a parameter in this study. To fulfill this purpose, the average rainfall data of last 30 years of two stations (Suri, Dumka) have been taken and district resource maps and the Ph.D. thesis of Bhattacharyya (2014) have been consulted. For producing the data layers ( Figure  3g), soil erosivity has been calculated using the modified Fournier index (Eq. (6)) (1960). The parameter of water coverage has been considered in this study because of the presence of large reservoir river projects like Masanjore or Tilpara in the area. The data layer ( Figure 3m) has been prepared by categorizing the basin area into water logged and nonwater-logged areas and demarcating the reservoir area from the SOI toposheet. Variance inflation factor (VIF) statistics have been computed for this. V€ or€ osmarty and Dobos (2020) reported value of <5 does mean there is no co-linearity problem. The result shows that the VIF score is within permissible limits, therefore, no significant multi-co-linearity exists among selected factors.
where As ¼ specific catchment area in meters, r ¼ slope gradient in degrees where T B ¼ At satellite temperature, k ¼ wavelength of emitted radiance in meters, q ¼ h Ã c=r, h ¼ Planck's constant, r ¼ Boltzmann constant, and c ¼ velocity of light e ¼ Land surface emissivity.
where b IR ¼ Infra-Red band brightness value, b R ¼ Red band brightness value where Pi ¼ Precipitation of the month I, and P ¼ Mean annual Precipitation.

Modelling gully erosion susceptibility
3.4.1. Selection of training sites for PBM and SPBM For selecting the training sites for the gully erosion susceptibility cluster modelling, the Gully affected sites are considered '1', and non-affected sites are considered '0'. 3658 points from both the gully and non-gully erosion points, 70% (2560 points) of data have been used for modeling, and 30% (1098 points) were kept for validation. The training sites for gully affected and non-affected points have been recognized using SOI toposheet with 1:50,000 scale, high-resolution Google Earth imagery with 2.62 m. spatial resolution, and field investigation. Since the gully-affected sites were mostly concentrated in the upper catchment of the study area, both the training and testing sites have been taken from these regions. After delineating affected regions with the help of SOI toposheets, Google earth imagery and field investigation, training sites have been randomly extracted as point feature from those demarcated gully erosion prone areas. Here distance from gully was not considered. In this context, studies like S€ uzen and Doyuran (2004) suggested considering equal numbers of vulnerable and non-vulnerable points for better performance of the employed models. Following this principle, equal training sites were extracted from gully affected and non-gully affected areas. The broken ground and erosion prone surface areas have been kept aside while training sites extraction for non-gully affected areas.

Factor cluster modelling technique
After the Selection of the sample dataset, three ML classification algorithms, random forest (RF), AdaBoosting (ADB), and support vector machine (SVM), have been applied to the training dataset. During this training process, 1 has been assigned for gully affected training sites, and 0 has been assigned for non-gully affected training sites. In the sample point based modelling (SPBM), the output of each cluster was mapped with the help of interpolation technique using inverse distance weighting (IDW) in the ArcGIS environment. Whereas, in PBM, the classification algorithm has been run on the entire dataset, including all pixels of the study area. Although the ML techniques have run upon the basic regression techniques, the ML techniques have superiority over traditional statistical methods for delineating the database's complex behavior, especially the geospatial phenomena (Chowdhuri et al. 2021;Panahi et al. 2021). Therefore, we have employed the three ensemble machine learning techniques: Random Forest (RF), AdaBoosting (ADB), and Support Vector Machine (SVM) to explore the gully erosion susceptible zones for the present study area. Zhou et al. (2020), Pal and Paul (2021), and Pal and Debanshi (2021) have successfully employed these ML techniques to inspect and map various geospatial phenomena in their studies. A detailed methodological analysis of the SPBM and PBM has been discussed below.
In the PBM, the prediction of all the pixels of the basin area was in binary expression i.e. 1 for gully erosion susceptible and 0 for non-susceptible, and eventually, the factor cluster wise gully erosion prediction mapping has been done in binary expression. But in the case of SPBM, the interpolated map, the susceptibility value varied 0 to 1, where 0 indicated non-susceptible and 1 showed susceptible. So the maps of gully erosion prediction derived from each factor cluster have been converted to a binary map considering 0.5 as the threshold. In this way after producing the factor cluster wise gully erosion prediction maps for both PBM and SPBM approach the maps have been spatially integrated in order to yield gully erosion susceptibility. Since there were four maps containing 0 and 1 values, the values of gully erosion susceptibility maps varied from 0 to 4. In this map 0 represents those areas where all four clusters negated gully erosion and 4 represents those parts where all four clusters predicted gully erosion. Based on the agreement of the cluster on gully erosion prediction, the susceptibility has been categories into five classes, viz. extremely susceptible, highly susceptible, moderately susceptible, least susceptible and relatively safe. The susceptibility value 4 has been categorized as extremely susceptible followed by susceptibility value 3 which denotes prediction agreement of 3 clusters has been categorized as highly susceptible. The susceptibility values 2 and 1 denote the agreement of two and one clusters and have been categorized respectively as moderately susceptible and least susceptible. The 0 susceptibility denotes the agreement of no cluster regarding the prediction of gully erosion and therefore, has been labeled relatively safe.

Random Forest (RF)
In unsupervised ensemble learning, the random forest has proven to be a reliable and robust algorithm for classification and regression (Schonlau and Zou 2020;Meshram et al. 2021). The random forest algorithm uses the bootstrap aggregation technique while growing multiple decision trees (Syam and Kaul 2021). Therefore, each tree can predict and 'vote' independently at the time of supply of the input data for the corresponding classes a ð Þ, X Y rf ¼ majorityvote X y a ð Þ Ù È É Y 1 , where X y a ð Þ Ù the class's prediction of the bth random forest (Hartini et al. 2021). The voting capability of RF reduces the noise and improves the robustness and prediction capability of the input data (Zhou et al. 2020). RF reduces the dimensionality in the database by using a built-in feature selection system that can regulate multiple parameters without removing some of them (Tella et al. 2021). By estimating the increase in prediction error in the dataset, RF can compute the variable importance scores of each constituent tree and also the entire database (Pal and Paul 2020; Elavarasan and Vincent 2021). Studies like Chan and Paelinckx (2008) and Pal and Mather (2003) reported that RF (uses bagging) is more sensitive to noise than the algorithms based on boosting technique. Finally, RF uses proximity between pairs of variables to locate the outliers (Norouzi and Moghaddam 2020). To calculate the proximity between two samples, RF calculates the frequency of appearance of the sample at the same terminal node (Wangchuk and Bolch 2020). After the computation of proximity for each pair of causes and the formation of trees, the variables are normalized by dividing with numbers of trees (Nahkala et al. 2021). For the present study, a random forest algorithm from the Scikit-learn ML library from python 3.6 has been applied for PBM, and the Random Forest tool from WEKA has been used for SPBM.

AdaBoosting (ADB) classification
AdaBoosting, or adaptive boosting (ADB), is a tree-based ensemble algorithm that uses an iterative process to generate weight at each sample for both classification and regression (Mehmood and Asghar 2021;Pal and Paul 2021). ADB consecutively grows multiple interdependent classifiers and generates multiple classifiers to 'vote' them (Xu et al. 2021). ADB generates equal weights for each dataset. In the iteration process, the weights of each unclassified dataset have been increased, whereas the weight for correctly classified samples is gradually decreased . These weights indicate the classifier's accuracy and are a function of the overall weights of the correctly classified samples . The algorithm's accuracy depends on the given weight of correctly classified samples. Unlike RF, ADB is sensitive to noise in the training data; therefore, it is necessary to preprocess the training data to obtain optimum output (He et al. 2020). Python-based AdaBoost Classifier from ensemble learning tool from Scikit-learn ML library is used to classify factor clusters in this present study, whereas AdaBoostM1 tool in WEKA is used for SPBM classification.

Support vector machine (SVM)
SVM is a popular robust supervised classification technique initially developed for binary classification problems (Cortes and Vapnik 1995). In this classification, several classes can be determined as an arrangement of n-class or n(n À 1)/2or binary classification techniques by applying different selection schemes that lead to a final decision (Ul Din and Mak 2021). SVM generates an optimum hyperplane in a given training set that elucidates a complete separation of linearly different clusters in a hyper-surface (Talukdar et al. 2021). During the training, SVM searches for a hyper-plane that best separates the dataset into two classes (binary), such as 1 and À1. These instances can be identified as h 1 and h À1 or support vectors (Deiss et al. 2020). The hyperplane between h 1 and h À1 is a 0 or non-linear convex programming problem. The pairs of parameters of the separating hyperplane determination can be determined by solving the following optimization problem: In respect to 1Àn i Ày i w Á where penalty parameter (C) and margin of tolerance (n) needed to be tuned to accuracy and efficiency of the model, for this study, the SMOreg tool from WEKA has been used for SPBM, and SVC algorithm from Python has been used for PBM classification in Python. GridsearchCV tool from the Scikit-learn ML library has been used for hyperparameter selection.

Hyperparameter optimization
The selection of a proper set of hyperparameters is necessary to extract maximum performance and accuracy from the model and is considered an integral part of machine learning-based modelling (Kotthoff et al. 2019). Studies like Khalid et al. (2020) and Victoria et al. (2021) reported that high accuracy and performance could be achieved by tuning hyperparameters. There are various techniques to tune hyperparameters like randomized search, grid search, Bayes optimization, halving grid and random search, hyper Opt-Sklearn, kNN, or IBK. For this, in the present study, kNN or IBK techniques have been used to tune the hyperparameter in WEKA for SPBM, whereas grid search optimization and the K-fold cross-validation technique have been used to adjust the hyperparameter and validate the predicted dataset. Each cluster model has run a 10-fold iterative process for every 340 candidates totaling 2400 fits. In Python, the GridSearchCV tool from the SK-learn ML library is employed for this present work.
The final output layers of pixel based modelling (PBM) and sample point based modelling (SPBM) have been classified into five equidistant susceptible zones: relative safe, least susceptible, moderately susceptible, highly susceptible, and extremely susceptible classes. The areas emerged as low to very low prone to gully erosion, categorized as the least susceptible and relative safe zones. In contrast, areas with increasing gully erosion effects are moderately, highly, and extremely susceptible zones.

Accuracy assessment of the PBM and SPBM
To affirm the performance of the employed models, six statistical matrices, namely percentage of correctly classified data, receiver operating characteristics (ROC) curve derived area under curve (AUC), Precision, Sensitivity, F1-score, and Matthew's correlation coefficient (MCC), have been calculated based on 1098 sample points. The validation techniques have been calculated based on four matrices which are true positive (TP), false negative (FN), true negative (TN), and false-positive (FP). Based on these matrices, the specificity and sensitivity of the models have been calculated. The ROC shows the distinguishing capability of the model between classes and is constructed by plotting sensitivity (Y-axis) against 1-specificity (X-axis), where AUC denotes the degree of separability (Narkhede 2018). Sensitivity or recall has been calculated using TP and FN, and specificity has been calculated using TN and FP. In MCC, all four matrices have been used for better reliability and performance (Chicco et al. 2021). The MCC value ranges from À1 to 1, where À1 indicates the least accuracy in the prediction and vice-versa. The F1score has been calculated against the value of precision and sensitivity. Therefore, it gives a better result, even better in uneven class distribution (Chicco and Jurman 2020). The value of ROC-AUC, Precision, and F1-score, ranges between 0-1, where absolute 0 or near to 0 indicates less reliability on the model and 1 indicates maximum reliability on the employed models (Yao and Shepperd 2020). All four matrices have some superiority and flaws against another; therefore, a score of any single model cannot be taken as a good indicator for model validation. The calculation for the stated matrices is as follows (Eqs. (9)-(12)).
4. Results Figure 4 displays the comparative views of pixels and sample points based on ensemble ML models for gully erosion susceptibility in the upper Mayurakshi river basin. Each map has been classified into five gully erosion susceptible zones, and the area under each class is shown in Table 1. About 15-20% of the area has been detected as a very highly susceptible gully erosion zone as per pixel based ML models. On the other hand, this ranges from 1-26% in the case of sample point based ML models. From these statistics, it can be stated that sample point based modelling reflected a greater degree of prediction uncertainty. The gap between pixel and point based very highly susceptible area have been found very high in the ADB model and found low in the case of the RF model. Such inconsistency in prediction is caused for the question of the reliability of sample point based models. Geographically, the upper parts of the study region are high to very highly susceptible to gully erosion. Highly erodible coarse-textured laterite soil, relatively sloppy Chottanagpur plateau fringe, and highly dense 1st and 2nd orders streams are mainly responsible for greater gully erosion susceptibility (Pal 2016;Pal and Debanshi 2018; Debanshi and Pal 2020). Areal patches of very highly erosion susceptibility look relatively continuous in the case of point-based modelling, and this geolocation is deliberately different from pixel based ML model output. So, justification of the suitability of the model's application for such work is essential. Moreover, it is also necessary to justify the suitability of the pixel and point based ML models. For that accuracy, a test is highly crucial. Table 2 shows the accuracy and performance level of the applied models with the help of some statistical techniques like AUC-ROC, precision, sensitivity, F1 score, and MCC. In the case of sample point based ML models, AUC and MCC range from 0.53-0.74 and 0.48-0.69, whereas it ranges from 0.81-0.91 and 0.78-0.91, respectively, in the case of pixel based ML models. Sensitivity, precision, and F1 score also represent the same trend (Table  2). Based on these statistics, it could be inferred that pixel based gully erosion susceptibility models are superior to sample point based models. Obviously, in the case of sample point based models, the sample size may influence the accuracy level; however, at any sample size, it will not be as precise as the pixel based models. Among the tree applied ensemble ML models, RF was the best representative since the accuracy level, as per all the test statistics, was high in pixel and point-based models. AUC-ROC value shows a very good to excellent agreement between training and test data set in case of all the models applied for predicting pixel based susceptibility. However, RF has been recognized as the best suited. The agreement between the test and training dataset was unacceptable to fairly good in point based ML models. To assess the stability of the validation result, we have performed the validation technique with 11 different sample sites. However, the result did not vary significantly, and the best performance result has been mentioned here (Table 2).

Gully erosion susceptible models PBM and SPBM
Field data regarding extension of the gully, widening, and deepening of the gully has been collected from thirty-seven sites to validate the ML based gully erosion, susceptible models. Gully erosion rate have been correlated with the degree of gully erosion susceptibility. Correlation coefficients between the degree of susceptibility and head ward extension of the gully, gully widening and deepening were 0.65-0.77, 0.58-0.62, and 0.73-0.81, respectively, in case of pixel based models, and all these values are statistically significant at 0.05 to 0.1 level of significance. The highest correlation has been found in case of the RF model. In the case of sample point based ML models, correlation coefficient values range from 0. 12-0.35, 0.13-0.29, and 0.17-0.34, and the significance level is <0.05.  topographical factor cluster. The nature of susceptible areas has been determined based on the considered factors under each factor cluster. Point based models show continuous patches of susceptible areas, whereas its areal extension is highly fragmented. All the factor clusters based model has recognized comprehensive parts of the upper catchment area under gully erosion susceptibility. This proportion is susceptible high in the case of topographical factor cluster. Both pixel and point based factor cluster models have recognized upper parts of the study region is under soil erosion susceptibility; however, in terms of geolocation, these are not always coinciding. So, it is necessary to justify whether the pixel based model is superior to the sample point based model. Here one thing is clear if the factor cluster based model is validated with the testing data set, it may show relatively poor agreement in most cases. However, considering the relative difference in statistical values of the applied measures, it could be easy to infer the suitability of the models.  agreement between training and testing datasets. Here also RF model has appeared the best suited than other models.

Discussion
The study exhibited that about 20% area at the upper part of the study area is highly susceptible to gully erosion. The predictive ability of the RF model is the most credible since its agreement between training and the testing dataset is very high. The predictive ability of pixel based models is far better than sample point based models. Why is a higher degree of gully erosion susceptibility found in the upper part of the study region? Why predictive ability of pixel based models is higher than sample size based models. May sample size be an issue in improving the accuracy of the models?
The upper reach of every river basin is characterized by higher 1st and 2nd order stream frequency and density and is prone to gully erosion with different intensities (Das 2014). New gully formation and extension of the gully are very general mechanisms in such an area. But in a river basin where highly erodible coarse texture soil is available, such a rate is very high. In the Mayurakshi river basin, the lateritic composition of soil enriched with silica, aluminium, iron etc., is highly prone to chemical weathering and consequent erosion (Ghosh et al. 2015). A field study also revealed that gully formation and extension are very prominent even in the forested lateritic region. In the region where hard lateritic duricrust is there, the rate of gully erosion is relatively less. In the present case study, a few lateritic patches are where hard cap lateritic mass is armoured by the fragile laterite soil. Poor undergrowth vegetation in the Sal (Shorea robusta) forest available in a wider part of the upper catchment also insists intensive gully erosion. The growing rate of stone collection from the hilly region, settlement growth in the forested region, and illegal deforestation are all responsible for creating an ambient condition of highly intensive soil erosion. Some forested areas and their peripheral part are also prone to high gully erosion susceptibility. Higher drainage frequency and density in its upper part coupling with a relatively steeper slope and fragile laterite soil firmly control the gully erosion. Pal and Debanshi (2018) also condemned the fragility of laterite soil as a dominant factor for the high rate of gully erosion in this area. Pixel based modeling in Python software provided more accurate results than sample point based modeling in this study. There is no need for interpolation while every pixel has been considered for analysis. So, here total pixel is the total sample points. But in the case of sample point based modeling, unknown sites have been predicted using the interpolation method based on known values of the training sites. Interpolation often generalizes the trend of value gradient. It may not always represent reality properly (Abdulmanov et al. 2021). During interpolation, several issues are associated with making the map more precise. Sample point size for interpolation is the first and foremost issue that can dictate the precision of the interpolated map. The sample size is indirectly related to the degree of accuracy of the interpolated map (Simpson and Wu 2014). In this study, 3658 sample sites have been taken for analysis.
Ensemble machine learning models have been applied for predicting gully erosion susceptibility considering its robustness and high precision level. In the present study, all models have been found acceptable, but precision is high in the case of the RF model. In the RF algorithm, a built-in features selection system is used to reduce dimensionality without removing any data in the database. Therefore, there is very little chance of data loss or outlier intervention in the dataset, making this algorithm more reliable and errorfree than the other algorithm (Zhou et al. 2020). The RF (which uses bagging) is more sensitive to noise than the algorithms based on boosting technique; therefore, it generates more accurate predictions by controlling the noise in the dataset (Pal and Mather 2003;Chan and Paelinckx 2008). The multi-model approach has been adopted here to enhance the reliability of the comparison. Multiple models with the same aim cross-validate themselves and, therefore, are more reliable in the modeling approach . Moreover, such an approach helps justify the best representative model of all. Considering these advantages, in recent times so, many scholars like Kadavi et al. (2018), Mosavi et al. (2018), Hong et al. (2018), Liu et al. (2017), and Jamali (2019) have taken this line of thinking.
The identification of extreme gully erosion susceptible zones has its implication from the management point of view. These areas deserve special care for resisting gully erosion and soil loss. There may be some ways to check gully erosion in this region theoretically; however, check damming is an effective measure. For the last seven decades, this engineering mechanism has been highly relied on. Small scale check damming can break down the slope and store the regional soil. A good number of check dams have beeninstalled in the catchment of the Massanjore dam over the Mayurakshi river to check soil erosion and thereby lengthen the longevity of the dam (Pal 2016;Debanshi and Pal 2020). So, this could be a good way. Renovation of existing dams and installing a new dam in the deserving sites may check the rate of soil erosion. Further research is necessary to determine the deserving sites for new check dam sites. As far as the soil type is concerned, it is pretty challenging to arrest such erosion but check damming may slow down such rate. Since the natural vegetation of this region is tall with lesser canopy density, that cannot contribute much to protecting the soil from splash erosion. So vegetation area should be increased, emphasizing dwarf varieties of plants having greater canopy density, which will reduce fall velocity and provide good resistance to gully erosion.
The present study has compared the pixel and point based ensemble ML models for gully erosion susceptibility. This approach is very new, and the finding of this may encourage scholars to do with software environments like Python, where every pixel is a sample. It also proved that sample point-based modeling, which is done frequently, has considerable noise, and such output may mislead spatial planning. It is further to be mentioned that here only the output of particular sample size has been generated; if instead of this, the same modeling process was run using different sample sizes and accuracy was assessed, the result was more easily convincing. So the analysis of sample size effect on model output is a good scope of work.

Conclusion
The present study inferred that pixel based ensemble models' predictability is higher than sample point based models. The RF model shows the highest precision for predicting gully erosion susceptibility in pixel and point based modeling. Coarser laterite soil dominated fertile soil in the upper part of the study region is more susceptible than its counterparts. So, the study recommends using pixel based ensemble modeling to predict susceptibility instead of sample point-based modeling. The conventional practice of modeling with sample point data using interpolation technique generalizes the spatial character and thereby reduces the precision of the model. An increase in the undergrowth and check damming may be some fruitful measures to slow down the gully erosion rate. In such a case, implementing management strategies with limited resources require greater accuracy and a higher degree of spatial differentiation for prioritizing the area of action. Since the study demonstrated the gully erosion susceptible region, this information will be useful to the decision-makers for protecting gully erosion and soil resources.