Machine Learning-Based and 3D Kinematic Models for Rockfall Hazard Assessment Using LiDAR Data and GIS

: Rockfall is one of the most hazardous phenomena in mountainous and hilly regions with high and steep terrain. Such incidents can cause massive damage to people, properties, and infrastructure. Therefore, proper rockfall hazard assessment methods are required to save lives and provide a guide for the development of an area. The aim of this research is to develop a method for rockfall hazard assessment at two di ﬀ erent scales (regional and local). A high-resolution airborne laser scanning (ALS) technique was utilized to derive an accurate digital terrain model (DTM); next, a terrestrial laser scanner (TLS) was used to capture the topography of the two most critical areas within the study area. A staking machine-learning model based on di ﬀ erent classiﬁers, namely logistic regression (LR), random forest (RF), artiﬁcial neural network (ANN), support vector machine (SVM), and k-nearest neighbor (KNN), was optimized and employed to determine rockfall probability by utilizing various rockfall conditioning factors. A developed 3D rockfall kinematic model was used to obtain rockfall trajectories, velocity, frequency, bouncing height, kinetic energy, and impact location. Next, a spatial model combined with a fuzzy analytical hierarchy process (fuzzy-AHP) integrated in the Geographic Information System (GIS) was developed to assess rockfall hazard in two di ﬀ erent areas in Ipoh, Malaysia. Additionally, mitigation processes were suggested and assessed to provide a comprehensive information for urban planning management. The results show that, the stacking random forest–k-nearest neighbor (RF-KNN) model is the best hybrid model compared to other tested models with an accuracy of 89%, 86%, and 87% based on training, validation, and cross-validation datasets, respectively. The three-dimension rockfall kinematic model was calibrated with an accuracy of 93% and 95% for the two study areas and subsequently the rockfall trajectories and their characteristics were derived. The assessment of the suggested mitigation processes proves that the proposed methods can reduce or eliminate rockfall hazard in these areas. According to the results, the proposed method can be generalized and applied in other geographical places to provide decision-makers with a comprehensive rockfall hazard assessment. study was conducted to select the optimized hybrid machine learning model (stacking) based on various machine learning classiﬁers KNN. Grid search The rockfall characteristics (kinetic energy, frequency, bounce height, and impact locations) derived through raster modelling were used to produce the rockfall hazard maps. In order to produce rockfall hazard maps, a developed spatial model was employed taking into account the e ﬀ ect of each factor based on the fuzzy-AHP outcomes. After collecting the opinions of the experts, the geometric mean was applied and the ﬁnal pairwise matrix was prepared. The analysis of the pairwise matrix based on fuzzy-AHP reveals that the value consistency ratio (CR) is 0.03 while the values of consistency index (CI) and randomness index (RI) are 0.019357 and 0.58, respectively. The result of fuzzy-AHP shows that kinetic energy has a high inﬂuence on rockfall hazard in comparison with frequency, impact location, and bounce height, respectively. The derived weights of rockfall energy, frequency, impact location, and bounce height were 0.48, 0.3, 0.1, and 0.12, respectively. The generated maps of rockfall hazard are illustrated in Figure 10.


Introduction
Rockfall is a type of landslide that can be defined as blocks detaching from a cliff face and heading rapidly downslope with different motion modes encompassing flying/falling, bouncing, rolling, and sliding [1,2]. These blocks vary in size ranging from decimeters (boulder fragments) to meters (big rocks) [3]. Such a phenomenon frequently occurs in mountainous and hilly areas with high and steep terrain such as Malaysia which can cause damage to human or material [4]. Rockfall varies temporally and spatially, hence it is challenge to predict such events. Consequently, a proper rockfall hazard assessment is required to reduce or eliminate damages and save lives and properties through providing a comprehensive information of rockfall probability to the developers of urban areas.
The use of a digital terrain model (DTM) is a fundamental step in each rockfall hazard assessment. In recent decades, Light Detection and Ranging (LiDAR) techniques, both terrestrial laser scanner (TLS) and airborne laser scanner (ALS), have been increasingly employed in geo-hazard studies to create an accurate DTM [5][6][7]. This is because of its ability to produce high-resolution data in both horizontal and vertical directions. However, LiDAR raw data include both ground and non-ground points. Therefore, data pre-processing, filtering and interpolation processes are required to remove the outliers (due to the presence of LiDAR system error and ray multipath) and the non-ground points (man-made and natural features) [2]. DTM is the source of most significant rockfall conditioning factors such as altitude, slope, aspect, curvature, etc. [3]. Nevertheless, other relevant conditioning factors such as hydrological, geological, anthropogenic, and environmental factors are also significant for rockfall probability mapping [3]. Remote sensing (RS)-based data and Geographic Information System (GIS) tools, together with machine-learning approaches, have been popularly used in numerous studies. In recent studies, soft computing, heuristic, statistic, and deterministic models have been used to evaluate rockfall in terms of source identification [8][9][10], kinematic modeling [11][12][13], hazard assessment [1,14,15], and risk assessment [16][17][18]. The numerical modeling predicts either the blocks trajectories using Newtonian mechanics or the rockfall runout zone based on empirical measurements [2]. The trajectories modeling represents the rockfall physics besides kinetic energy and bounce height of a moving block, which are highly significant in the design of mitigation processes [3]. Several models have been developed for rockfall trajectory simulation. These models can be classified into three groups based on the properties of their simulation and terrain. First, the rockfall trajectories simulation can be performed utilizing either a probabilistic or deterministic approaches. Second, the rockfall kinematic can be simulated via a rigid body or a lumped mass method. Third, the topography of a slope can be represented by either a 3D terrain or a 2D profile [3].
In comparison with analytic and expert opinion-based methods, machine-learning approaches are considered more effective for landslide spatial prediction, including rockfall [19]. The essential conception of these approaches is the utilization of machine-learning classifiers for assessing rockfall probability through analyzing the spatial relationship between previous landslides incidents (i.e., inventory) and the conditioning factors [9]. Various machine-learning classifiers have been developed and employed for generating landslide probability maps in different areas worldwide. This includes logistic regression (LR) [20,21], random forest (RF) [22,23], artificial neural network (ANN) [24,25], support vector machine (SVM) [26,27], and k-nearest neighbor (KNN) [28].
Although numerous studies have been performed to assess rockfall hazard, the use of machine-learning classifiers is still very rare in this field, especially hybrid models. In addition, most of rockfall hazard studies have been conducted either for local scale (small area) or regional scale (large area) but performing a study in two different scales at the same time has not been done. Such a study produces an inclusive evaluation of rockfall in terms of sources, trajectories, their characteristics, leading up to hazard assessment. On the other hand, rockfall trajectories, velocity, frequency, kinetic energy, bouncing height, and impact locations were not considered in most of rockfall hazard assessment studies. Therefore, the current research optimizes a hybrid machine learning model based on various classifiers (LR, RF, ANN, SVM, KNN) for producing a rockfall probability map and thus identifying rockfall sources of the whole study area (regional scale), and employed a developed 3D rockfall kinematic model to assess rockfall trajectories and their characteristics for the two small areas within the study area (local scale). It also develops a spatial model based on fuzzy analytical hierarchy process (fuzzy-AHP) integrated into GIS to generate the final rockfall hazard maps. In addition, mitigation processes were suggested and assessed in these areas.

