Flash flood susceptibility mapping using LiDAR-Derived DEM and machine learning algorithms: Ljuboviđa case study, Serbia

2.1 Study area

The study area covers approximately 537.55 km² between the latitudes 44°08′ to 44°21′N and the longitudes 19°20′ to 19°45′E and is located in the western part of the Republic of Serbia (Figure 1). The study area includes parts of Kolubara District and Mačva District. The elevation of the study area ranges from 162 to 1250 m a.m.s.l, with an average altitude of 584.3 m, and with a standard deviation of 219.7 m. There are plains (Kašica, Donja i Gornja Ljuboviđa, Gračanica) and mountains (Bobija, Medvednik, Ranjenica, Nemić) in the study area, making it diverse in terms of the complexity of the topography.

Figure 1

Location of the study area (parts of Kolubara and Mačva District).

In the past, the test area has experienced larger, more intense flash floods. The largest ones occurred on November 18, 1979, May 30, 1985, June 21, 2001, May 14, 2014, and June 22, 2020 [25]. The total annual rainfall in the different parts of the study region ranges from 913.34 to 1205.33 mm. According to records obtained from the Serbian Meteorological Service [26], the maximum rainfall occurs between March and June.

2.2 Data used

A crucial stage in the mapping of flash flood susceptibility is the creation of the inventory map. A flash flood inventory map shows one or more instances of previous floods in a particular location. For flood susceptibility mapping, it is crucial to have an inventory map that shows past flash flood occurrences. Making accurate predictions about future flash flood events in a location requires evaluating the historical records of flash floods. In this study, we used aerial images, field surveys, and high-resolution WorldView 2 imagery to map the flooded area. The study’s flash flood inventory data, which includes 105 recorded flash floods, spans the years 2012–2024.

A further important subject that has been studied by numerous scientists is the flash flood conditioning factors. In order to assess the risk of flash floods, this study used a number of factors, including elevation, slope, aspect, plan curvature, topographic wetness index (TWI), stream power index (SPI), normalized difference vegetation index (NDVI), rainfall, land cover/use (LULC), lithology, and the distance to the stream and roads. Details about the origins of the flash flood conditioning factors are provided in Table 1.

In order to identify morphometric characteristics and indices, the research of a flash flood basin often starts with topographic maps and DEM. Because of the significant impact that the basin’s morphometric features have on the precipitation–runoff connection within the basin, their interpretation can help explain how susceptible a particular basin is to flash floods. Consequently, the geomorphology of the basin, which regulates the hydrological response during periods of intense precipitation, is quantitatively reflected by morphometric parameters.

A Leica ALS80 LiDAR system was used to gather LiDAR data from May 24 to July 24, 2022, at a height of 1,500 m above the ground. The data were obtained from the survey using the Piper Seneca V aircraft and the Leica ALS80 laser scanner (Figure 2). For a single laser pulse, the laser scanner was set up to record up to four returns. The point spacing is 1.0 m on average [27]. For this study, LiDAR data were used, with a vertical accuracy of 0.15 m and a horizontal accuracy of 0.22 m. Delivered in binary LAS file format, all LiDAR data are divided into ground and non-ground points. The classification was done in the TerraSolid software, so for the generation of a DEM with a resolution of 1 m, only the ground class was used, and small irregular shapes were manually cleaned (Figure 3a).

Figure 2

Piper Seneca V aircraft and the Leica ALS80 laser scanner.

Figure 3

Conditioning factors related to flash floods: (a) the elevation; (b) the slope; (c) the aspect; (d) plan curvature; (e) the TWI; (f) the SPI; (g) the rainfall; (h) lithology; (i) NDVI; (j) distance to streams; (k) land use/cover; and (l) distance to roads.

The maximum magnitude of the slope change for each raster cell and its eight adjacent cells is represented by the slope angle in the DEM. The accelerated rate of particle movement and transport is impacted by the speed, kinetic energy, and tangential stresses of surface runoff, which all rise with increasing slope. In contrast, the dynamics of the process are considerably softer on gentler slopes. Long, steep slopes have a high potential for surface runoff generation, which subsequently causes gullies and furrows to form. With a range of values from 0 to 70.72°, the average terrain slope in the studied area is 16.09° (Figure 3b). Strong erosion, together with heavy soil mass washing and movement, are characteristics of this type of slope.