Study Area
Ipoh is one of the most important cities in Malaysia as it contributes to the economic and tourism sectors [29]. Due to the geological setting of Ipoh and the existence of active faults, rockfall incidents frequently occur in this area. Such events are mostly triggered by heavy rainfall throughout the year in addition to man-made slope cutting. Limestone hills are formed by a natural dissolution of carbonate via water that produce unique features which represent an amazing natural beauty to the landscape. Nevertheless, due to the existence of fractures and extensive joints within these hills, they can also cause threats to people, properties, and transportation corridors. Limestone hills with sub-vertical to vertical slopes and various morphological features such as caves are prominent in Ipoh. Granite and marble are also existing within the study areas. The temperature in Ipoh ranges from 23 to 35 • C with a high humidity (~84.1%) all around the year whereas the annual rainfall is about 2190 mm. The area is situated between the northwest corner (4 • 41 00"N, 101 • 03 00"E) and the southeast corner (4 • 30 30"N, 101 • 10 00"E) ( Figure 1a). Ipoh was chosen as the study area for the application of the proposed hybrid machine-learning model to determine the rockfall probability and thus the source areas for the entire study area (regional scale). After identifying rockfall sources, the most hazardous areas were chosen (Lang [ Figure 1b] and Rapat [ Figure 1c]) to perform the 3D rockfall kinematic modeling and produce rockfall hazard maps (local scale). These areas contain some temples which attract many visitors and tourists every day.

Materials and Methods
This section presents the proposed methodology in terms of input data, rockfall source identification, trajectories and their characteristics, hazard, and mitigation processes. The optimization of the proposed hybrid model based on various machine-learning classifiers and the integration of fuzzy-AHP with the developed spatial model for rockfall hazard are also illustrated here.

Data
The main datasets are LiDAR, rockfall inventory, and various spatial layers. An airborne laser scanner (ALS) device was employed to capture point clouds of millions of vector points (x,y, and z) for the study area and surroundings with flight height of 1000 m and a frequency rate of 25,000 Hz in 2015. The ALS produced a point cloud with point density of 9 points per square meter with rootmean-square-error of 0.15 m vertically and 0.30 m horizontally. The captured raw dataset was preprocessed to eliminate noises, outliers, thus generating a terrain model that is utilized to derive the conditioning factors of rockfall. The raw dataset of LiDAR includes ground and non-ground points, and therefore a filtering algorithms was applied to remove the non-ground points and deriving an

Materials and Methods
This section presents the proposed methodology in terms of input data, rockfall source identification, trajectories and their characteristics, hazard, and mitigation processes. The optimization of the proposed hybrid model based on various machine-learning classifiers and the integration of fuzzy-AHP with the developed spatial model for rockfall hazard are also illustrated here.

Data
The main datasets are LiDAR, rockfall inventory, and various spatial layers. An airborne laser scanner (ALS) device was employed to capture point clouds of millions of vector points (x, y, and z) for the study area and surroundings with flight height of 1000 m and a frequency rate of 25,000 Hz in 2015. The ALS produced a point cloud with point density of 9 points per square meter with root-mean-square-error of 0.15 m vertically and 0.30 m horizontally. The captured raw dataset was pre-processed to eliminate noises, outliers, thus generating a terrain model that is utilized to derive the conditioning factors of rockfall. The raw dataset of LiDAR includes ground and non-ground points, and therefore a filtering algorithms was applied to remove the non-ground points and deriving an accurate DTM. The rest point clouds were transformed into a DTM via an interpolation technique. In this research, a simple morphological filtering (SMRF) process [30], and kriging interpolation methods were employed to generate the DTM.A terrestrial laser scanner (TLS) with an accuracy of 3 mm was utilized to capture the slope topography from easily reachable locations in the two selected areas (Lang and Rapat) ( Figure 2). Several scans were conducted from different positions with overlapping to cover the study area and to avoid shadow areas. The TLS generate multi point clouds of millions of vector points with an accuracy of 3 cm. The main goals of the TLS observation were to generate a 3D terrain model of entire slopes (for rockfall kinematic modelling) and to characterize the rocks' geomechanical properties (the shape and volume of falling boulders). The derived point clouds were registered via using a geodetic Global Navigation Satellite System (GNSS) device ( Figure 2) to determine the locations of a number of scanning positions and tie points. The rockfall inventory data was gathered from various sources encompassing of field observations (Figure 2), historical records, and high-resolution aerial photos (0.1 m). The in-situ measurements were performed utilizing a precise GNSS device that receives real-time corrections from the Survey Department of Malaysia. Ninety-five samples of rockfall, with their related features (Figure 1a), and 95 samples of non-rockfall were prepared for probability assessment. The non-rockfall samples were selected randomly based on land use, geology, and a slope threshold. The inventory data was split into two groups: training (70%) and testing (30) [31]. The training dataset was utilized to train the models and the remaining dataset was utilized for model validation and accuracy assessment. Several GIS datasets such as land-use, geological, rainfall, vegetation, roads, faults, and streams were also included to prepare Remote Sens. 2020, 12, 1755 5 of 23 rockfall conditioning factors. These factors were obtained from different official agencies and some of them were registered and transformed into an identical coordinate system.
Ninety-five samples of rockfall, with their related features (Figure 1a), and 95 samples of non-rockfall were prepared for probability assessment. The non-rockfall samples were selected randomly based on land use, geology, and a slope threshold. The inventory data was split into two groups: training (70%) and testing (30) [31]. The training dataset was utilized to train the models and the remaining dataset was utilized for model validation and accuracy assessment. Several GIS datasets such as landuse, geological, rainfall, vegetation, roads, faults, and streams were also included to prepare rockfall conditioning factors. These factors were obtained from different official agencies and some of them were registered and transformed into an identical coordinate system.

Identification of Rockfall Source Areas
The LiDAR dataset includes both ground (bare earth) and non-ground (man-made and natural) points. Therefore, after removing outliers and noises, the filtering process was applied to remove

Identification of Rockfall Source Areas
The LiDAR dataset includes both ground (bare earth) and non-ground (man-made and natural) points. Therefore, after removing outliers and noises, the filtering process was applied to remove undesired features and then generated the final DTM using an interpolation algorithm. In this research, a SMRF process [30], and kriging interpolation methods were employed to generate the DTM. Consequently, various rockfall conditioning factors were derived including morphological factors such as altitude, slope, aspect, and curvature. Hydrological factors such as stream power index (SPI), sediment transport index (STI), topographic roughness index (TRI) and topographic wetness index (TWI) were also extracted from the generated DTM. A total number of 15 rockfall conditioning factors were prepared for optimizing the hybrid machine learning models. The rockfall conditioning factors comprise morphological and hydrological factors, land use, geology, rainfall, vegetation density, and distance to road, faults, and stream. Next, georeferencing process was applied to convert the datasets that were obtained from different sources into an identical data format while the inventory dataset was subjected to the removal of missing values. The GIS database, inventories, and LiDAR dataset were finally registered into the same coordinate system (WGS84). Since the sampling points of rockfall inventory are subjected to study relationships between various conditioning factors, the variance inflated factor (VIF) method was employed to evaluate the multicollinearity of the factors [32]. The inventory dataset was divided into two subsets one to train the model (70% of the inventory samples) and the other to validate the model (30% of the inventory samples). In addition, the inventory was used to apply a cross validation method (5-fold) to test the accuracy of the proposed hybrid models. A detailed comparative study was conducted to select the optimized hybrid machine learning model (stacking) based on various machine learning classifiers such as LR, RF, ANN, SVM, and KNN. Grid search method was applied to select the optimized hyper-parameters of each machine-learning classifier. The selection of the best hybrid model was based on the standard accuracy metrics [23] such as overall accuracy, receiver operating characteristics (ROC) curve, and cross-validation area under the curve (CV-AUC) [23]. Rockfall probability was then produced using the best hybrid model. These processes were implemented through Python language. To identify rockfall source areas, the probability map was intersected with the slope raster which was reclassified based on a slope threshold derived from the inventory dataset. Figure 3 shows the workflow of the process determining the rockfall probability and source areas.

3D Rockfall Kinematic Modeling
After identifying the rockfall source areas for the whole study area (regional scale), the most hazardous areas within the analysis area were selected to evaluate rockfall kinematics in terms of trajectories, velocity, frequency, bouncing height, kinetic energy, and impact location. These areas are Lang and Rapat that attract many tourists daily. Rockfall probable trajectories were modeled from all probable source samples utilizing a developed 3D probabilistic model to assess their spatial distribution and characteristics. The most significant dataset for rockfall modeling are geometric parameters (topography, source identification, outcropping materials, rock shape, and volume), and mechanical parameters (friction angle and coefficient of restitution). In this research, the rock shape and volume were determined based on field measurement and TLS dataset. Three-dimensional back analysis was conducted to calibrate the mechanical parameters by determining the best agreement between the field observation and rocks distribution obtained from the modeling. An initial velocity of 0.5 m/s was assigned to the falling block. The 3D rockfall modeling was performed using the RocPro3D [33] considering rock size and shape via employing the rigid body approach. In order to mimic the uncertainty associated with an initial falling direction, 10 blocks with different directions (3° interval) were thrown in each seeder points that distributed along the source line with an equal

3D Rockfall Kinematic Modeling
After identifying the rockfall source areas for the whole study area (regional scale), the most hazardous areas within the analysis area were selected to evaluate rockfall kinematics in terms of trajectories, velocity, frequency, bouncing height, kinetic energy, and impact location. These areas are Lang and Rapat that attract many tourists daily. Rockfall probable trajectories were modeled from all probable source samples utilizing a developed 3D probabilistic model to assess their spatial distribution and characteristics. The most significant dataset for rockfall modeling are geometric parameters (topography, source identification, outcropping materials, rock shape, and volume), and mechanical parameters (friction angle and coefficient of restitution). In this research, the rock shape and volume were determined based on field measurement and TLS dataset. Three-dimensional back analysis was conducted to calibrate the mechanical parameters by determining the best agreement between the field observation and rocks distribution obtained from the modeling. An initial velocity of 0.5 m/s was assigned to the falling block. The 3D rockfall modeling was performed using the RocPro3D [33] considering rock size and shape via employing the rigid body approach. In order to mimic the uncertainty associated with an initial falling direction, 10 blocks with different directions (3 • interval) were thrown in each seeder points that distributed along the source line with an equal interval (1 m).
Raster modeling was performed to obtain the rockfall frequency, bouncing height, kinetic energy, and impact locations. These characteristics were determined according to the outcomes of the 3D rockfall modeling (trajectory and velocity), and the DTM grid cells. The map of trajectories' frequency corresponds to a grid mapped onto the DTM, for which each cell gives the ratio of trajectories number that cross it in comparison with the total number or rockfall trajectories. The map of impacts was calculated for which each cell gives the number of impacts (dissipative impacts) that it contains. On the other hand, rockfall bouncing height was calculated for each cell and the highest height was assigned to a particular cell. The derived rockfall mass and velocity were used to calculate the kinetic energy and the highest value was assigned to each cell. Figure 4 illustrates the workflow for deriving rockfall trajectories and their characteristics.
Remote Sens. 2020, 12, x FOR PEER REVIEW 8 of 24 frequency corresponds to a grid mapped onto the DTM, for which each cell gives the ratio of trajectories number that cross it in comparison with the total number or rockfall trajectories. The map of impacts was calculated for which each cell gives the number of impacts (dissipative impacts) that it contains. On the other hand, rockfall bouncing height was calculated for each cell and the highest height was assigned to a particular cell. The derived rockfall mass and velocity were used to calculate the kinetic energy and the highest value was assigned to each cell. Figure 4 illustrates the workflow for deriving rockfall trajectories and their characteristics.