Aspect is likewise regarded as a secondary component that affects soil humidity and floodwater flow direction. Additionally, the slope’s curvature affects water flow because regions with zero curvature are often more vulnerable to floods than regions with positive or negative curvature values. DEM was also used to create an aspect map with nine classes – flat, north, northeast, east, southeast, south, southwest, west, and northwest – that show the direction of the terrain (Figure 3c).

The term “terrain curvature” refers to the shape of a slope and its influence on the genesis and development of erosion processes. For analysis purposes, terrain curvature is calculated using DEM. It is expressed in rad/m, but due to small values and for easier presentation of the results, it is multiplied by 100. In the horizontal plane, plan curvature represents the curvature of the contour line and is perpendicular to the direction of the greatest slope (Figure 3d). Water movement’s convergence and divergence are depicted by plan curvature. Whereas convergent slopes create areas of water and sediment buildup, divergent slopes provide runoff that influences how intense the washing and dredging operations are. Positive and negative plan curvature values are also possible. Laterally convex surfaces have positive values, whereas laterally concave surfaces have negative values. It was established that convex slopes predominate over concave ones in the studied area. These surfaces show that erosion processes (scouring, furrowing, and dredging) have begun and that surface runoff is occurring quickly.

Another crucial hydrological element of flash flood susceptibility analysis is the TWI, which captures the patterns of soil moisture associated with floodplains and measures the likelihood that water will collect at any point within a watershed due to gravity. TWI describes the effect of topography on the location and size of saturated areas of runoff generation. It is defined as [28]:

where CA is the catchment area and β is the slope angle in degrees.

Water pooling and downstream movement parameters can be thoroughly assessed thanks to the TWI, which incorporates data on flow dynamics and gravitational influence across various watershed zones (Figure 3e). In mapping flash flood susceptibility, regions with higher SPI values are more vulnerable to flash floods because of the quick erosion and channel creation that might result from strong water flow. SPI is a measure of the erosive power of flowing water based on the assumption that discharge is proportional to specific catchment area. SPI was calculated based on the formula given by Moore et al. [29]

where CA is the specific catchment’s area and β is the local slope gradient measured in degrees. Because of the extreme energy of the water flow, high SPI locations are therefore recognized as having a greater danger of abrupt and catastrophic flooding (Figure 3f).

Heavy rainfall is the primary cause of the interrelated and nearly simultaneous processes of runoff, washing, and soil erosion. Their byproducts, a vast quantity of water and sediment, then join the hydrographic network and continue to move as a two-phase fluid. The water flow and sediment transport in torrential flood waves comprise the largest portion of the overall annual water runoff and sediment transport, which is influenced by the water power in torrential processes caused by heavy rainfall. Then, soil erosion is the most severe process that degrades land resources by removing the soil’s top layer. The total annual rainfall in the different parts of the study region ranges from 913.34 to 1205.33 mm (Figure 3g).

This study also took into account a few environmental parameters, like LULC, NDVI, lithology, and distance to streams and roads. Geology is generally accepted to control the drainage pattern and the method by which rainfall is runoff. Flash flood probability is significantly influenced by lithology. The amount of water infiltration is regulated by this component, which has a significant effect on the creation of surface runoff (Figure 3h).

Additionally, NDVI is calculated by dividing the difference between the red (RED) and near-infrared reflectances by their sum, and it shows the health of the vegetation [30]. The NDVI scale goes from +1 to −1, with +1 denoting dense green-leafy vegetation and −1 denoting aquatic bodies. Vegetation has a detrimental effect on flooding, as was previously established. NDVI was calculated from the Sentinel 2 images collected on June 20, 2022. Cloud cover was 3%. This means that only a small portion of the sky, 3% to be exact, is obscured by clouds. The ArcGIS program was utilized to classify the Sentinel 2 data according to the NDVI (Figure 3i).

Flood susceptibility is significantly influenced by LULC. Different land cover types respond to floods in different ways. Densely vegetated areas, for instance, are less likely to flood because they impede water flow, but non-vegetated places, like bare lands or urban areas with impervious surfaces, accelerate rainwater runoff. Ten LULC classes are included in the CORINE land cover 2018 dataset, which we used for this analysis (Figure 3k and with a description in Table 2).