Rockfall Hazard Assessment
A spatial model was developed to assess rockfall hazard based on the obtained characteristics of kinetic energy, bounce height, impact location, and frequency ( Figure 5). Bounce height, kinematic energy, and impact location expose the dynamic interaction between the blocks and topography during their kinematic process. These factors locate the potential destruction that a rockfall can cause to urban areas, tourist places, and transportation corridors. Therefore, they are important factors in rockfall hazard assessment and in the rockfall barrier design. On the other hand, the frequency is one of the most significant factors in rockfall hazard assessment as it reflects the trajectory's spatial distribution and density over the study area.
However, each factor has altered varying influence on rockfall hazard, therefore each factor must be weighted appropriately. The fuzzy-AHP as a qualitative approach can accomplish well for such applications. The fuzzy-AHP was employed in this research to systematically allocate weightage to hazard factors. A numerical relational scale was used to compare between two parameters. Fuzzy logic permits the input of vague, imprecise, and ambiguous data. Consequently, it is widely utilized in decision-making processes. In order to obtain a proper weight of each factor, a questionnaire was

Rockfall Hazard Assessment
A spatial model was developed to assess rockfall hazard based on the obtained characteristics of kinetic energy, bounce height, impact location, and frequency ( Figure 5). Bounce height, kinematic energy, and impact location expose the dynamic interaction between the blocks and topography during their kinematic process. These factors locate the potential destruction that a rockfall can cause to urban areas, tourist places, and transportation corridors. Therefore, they are important factors in rockfall hazard assessment and in the rockfall barrier design. On the other hand, the frequency is one of the most significant factors in rockfall hazard assessment as it reflects the trajectory's spatial distribution and density over the study area.
Remote Sens. 2020, 12, x FOR PEER REVIEW 9 of 24 created containing the focal target (rockfall hazard) and the four selected criteria (rockfall energy, frequency, impact location, and bouncing height). The questionnaire was sent to many experts in the field of geohazards from different countries (as shown in Table 1), to mimic the uncertainty associated with such qualitative information. Consequently, fuzzy-AHP was employed to obtain the contribution weight of each rockfall characteristic's factor.   However, each factor has altered varying influence on rockfall hazard, therefore each factor must be weighted appropriately. The fuzzy-AHP as a qualitative approach can accomplish well for such applications. The fuzzy-AHP was employed in this research to systematically allocate weightage to hazard factors. A numerical relational scale was used to compare between two parameters. Fuzzy logic permits the input of vague, imprecise, and ambiguous data. Consequently, it is widely utilized in decision-making processes. In order to obtain a proper weight of each factor, a questionnaire was created containing the focal target (rockfall hazard) and the four selected criteria (rockfall energy, frequency, impact location, and bouncing height). The questionnaire was sent to many experts in the field of geohazards from different countries (as shown in Table 1), to mimic the uncertainty associated with such qualitative information. Consequently, fuzzy-AHP was employed to obtain the contribution weight of each rockfall characteristic's factor.

Mitigation Process Suggestions
Rockfall mitigation processes are an indispensable tool to alleviate rockfall hazard. Numerous efforts have been conducted to supply developers with a detailed design guideline of rock slope for mitigation processes [34]. Various rockfall protection measures were developed intending to detect or suspending rockfall from its trajectory, and to resist the probable or modelled influence energy resulting from rockfall modelling. Rockfall is a worldwide phenomenon and cause severe destruction to lives and properties. Thus, rockfall protection is an important matter when developing a new urban area or industrial facility and infrastructure, or when defending exposed valuable targets in exposed areas. The rockfall protection includes hazard assessment, identification of protection solutions, the design of structural countermeasure [35], and finally the maintenance program definition. Such a process requires an accurate assessment of rockfall blocks distribution and probability in the potential sources, probable rockfall trajectories, dynamic quantification of rockfall bouncing height, velocity, impact location, and kinetic energy, and intensity of impacts.
In the current research, protection processes were suggested and assessed based on the obtained results of rockfall frequency, kinetic energy, bouncing height, and impact location. In the Lang area, a barrier with (3 m) height was suggested (Figure 6a), taking into account the rockfall characteristics and the possibility of constructing such a mitigation process. The barrier was assessed in terms of the kinetic energy, velocity, bouncing height, and the time of reaching of the falling rocks. This was done to assess the performance and the correctness of the suggested locations. In the Rapat area, the existing fence (yellow line in Figure 6b), with a height of (2 m) was assessed to check the importance of this fence in reducing the rockfall hazard. This was done in terms of the rockfall kinetic energy, velocity, bouncing height, and the time of reaching of the falling rocks. In the Rapat area, it is difficult to suggest a barrier because there are temples at the slope foot. Therefore, alternative mitigation processes were suggested to reduce the rockfall kinetic energy thus rockfall hazard. This included rock sheds and a catch net.
• Qualitative assessment with sensitivity

Mitigation Process Suggestions
Rockfall mitigation processes are an indispensable tool to alleviate rockfall hazard. Numerous efforts have been conducted to supply developers with a detailed design guideline of rock slope for mitigation processes [34]. Various rockfall protection measures were developed intending to detect or suspending rockfall from its trajectory, and to resist the probable or modelled influence energy resulting from rockfall modelling. Rockfall is a worldwide phenomenon and cause severe destruction to lives and properties. Thus, rockfall protection is an important matter when developing a new urban area or industrial facility and infrastructure, or when defending exposed valuable targets in exposed areas. The rockfall protection includes hazard assessment, identification of protection solutions, the design of structural countermeasure [35], and finally the maintenance program definition. Such a process requires an accurate assessment of rockfall blocks distribution and probability in the potential sources, probable rockfall trajectories, dynamic quantification of rockfall bouncing height, velocity, impact location, and kinetic energy, and intensity of impacts.
In the current research, protection processes were suggested and assessed based on the obtained results of rockfall frequency, kinetic energy, bouncing height, and impact location. In the Lang area, a barrier with (3 m) height was suggested (Figure 6a), taking into account the rockfall characteristics and the possibility of constructing such a mitigation process. The barrier was assessed in terms of the kinetic energy, velocity, bouncing height, and the time of reaching of the falling rocks. This was done to assess the performance and the correctness of the suggested locations. In the Rapat area, the existing fence (yellow line in Figure 6b), with a height of (2 m) was assessed to check the importance of this fence in reducing the rockfall hazard. This was done in terms of the rockfall kinetic energy, velocity, bouncing height, and the time of reaching of the falling rocks. In the Rapat area, it is difficult to suggest a barrier because there are temples at the slope foot. Therefore, alternative mitigation processes were suggested to reduce the rockfall kinetic energy thus rockfall hazard. This included rock sheds and a catch net.

Rockfall Source Areas
The multicollinearity of the conditioning factors was assessed utilizing the VIF process to analyze the correlations among these factors. Table 2 illustrates the variance inflation factor (VIF) values determined among the conditioning factors based on rockfall data samples. The highest VIF

Rockfall Source Areas
The multicollinearity of the conditioning factors was assessed utilizing the VIF process to analyze the correlations among these factors. Table 2 illustrates the variance inflation factor (VIF) values determined among the conditioning factors based on rockfall data samples. The highest VIF is 3.268 for the Sediment Transport Index (STI) factor. To identify the best fit hybrid machine learning model, fifteen various models were compared based on various machine learning approach including LR, RF, ANN, SVM, and KNN. Table 3 demonstrates the results of the comparative experiments, illustrating various accuracy metrics of the models. Generally, the stacking RF-KNN outperform the single and hybrid models. The hybrid model achieved the best training accuracy of 89.30%, validation accuracy of 86.50%, and 0.875 of 5-fold CV-AUC. The stacking RF-KNN model prophesied the rockfall probability ranging from 0 (no probable occurrence) to 1 (very high probable of occurrence). Nevertheless, to facilitate the interpretation of the produced maps, the values of probability were reclassified into five classes using a quantile approach ranging from very low to very high (Figure 7a). The reclassified probability map shows that Lang and Rapat areas are highly probable to encounter rockfall incidents.
The analysis of the inventory dataset revealed that rockfall incidents have occurred where the terrain is steep with a slope angle ranging from 57 • to 81 • . Consequently, a slope threshold was set as 57 • and the slope raster was reclassified based on this threshold Figure 7b, shows the results of probable rockfall sources (red line).
approach ranging from very low to very high (Figure 7a). The reclassified probability map shows that Lang and Rapat areas are highly probable to encounter rockfall incidents.
The analysis of the inventory dataset revealed that rockfall incidents have occurred where the terrain is steep with a slope angle ranging from 57° to 81°. Consequently, a slope threshold was set as 57° and the slope raster was reclassified based on this threshold Figure 7b, shows the results of probable rockfall sources (red line).

Rockfall Trajectories and Raster Modeling
Many rocks can be released at random from each seeder point with diverse initial direction in a specific shape and size and particular space. Few thousands of rocks are modeled from all the polyline sources in the two sub-study areas. A one meter distance interval was set for the seeders along the polyline and tens of blocks are thrown at each location. Figure 7 demonstrates the modeling result of rockfall trajectories.
In the Rapat area, 4300 rockfall trajectories were modeled to assess their kinematic behavior and spatial distribution. In order to evaluate the extent of the modeled trajectories without taking into account the existing obstacle, the modeling was conducted without considering the situated obstacle (yellow line in Figure 8b).

Rockfall Trajectories and Raster Modeling
Many rocks can be released at random from each seeder point with diverse initial direction in a specific shape and size and particular space. Few thousands of rocks are modeled from all the polyline sources in the two sub-study areas. A one meter distance interval was set for the seeders along the polyline and tens of blocks are thrown at each location. Figure 7 demonstrates the modeling result of rockfall trajectories.
In the Rapat area, 4300 rockfall trajectories were modeled to assess their kinematic behavior and spatial distribution. In order to evaluate the extent of the modeled trajectories without taking into account the existing obstacle, the modeling was conducted without considering the situated obstacle (yellow line in Figure 8b). In the Lang area, the highest kinetic energy estimated was 678.05 kJ for a rock diameter of 0.75 m. The difference in energy distribution can be recognized in Figure 9a. In the Rapat area, the highest energy value was observed as 575.89 kJ J for a rock diameter of 0.75 m (Figure 9a).
Rockfall frequency in each site is illustrated in Figure 9b. In the Lang area, the frequency raster was classified into five classes including very-low (grey colour), low (green colour), moderate (blue colour), high (yellow colour), and very-high (red colour) (Figure 9b). Each colour reflects the percentage of the rockfall trajectories for a certain raster cell. The frequency distribution is diverse along the analysis area. For the Lang site, the highest bounce height estimation was 92.7 m with rock diameter 0.75 m while the lowest height was 0.3 m. However, 49% of the falling boulders have a bouncing height ranging from (0-2 m), while 20% of the falling rocks have a bouncing height ranging from (2-5 m). Such information is significant for supporting the mitigation processes design. For the Rapat area, the rockfall bouncing height over the study area was also estimated and demonstrated in Figure 9c. The south-eastern part of the study area and the intersection of the main and secondary roads are exposed to the highest bouncing height of 52.1 m. However, 42% of the falling rocks have a bouncing height ranging from 0-2 m, whereas 19% of them have a bouncing height ranging from 2-5 m. On other hand, 39% of the falling rocks in this scenario have a bouncing height greater than 5. Figure 9d illustrates the distributions of the impact locations in the Lang and Rapat areas. The highest impact points number was 114. For Rapat, the highest value was 135 and the highest impact locations were observed near to the rockfall sources, particularly in the south-eastern part of Rapat. In the Lang area, the highest kinetic energy estimated was 678.05 kJ for a rock diameter of 0.75 m. The difference in energy distribution can be recognized in Figure 9a. In the Rapat area, the highest energy value was observed as 575.89 kJ J for a rock diameter of 0.75 m (Figure 9a).
Rockfall frequency in each site is illustrated in Figure 9b. In the Lang area, the frequency raster was classified into five classes including very-low (grey colour), low (green colour), moderate (blue colour), high (yellow colour), and very-high (red colour) (Figure 9b). Each colour reflects the percentage of the rockfall trajectories for a certain raster cell. The frequency distribution is diverse along the analysis area. For the Lang site, the highest bounce height estimation was 92.7 m with rock diameter 0.75 m while the lowest height was 0.3 m. However, 49% of the falling boulders have a bouncing height ranging from (0-2 m), while 20% of the falling rocks have a bouncing height ranging from (2-5 m). Such information is significant for supporting the mitigation processes design. For the Rapat area, the rockfall bouncing height over the study area was also estimated and demonstrated in Figure 9c. The south-eastern part of the study area and the intersection of the main and secondary roads are exposed to the highest bouncing height of 52.1 m. However, 42% of the falling rocks have a bouncing height ranging from 0-2 m, whereas 19% of them have a bouncing height ranging from 2-5 m. On other hand, 39% of the falling rocks in this scenario have a bouncing height greater than 5. Figure 9d illustrates the distributions of the impact locations in the Lang and Rapat areas. The highest impact points number was 114. For Rapat, the highest value was 135 and the highest impact locations were observed near to the rockfall sources, particularly in the south-eastern part of Rapat.