Finally, using a digital topographic map and the Euclidean distance tool in ArcGIS software, distance to streams and roads data layers were created (Figure 3j and l). The scale of the used digital topographic maps is 1:25,000 (DTM25). The flash flood conditioning factors are shown in Figure 3 and Table 3, and in order to preserve uniformity pixel-by-pixel, all of the factors, including the flood map, were resampled to 10 m.

2.3 Methods

A specific methodology for assessing flash flood susceptibility was developed, considering geographic and temporal constraints. A generalized flowchart is provided in Figure 4, and the following workflows were adopted to achieve the study’s objectives:

1. A detailed flood inventory was prepared using WorldView 2 images, aerial images, and field investigations, focusing on historical flood locations.

2. Based on prior relevant research and the geo-environmental characteristics of the study area, 12 predictor variables associated with flood susceptibility were identified and subsequently prepared.

3. Using the well-known scikit-learn (train_test_split) Python module, the flood inventory dataset was split into 70 and 30% for training and validation, respectively.

4. The multicollinearity of flash flood predictor variables was tested to prevent bias in the flood mapping method, and a popular grid search technique was adopted to adjust the hyperparameters.

5. The current flood prediction employed three modern machine learning techniques: SVM, CART, and RF. However, only the most effective ML: RF was utilized for the mapping of future flood predictions.

6. The accuracy of the model’s predictions was evaluated using statistical measures such area under the curve (AUC).

Figure 4

Flowchart for analysis of flash flood susceptibility based on machine learning models.

2.4 Multicollinearity analysis

Flash flood susceptibility mapping has been predicted using a variety of models based on different geoenvironmental parameters. Finding the two or more elements that have a high correlation with one another is essential among all of these conditioning components. A larger correlation between two or more conditioning elements than the model makes the evaluation less valid and automatically lowers the output result’s accuracy. Thus, the identification of these strongly correlated components is aided by multicollinearity analysis. Multicollinearity is essentially a statistical analysis of more than two variables that shows a linear relationship between them. For better outcome prediction, strong multicollinearity conditioning components must be eliminated from the models. Multicollinearity has typically been analyzed using the variance inflation factor (VIF) and tolerance (TOLERANCE) parameters. The following formula can be used to determine a multi-collinearity analysis’s TOLERANCE and VIF [31]:

where R j 2 represents the regression value of j on other variables in a dataset.

2.5 Machine learning method used in flash flood susceptibility modeling

SVM is a supervised machine learning algorithm, and it is used in classification and regression modeling. SVM is utilized in regression and classification modeling. The model, which is founded on the optimization principle, attempts to divide the several classes by fitting a hyperplane on the training data set [32]. As far away from the nearest data points from each class as feasible is how the hyperplane is oriented. The term “support vectors” refers to these nearest positions [33]. An alternative term for a hyperplane is a decision boundary.

The RF is one supervised machine learning algorithm used in classification and regression. The method comprises a set of tree-structured classifiers. Each tree assigns a unit vote to the input vector (x) based on the most frequent class. By dividing each node into a random subset of input characteristics or predictive variables rather than the best variables, a RF classification lowers the generalization error. Additionally, an RF uses bootstrap aggregation or bagging to force the trees to grow from various subsets of training data to boost the trees’ diversity [34].

One well-known and popular decision tree algorithm is CART. The CART model is applied to problems involving both classification and regression. The binary tree is created by this paradigm by dividing the input into binary components [22]. The root node of this algorithm initially contains all of the data. Next, among the modeling traits in question, those that best distinguish the target class are chosen and added to the root node.

In this study, a set of candidate hyperparameters was chosen using grid search. Every potential combination of the chosen values is represented by a model. Performance is then assessed by testing each combination with 12-fold cross-validation. Table 4 displays the hyperparameters derived from the grid search.

2.6 Methods of validation and accuracy assessment

To measure the predictive outcome, flash flood susceptibility maps must be validated and their accuracy evaluated. Therefore, the three machine learning output findings were evaluated in this study using a variety of statistical indices and the area under the receiver operating characteristic (AUC) curves. Sensitivity (SST), specificity, positive predictive values, and negative predictive values were employed in statistical indices. Every machine learning model produced a better outcome if the value of these statistical indices was higher, and vice versa. The following formulas can be used to determine the four statistical indices that were employed in this investigation [35,36]:

where TP represents true positive, TN represents true negative, FP represents false positive, and FN represents false negative.

3.1 Analysis of multicollinearity

VIF measures the degree to which multicollinearity increases the variance of an estimated regression coefficient. High VIF values can lead to unstable coefficient estimations, making it challenging to determine the precise effects of each independent variable on the dependent variable. To improve the prediction of results, strong multicollinearity conditioning components (VIF values more than 5) must be eliminated from the models. Taking into account the VIF and TOLERANCE limit, the multicollinearity test of 12 flash flood causative components was conducted for this analysis (Table 5). Slope (2.64) and distance to road (0.87) were linked to the highest and lowest VIF values, which ranged from roughly 0.87 to 2.64. The range for TOLERANCE was from 0.29 to 0.98. Rainfall (0.29) and aspect (0.98) were linked to the greatest and lowest TOLERANCE limit readings. Multicollinearity was not an issue with the factors chosen to estimate the risk of flash floods. Thus, the flash flood vulnerability in the current study area was modeled using all 12 explanatory variables [37,38].

3.2 Flash flood susceptibility modeling

The percentages of each model vary, according to the analysis of the models used in this investigation. In order to detect and forecast flood-prone areas, these variances rely on the level of flood susceptibility in the region [39]. On raster maps, the flood susceptibility maps produced by each model showed values between 0.0 and 1.0. In order to preserve uniformity pixel-by-pixel, all of the conditioning factors, including the final flash flood maps, were resampled to a resolution of 10 m. Every pixel was given a unique value that represented how vulnerable it was to flooding. In particular, more flooding susceptibility is indicated by a higher number on the scale and lower flooding susceptibility is indicated by a lower value [40]. Five susceptibility classifications were created from the resulting models (Table 6 and Figure 5). For the RF model, these categories display the following percentages: 15.49, 16.04, 15.67, 23.10, and 29.70% (Figure 5c). Similar estimates were made for the SVM models very low, low, moderate, high, and extremely high percentage classes, which were 18.54, 16.74, 19.36, 21.24, and 24.12%, respectively (Figure 5a). Furthermore, the CART models’ flood-predicting percentage classes were, in order, 16.45, 14.24, 22.01, 20.78, and 26.52% (Figure 5b). In Figure 5, it can be seen that over 95% of the populated places (given in yellow) of the test area are located in high and extremely high zones susceptible to flash floods.

Figure 5

Flash flood susceptibility maps generated using machine learning models: (a) SVM, (b) CART, and (c) RF.

The Ljuboviđa basins’ flood-prone regions are primarily found in the west and north-west regions of the study area, according to the results of the percentage susceptibility classes that each model projected. These regions trace the east-west main flow of the Ljuboviđa River. A steep slope between high and low places, which makes it easier for water to concentrate during periods of severe rainfall, is what defines flood zones. Additionally, this slope makes it easier for water to flow through the bank undercutting [39,41]. Another factor contributing to floods in these locations is the absence of vegetative cover. However, the study area’s west and center sections, as well as the non-floodable zones, are distinguished by extremely dense vegetation, which prevents extended flooding in these areas.

3.3 Evaluating various ML algorithms’ performance

A crucial step in creating and refining machine learning models for various phenomena is accuracy assessment. The majority of recent machine learning research uses this curve and AUC to assess the classifiers’ performance [42,43].

An asymptotic 95% confidence interval for each of the three models used is displayed by ROC curve analysis. The CART model scored an average AUC value of 0.831 (with a confidence interval between 0.724 and 0.938), whereas the RF model was found to be the best model in this investigation with an average AUC value of 0.854 (with confidence values ranging between 0.761 and 0.947). Furthermore, the SVM model exhibited an average AUC of 0.802 (with a confidence CI of between 0.703 and 0.901), which is shown in Figure 6 and Table 7. The standard errors for the RF, CART, and SVM models were 0.047, 0.052, and 0.050, respectively (Table 7). The RF model was thus selected as the most effective model for predicting regions that are vulnerable to flooding.

Figure 6

AUC of models run on the validation dataset.