Rockfall Hazard Maps
The rockfall characteristics (kinetic energy, frequency, bounce height, and impact locations) derived through raster modelling were used to produce the rockfall hazard maps. In order to produce

Rockfall Hazard Maps
The rockfall characteristics (kinetic energy, frequency, bounce height, and impact locations) derived through raster modelling were used to produce the rockfall hazard maps. In order to produce rockfall hazard maps, a developed spatial model was employed taking into account the effect of each factor based on the fuzzy-AHP outcomes. After collecting the opinions of the experts, the geometric mean was applied and the final pairwise matrix was prepared. The analysis of the pairwise matrix based on fuzzy-AHP reveals that the value consistency ratio (CR) is 0.03 while the values of consistency index (CI) and randomness index (RI) are 0.019357 and 0.58, respectively. The result of fuzzy-AHP shows that kinetic energy has a high influence on rockfall hazard in comparison with frequency, impact location, and bounce height, respectively. The derived weights of rockfall energy, frequency, impact location, and bounce height were 0.48, 0.3, 0.1, and 0.12, respectively. The generated maps of rockfall hazard are illustrated in Figure 10.
consistency index (CI) and randomness index (RI) are 0.019357 and 0.58, respectively. The result of fuzzy-AHP shows that kinetic energy has a high influence on rockfall hazard in comparison with frequency, impact location, and bounce height, respectively. The derived weights of rockfall energy, frequency, impact location, and bounce height were 0.48, 0.3, 0.1, and 0.12, respectively. The generated maps of rockfall hazard are illustrated in Figure 10.
In the Lang area, the highest rockfall hazard was estimated in the southern portion of the site at the secondary road and the urban area (Figure 10a). Moreover, high rockfall hazard was observed just in front of a temple and behind the existing fence (brown line in Figure 10a). On the other hand, the rockfall hazard on the main road at the northern part was very-low to low. Whereas, the rockfall hazard in the area just close to the fuel station was from very-low to moderate. However, 25% of the site has a moderate hazard, while 20% of the site has high rockfall hazard and 11% of the study area may be exposed to very high rockfall hazard. Figure 10b shows the rockfall hazard assessment in the Rapat site. The distribution of the highest rockfall hazard can be seen from the intersection of the main road with a secondary road to the southeastern part of the site. The rockfall hazard along the main road is ranging from very low to high. However, the percentage of the area that may face moderate rockfall hazard is 17%; whereas 19% of the study area prone to a high degree of rockfall hazard. On the other hand, 9% of the study area may face very high rockfall hazard.

The Assessment of the Suggested Mitigation Approaches
In the Lang area, a barrier with the height of 3 m was proposed based on the rockfall characteristics (mentioned above) and the applicability of applying a particular approach. The proposed barrier could effectively protect the site with an accuracy of 97%. In the Rapat area, the existing fence was able to stop all rockfall trajectories distributed between the intersection area of the In the Lang area, the highest rockfall hazard was estimated in the southern portion of the site at the secondary road and the urban area (Figure 10a). Moreover, high rockfall hazard was observed just in front of a temple and behind the existing fence (brown line in Figure 10a). On the other hand, the rockfall hazard on the main road at the northern part was very-low to low. Whereas, the rockfall hazard in the area just close to the fuel station was from very-low to moderate. However, 25% of the site has a moderate hazard, while 20% of the site has high rockfall hazard and 11% of the study area may be exposed to very high rockfall hazard. Figure 10b shows the rockfall hazard assessment in the Rapat site. The distribution of the highest rockfall hazard can be seen from the intersection of the main road with a secondary road to the south-eastern part of the site. The rockfall hazard along the main road is ranging from very low to high. However, the percentage of the area that may face moderate rockfall hazard is 17%; whereas 19% of the study area prone to a high degree of rockfall hazard. On the other hand, 9% of the study area may face very high rockfall hazard.

The Assessment of the Suggested Mitigation Approaches
In the Lang area, a barrier with the height of 3 m was proposed based on the rockfall characteristics (mentioned above) and the applicability of applying a particular approach. The proposed barrier could effectively protect the site with an accuracy of 97%. In the Rapat area, the existing fence was able to stop all rockfall trajectories distributed between the intersection area of the main road and the secondary one to the southeast of the study area. Thus, it can protect the secondary road. However, the hazard still exists for the temple area which is located between the hill and the fence. Therefore, a rockfall shed was proposed that can protect the area between the main entrance and the temple entrance because this is the most significant part where tourists visit. For the main road, due to it being near the hill, a barrier cannot protect this part because not enough space is there to apply such a process. Therefore, a net location was proposed close to the source areas which can stop the small rocks and reduce the kinetic energy of the big rocks that can reduce the rockfall hazard degree.
The assessment of the proposed barrier at the Lang area with a rock diameter of 0.75 m is illustrated in Figure 11. The frequency of the kinetic energy along the proposed barrier is illustrated in Figure 11a. The highest kinetic energy was observed as 338 KJ with the lowest frequency of 3%. The highest frequency of the kinetic energy is 35% and 32% of rockfall energy ranging from 0-50 KJ and 50-100 KJ, respectively. Whereas, 29% of the reaching rocks have an energy ranging from 100-250 KJ. Regarding the rockfall velocity (Figure 11b), along the proposed barrier, the highest frequency is 35% of rockfall velocity ranging from 10-15 m/s. On the other hand, the lowest frequency was observed as 4% of rockfall velocity ranging from 30-36 m/s. For the bouncing height of the falling rocks along the proposed barrier, 50% of reaching rocks have height ranging from 0-0.5 m, and 25% of reaching rocks have rockfall bouncing height ranging from 0.5-1 m (Figure 11c). On the other hand, 2% of the falling rocks have bouncing height ranging from 2.5-2.7 m. Regarding the time factor, 9% of the falling rocks reach the proposed barrier within 4.5-5.5 s, and 4% of the reaching rocks reach within 10.5-12.5 s. However, 56% of the falling rocks reach the proposed barrier within 5.5-7.5 s (Figure 11d).
In the Rapat area, the existing fence was assessed and alternative mitigation processes were suggested. Figure 12 illustrates the assessment of the existing fence in terms of (a) kinetic energy, (b) velocity, (c) bouncing height, and (d) the reaching time factor. The figure shows that 37% of the reaching rocks along the existing fence have an energy ranging from 0-50 KJ, which is the highest frequency and the lowest energy range in this scenario, respectively. However, 1% of the reaching rocks has an energy ranging from 350-438 KJ, which is the highest kinetic energy along the fence in this scenario. Regarding the velocity of the reaching rocks, Figure 12b demonstrates that 21 of the reaching rocks reach the existing fence with velocity ranging from 0-5 m/s. On the other hand, 1% of the reaching rocks reach the fence with velocity ranging from 35-41 m/s. However, the highest frequency of the velocity is 47% with a velocity ranging from 5-10 m/s. The bouncing height of the reaching rocks along the existing fence was also assessed and presented in Figure 11c. It shows that the highest bouncing height is 7 m with a frequency of less than 7%, whereas 67% of the reaching rocks have bouncing height ranging from 0-1 m. However, the lowest frequency of bouncing height along the existing fence is 1% of bouncing height ranging from 4-5 m. Figure 12d illustrates the assessment of the existing fence in terms of the time factor. It clears that 21% of the reaching rocks reach the existing fence within 6-8 s, whereas, 1% of the rocks reaches the fence within 18-21.5 s. However, the highest frequency of the time factor is 33% of the time ranging from 10-12 s.
Remote Sens. 2020, 12, x FOR PEER REVIEW 17 of 24 Figure 11. The proposed barrier assessment at the Lang area in terms of (a) energy, (b) velocity, (c) height, and (d) time.
The bouncing height of the reaching rocks along the existing fence was also assessed and presented in Figure 11c. It shows that the highest bouncing height is 7 m with a frequency of less than 7%, whereas 67% of the reaching rocks have bouncing height ranging from 0-1 m. However, the lowest frequency of bouncing height along the existing fence is 1% of bouncing height ranging from 4-5 m. Figure 12d illustrates the assessment of the existing fence in terms of the time factor. It clears that 21% of the reaching rocks reach the existing fence within 6-8 sec, whereas, 1% of the rocks reaches the fence within 18-21.5 sec. However, the highest frequency of the time factor is 33% of the time ranging from 10-12 sec.

Rockfall Sources Identification
The multicollinearity of the conditioning factors was assessed utilizing the VIF process to analyze the correlations among these factors. According to the literature, a value of more than four is

Rockfall Sources Identification
The multicollinearity of the conditioning factors was assessed utilizing the VIF process to analyze the correlations among these factors. According to the literature, a value of more than four is considered high collinearity. Therefore, such factors must be removed from the modeling [32]. For the validation of the proposed hybrid models, an assessment was performed by comparing the base models with their stacking models. The stacking RF-KNN outperformed the single and hybrid models. Among the single models, RF accomplished the best result while LR achieved the worst accuracy. After optimizing the hybrid model based on different machine learning classifiers, the stacking RF-KNN model achieved the best performance based on the all accuracy metrics. Therefore, it was employed to produce the rockfall probability map of the whole study area using the conditioning factors and the inventory dataset. The stacking RF-KNN model predicted the rockfall probability ranging from 0 (no probable occurrence) to 1 (very-high probability of occurrence). The reclassified probability map shows that the Lang and Rapat areas are highly prone to rockfalls. In order to detect potential rockfall source areas, the reclassified slope raster was intersected with probability values resulted from the proposed hybrid model (stacking RF-KNN). In addition, the qualitative assessment of produced map shows that the rockfall inventories are typically situated in the regions of those detected as probable rockfall source areas indicating the robustness of the proposed hybrid model.

Rockfall Trajectories and Their Characteristics
The rockfall spatial distribution can be distinguished in both sites. For the Lang area, 6000 boulders were thrown with a spherical shape with a diameter of 0.75 m. The result shows that some rockfall trajectories extended to the urban area crossing the secondary road in the southern part of the study area, another rest just in front of a temple and beyond a situated fence (brown line in Figure 8a), and others crossed the main road in the north region. Nevertheless, some trajectories could reach up to the river that is located next to the slope toe.
In the Rapat area, many rockfall trajectories were seen extending along the main road at the northern region and some of them crossed the secondary road at other parts. Nevertheless, except for the trajectories at the northern part, most of the trajectories tend to rest within the temple's boundaries. Thus, the northern and the middle regions are exposed to rockfall hazard in comparison with other parts. In situ investigation and the inventory dataset have verified the probable areas that are predicted as a high potential of rockfall hazard. Many protection measures exist such as a barrier or fence, rockshed, net, fill, ditch, and berms. Nevertheless, utilizing each process relies on the applicability of a particular method and site condition [4]. In this research, diverse protection measures were suggested according to the obtained rockfall frequency, energy, bouncing height, and the impact locations. In addition, we considered mitigation measures in selecting barrier locations and criteria.
During flying or free-falling movement modes, the falling boulders gain the highest velocity and energy in comparison with other motion modes. However, a rock with bigger diameter has higher energy at a given velocity. This is because the kinetic energy is a function of rock velocity and volume. In the Lang area, the middle and southern regions are exposed to the highest energy. In Rapat area the highest energy was detected in the middle to the south-eastern region and on the main road in the northern region.
The inventory database was used to validate the estimated rockfall frequency. The result shows that the modelling agrees with the recorded rockfall events, as the inventory dataset reflects the realistic spatial distribution of rockfall occurrences. In the Lang area, the frequency distribution is diverse along the analysis area. The part between the hill and the proposed barrier (black line in Figure 9b) is estimated to encounter the highest frequency of rockfall, emphasizing the location correctness of the proposed barrier. In the Rapat site, the main road in the northern part is exposed to the highest frequency while in the south-eastern part, the highest frequency was observed close to the source areas.
Each rockfall trajectory was evaluated in terms of bounce height for each falling rock through calculating the highest bounce height in each raster cell. The highest bounce height was detected on a slope adjacent to the rockfall sources as falling boulders separate from a cliff with a high elevation. However, in the middle of the study area, some trajectories were observed to have high bouncing heights over the suggested obstacle (black line in Figure 9c). This is because of the suggested obstacle being next to the slope foot. Such information is significant for supporting the mitigation processes design. For the Rapat area, the south-eastern part of the study area and the intersection of the main and secondary roads are exposed to the highest bouncing height.
The other significant factor in rockfall hazard assessment is the impact point which reflect the exact location of a falling rock. This is because knowing the impact location and the kinetic energy enables the determination of the destruction and resistance degree and the rockfall vulnerability (degree of loss). Consequently, rockfall risk can be evaluated based on rockfall hazards and vulnerability. Impact location is also highly significant for the planning of protection solutions. For example, an obstacle should be situated adjacent to these locations because a falling block has the lowest kinetic energy with bounce height of zero during the impact phase. Therefore, an obstacle could protect the probable areas effectively at these locations. In both study areas, the highest number of impact points was estimated adjacent to rockfall sources. For the Lang site, the urban area was estimated to be highly affected by bouncing falling boulders (Figure 9d). However, the furthest side of the main road was noted to be completely safe from encountering impact locations. On the other hand, the highest impact locations number were also detected adjacent to the suggested obstacle position (black line), confirming the correctness of the proposed barrier location.

Rockfall Hazard Assessment
Rockfall hazard maps reflect the danger degree rather than single rockfall factors (frequency, bouncing height, kinetic energy, and impact location). For example, a given zone may have a moderate or high value of a hazard but a low value of a particular rockfall characteristics, and vice versa. The rockfall hazard maps were classified into five classes (very low, low, moderate, high, and very high), reflecting the degree of the hazard in a particular area.
In the Lang site, the existing temple in the northern part of the study site is close to the fuel station and the existing fence is prone to a hazard degree ranging from moderate to very high. The open land area in the middle of the site between the fuel station and the urban area is prone to a hazard degree ranging from very low to moderate. In the southern part of the site (close to a temple), the rockfall hazard ranges between moderate to high. For the two temples in the middle and southern parts of the site, the hazard degree ranges from high to very high. This is because of the closeness of these temples from the source areas.
In the Rapat area, the temples prone to rockfall hazards ranges from high to very-high. This is because these temples are located at the foot slope. The secondary road may encounter very low rockfall. This is because most of the falling rocks rest in the zone between the rockfall source areas and the existing fence.

Suggested Mitigation Approaches
Rockfall protection is a significant measure when defending exposed socio-economic related features or during the planning of a new urban development or industrial facilities and infrastructures in the rockfall vulnerable areas. The protection of rockfall involves hazard evaluation, identification of protection method, designing of a structural countermeasure, and finally the definition of a maintenance program [35]. Such processes require an accurate quantification of rockfall block size distributions and susceptibility in the probable source regions, probable rockfall trajectories, intensity and distribution of impacts, statistical variability and magnitude of involved dynamic and kinematic quantities (e.g., velocity, kinetic energy, bouncing height).
In this research, various protection processes were suggested and evaluated based on the obtained rockfall factors including kinetic energy, impact location, bouncing height, and frequency. In addition, the sites' conditions were also considered after applying a particular method. For an example, in the Lang area, a small portion in the middle of the site can face some rockfall trajectories due to its proximity to the barrier to the slope and the bouncing heights of some rockfall trajectories at that location. Although this area is open land, further attention can be paid to it and alternative mitigation process can be tested or a barrier can be constructed with a further distance from the source areas.
In the Rapat area, the existing fence could stop all rockfall trajectories distributed between the intersection area of the main road and the secondary one to the southeast part of the study area. Thus, it can protect the secondary road. However, the hazard still exists for the temple area especially between the hill and the fence. Therefore, rockfall sheds were proposed that can protect the area between the main entrance and the temple entrance because this is the most significant part where tourists visit. For the main road, due to it being near of the hill, a barrier cannot protect this part because not enough space is there to apply such a process. Therefore, a new location was proposed close to the source areas which can stop the small rocks and reduce the kinetic energy of the big rocks that can reduce the rockfall hazard degree.
In the Rapat area there is a fence existing within a distance from the rockfall sources with a height of 2 m. Therefore, we did not propose a barrier in this area. However, the existing fence was assessed and alternative mitigation processes were suggested. Figure 12 illustrates the assessment of the existing fence in terms of (a) kinetic energy, (b) velocity, (c) bouncing height, and (d) the reaching time factor.

Conclusions
Rockfall is one of the most hazardous phenomena in mountainous and hilly regions with high and steep terrain. Such incidents vary both spatially and temporally, thus it is challenging to assess these events accurately. Therefore, this research proposed a method for a comprehensive assessment of rockfall at two different scales (regional and local), in terms of source identification, trajectories, their characteristics, hazard, and mitigation processes. The main datasets used in rockfall hazard assessment are rockfall inventory, LiDAR, and GIS layers. The inventory dataset was prepared from different sources including field measurements. LiDAR techniques both aerial and terrestrial were used to capture the topography of the study areas. The ALS produced a point cloud with point density of 9 points per square meter with root-mean-square-error of 0.15 m vertically and 0.30 m horizontally. The TLS generated multi point clouds of millions of vector points with an accuracy of 3 cm. Filtering (to remove undesired or non-ground features) and interpolation processes were applied to generate the final DTM. A high-resolution and accurate DTM with a spatial resolution of (0.5 m) was produced using the LiDAR dataset. The LiDAR data, inventories, and the GIS database were finally registered and projected into the same projection system (UTM WGS84). Accordingly, various rockfall conditioning factors were derived including morphological and hydrological factors. Fifteen conditioning factors were prepared from different sources and the multicollinearity was assessed among these factors. The result shows that all of these factors are significant for rockfall probability assessment. A hybrid machine-learning model (stacking) were optimized based on diverse machine learning classifiers (LR, RF, KNN, SVM and ANN) using the inventory dataset and the prepared conditioning factors. The hyper-parameters of these classifiers were optimized using the grid search method. The stacking RF-KNN outperforms the other models with an accuracy of (89%), (86%), and (87%) based on training, validation, and cross-validation datasets, respectively. The rockfall probability map that produced through stacking RF-KNN was intersected with the reclassified slope raster based on a specific threshold to produce rockfall sources for the whole study area (regional scale).
The result shows that the Lang and Rapat areas are highly prone to rockfall hazard. Therefore, a developed 3D rockfall kinematic model was calibrated and employed to assess rockfall trajectories and their characteristics in terms of velocity, kinetic energy, frequency, bouncing height, and impact locations. Consequently, a developed spatial model combined with fuzzy-AHP was utilized to estimate rockfall hazard in the analysis areas. In addition, mitigation processes were suggested and assessed and the results prove that these processes can efficiently eliminate or reduce rockfall hazard. The proposed method can assist the decision-makers or the developer of urban areas through providing a detailed assessment of rockfall hazard in a particular area. The generated maps can be used to have a sustainable environment, avert more urbanization in hazardous regions, and to decrease the destruction and fatalities in case of a rockfall incident. Moreover, planners and governments can use the outcomes obtained through this research to distinguish safe areas for residents, support first responders in emergencies, and improve the strategies of urban planning. This information can reduce the necessity to conduct in situ investigation by agencies such as surveying departments. On the other hand, the obtained results can facilitate design in developing barriers and suggest different mitigation processes.