Rockfall Hazard Assessment in the Taihang Grand Canyon Scenic Area Integrating Regional-Scale Identiﬁcation of Potential Rockfall Sources

: Frequent rockfall events pose a major threat to the safe operation of the Taihang Grand Canyon Scenic Area (GCSA) in China. The traditional techniques for identifying potential rockfall sources and hazard assessment methods are often challenged in the alpine canyon landform. This study aims to establish an early identiﬁcation framework for regional potential rockfall sources applicable to the canyon region and to assess rockfall hazards in potentially hazardous areas using unmanned aerial vehicle (UAV) photogrammetry. Speciﬁcally, by incorporating high-precision topographic information and geotechnical properties, the slope angle distribution method was used for static identiﬁcation of potential rockfall sources. Moreover, SBAS-InSAR technology was used to describe the activity of potential rockfall sources. Finally, taking the key potentially hazardous area of the Sky City scenic spot as an example, the Rockfall Analyst tool was used to analyze the rockfall frequency, bounce height and energy characteristics based on the high-precision UAV 3D real scene model, and the analytic hierarchy process was introduced to achieve quantitative rockfall hazard assessment. The results show that the potential rockfall source areas in the Taihang GCSA is 33.47 km 2 (21.47%), mainly distributed in strips on the cliffs on both sides of the canyon, of which the active rockfall source area is 2.96 km 2 (8.84%). Taking the scenic spot of Sky City as example, the proposed UAV-based real scene modeling technology was proven to be able to quickly and accurately construct a 3D high-precision model of the canyon area. Moreover, the 3D rockfall simulation showed that the high-energy rockfall area was mainly distributed at the foot of the steep cliff, which mainly threatens the tourist distribution center below. The early identiﬁcation and quantitative evaluation scheme of rockfall events proposed in this study can provide technical reference for the prevention and control of rockfall hazards in similar alpine valley areas.


Introduction
A rockfall refers to the collapse of rock blocks on the surface of steep cliffs and slopes under the action of gravity, with volumes ranging from a few cubic centimeters to thousands of cubic meters [1][2][3][4]. A rockfall is an extremely rapid process and can travel long distances in the form of free-falling, bouncing, rolling and sliding, eventually depositing on talus cones or as individual debris [4,5]. Rockfalls are characterized by suddenness, unpredictability, rapidity, high energy and being catastrophic, being one of the hot issues in geological disaster research [5,6]. Due to its unique natural scenery and cultural landscape, Taihang Grand Canyon Scenic Area (GCSA) is considered one of the "top ten most beautiful grand canyons in China", attracting more than four million tourists every year. Due to the unique geomorphic conditions in the Taihang GCSA, frequent rockfall events pose a major threat to the safe operation of the scenic spot [7][8][9][10].
Identifying the potential rockfall source area is the first step to assessing the rockfall susceptibility and hazard and determining the prevention and control measures for rockfall events [11][12][13][14]. Due to the widespread distribution of potential rockfall points in the canyon area and most of them being located on the high, steep slopes, traditional manual survey methods are often limited in this situation. Technologies such as unmanned aerial vehicle (UAV) photogrammetry and 3D laser scanning have been successively applied to the early identification of rockfall hazards [15][16][17][18]. However, these techniques are generally not applicable to the regional-scale identification of potential rockfall source areas, and their application has certain requirements for the site environment. Therefore, their applications are seriously limited by the rugged and complex terrain in alpine canyon regions. In addition, some scholars have used multi-temporal imagery data acquired by optical satellites and aircrafts to identify potential rockfall source areas at the regional scale [11,[19][20][21]. Considering the high and steep topographic conditions as the primary intrinsic factor for rockfall events [22], researchers have explored various empirical and statistical models for identifying potential rock sources based on a high-precision digital elevation model (DEM). Guzzetti et al. [23] proposed that DEM units with a slope angle greater than 60 • are potential rockfall source areas, based on the statistics of historical rockfall locations. Frattini et al. [24] delineated the rockfall source area by constructing a 3D rockfall trajectory model and considered DEM units with slope angles greater than 37 • as potential rockfall source areas based on the experience of the study area. Previous studies have shown that a threshold angle of 45 •~6 0 • is sufficient to detect potential rockfall source areas [25,26]. It should be emphasized that the threshold angles are generally determined based on field observations, historical data and theoretical assumptions, so the threshold angles are usually different in different study areas [27]. Loye et al. [27] proposed the application of slope angle distribution (SAD) decomposition to identify regional potential rockfall sources, which showed good applicability in different geomorphic regions. Due to the advantages of high precision, working in all weather and a large scale, the interferometric synthetic aperture radar (InSAR) technology has been widely used for surface deformation monitoring, especially for the early identification of landslides [28]. Hosseini et al. [29] used TerraSAR-X Staring Spotlight InSAR data to achieve the detection of active rockfall source areas. Ground-based InSAR has the characteristics of a small scale, high flexibility and high temporal-spatial resolution, showing great potential in the deformation monitoring of key potentially hazardous areas [30][31][32].
Rockfall hazard assessment aims to identify and characterize potential rockfall events and estimate the corresponding occurrence frequency and damage degree [33]. For rockfall hazard assessment, the most commonly used method is to simulate the rockfall trajectory based on physical-mechanical models to obtain the relevant physical and mechanical parameters for rockfall hazard analysis [34,35]. At present, various mechanical models for rockfall trajectory simulation have been widely used in engineering practice, and corresponding mature commercial software solutions have been developed. As the most widely used statistical analysis software for rockfalls, RocFall 5.0 software (Rocscience Inc., Toronto, ON, Canada) calculates the 2D rockfall trajectory based on probability analysis, including motion characteristics such as the bounce height, speed and energy [36]. The RockyFor3D model combines a physics-based deterministic algorithm with a stochastic method to calculate the rockfall trajectory [37]. Lan et al. [38] proposed and built the first 3D spatially distributed rockfall analyst model (RA model) based on the ArcGIS platform, realizing the coupled analysis of the dynamic physical process and the spatial statistics of rockfalls. Using the high-performance parallel computing offered by GPUs, the program STONE established a 3D coupled discrete element model to realize the whole process simulation of large-scale rockfall movement [39]. Acquiring hazard factors such as the rockfall trajectory, kinematic energy and spatial probability and then conducting comprehensive multi-indicator analysis and hazard zoning is the mainstream method for Remote Sens. 2022, 14, 3021 3 of 28 evaluating the hazard of a single rockfall event. For example, Lan et al. [40] applied RA software to evaluate the rockfall hazard by comprehensively considering the characteristics of rockfall source areas and rockfall movement. Wang et al. [41] used RocFall software to carry out rockfall hazard assessment on the rock slope near the Cangjing Cave in Shijing Mountain. To the best of the authors' knowledge, the current rockfall hazard assessment is often limited to a rough regional assessment or a targeted single rockfall evaluation, lacking an integrated solution from regional-scale rockfall identification to single rockfall hazard assessment.
The main purpose of this study is to establish a technical framework for rockfall hazard assessment in the canyon region based on the accurate identification of potential rockfall source areas and activities at the regional scale. First, the geological and environmental background and the distribution and development characteristics of rockfall events in the Taihang GCSA are introduced. Second, the SAD method is used for static identification of the potential rockfall source areas at the regional scale, and the activity of potential rockfall source areas is identified using the InSAR technology, which facilitates the identification of key potentially hazardous areas. Third, taking the scenic spot of Sky City as an example, the RA 1.1 software (Lan H.X., Alberta, Canada) is used to simulate and analyze the characteristics of rockfall movements, and then a quantitative hazard assessment for rockfall events in key potentially hazardous areas in canyon regions is proposed. Finally, we draw our conclusions and discuss future work.

Environmental Geological Background
Taihang GCSA is located in the eastern part of Huguan County, Shanxi Province, China. The total area of the scenic area is 155.9 km 2 , including Baquan Gorge, Qinglong Gorge, Hongdou Gorge, Wuzhi Gorge and Wangmang Gorge ( Figure 1). trajectory, kinematic energy and spatial probability and then conducting comprehensiv multi-indicator analysis and hazard zoning is the mainstream method for evaluating th hazard of a single rockfall event. For example, Lan et al. [40] applied RA software to eva uate the rockfall hazard by comprehensively considering the characteristics of rockfa source areas and rockfall movement. Wang et al. [41] used RocFall software to carry ou rockfall hazard assessment on the rock slope near the Cangjing Cave in Shijing Mountain To the best of the authors' knowledge, the current rockfall hazard assessment is often lim ited to a rough regional assessment or a targeted single rockfall evaluation, lacking a integrated solution from regional-scale rockfall identification to single rockfall hazard as sessment.
The main purpose of this study is to establish a technical framework for rockfall haz ard assessment in the canyon region based on the accurate identification of potential rock fall source areas and activities at the regional scale. First, the geological and environmenta background and the distribution and development characteristics of rockfall events in th Taihang GCSA are introduced. Second, the SAD method is used for static identification o the potential rockfall source areas at the regional scale, and the activity of potential rock fall source areas is identified using the InSAR technology, which facilitates the identifica tion of key potentially hazardous areas. Third, taking the scenic spot of Sky City as a example, the RA 1.1 software (Lan H.X., Alberta, Canada) is used to simulate and analyz the characteristics of rockfall movements, and then a quantitative hazard assessment fo rockfall events in key potentially hazardous areas in canyon regions is proposed. Finally we draw our conclusions and discuss future work.

Environmental Geological Background
Taihang GCSA is located in the eastern part of Huguan County, Shanxi Province China. The total area of the scenic area is 155.9 km 2 , including Baquan Gorge, Qinglon Gorge, Hongdou Gorge, Wuzhi Gorge and Wangmang Gorge ( Figure 1). Under the influence of the rapid uplift of the Taihang Mountains and the intensifica tion of river undercutting erosion during neotectonic movement, an alpine and canyo landform with crisscrossing mountains and valleys was formed in the study area. Th topography of the study area is dominated by steep cliffs and steep slopes, with a Under the influence of the rapid uplift of the Taihang Mountains and the intensification of river undercutting erosion during neotectonic movement, an alpine and canyon landform with crisscrossing mountains and valleys was formed in the study area. The topography of the study area is dominated by steep cliffs and steep slopes, with an average slope of Remote Sens. 2022, 14, 3021 4 of 28 more than 40 • . The terrain of the study area is generally high in the west and low in the east, with elevations of 1300~1600 m and cutting depths of 400~700 m. The study area has a continental monsoon climate with an average annual precipitation of 543.8 mm. The precipitation is mainly concentrated from June to October, accounting for 75% of the annual precipitation.
The study area is located at the southern edge of the Taihang Mountains uplift belt in the south of the North China Plate, and the tectonics in the study area is dominated by the Mesozoic Yanshanian structures [42]. As shown in Figure 2, the faults in the study area mainly extend along the north-northeast direction, and a few faults extend along the northnorthwest direction. The main fault in the eastern part of the study area is the Huangyadi-Yangjiachi-Bada normal fault with a strike of 10 •~1 8 • . The fault group in the north of the study area belongs to the Xianli tectonic magmatic belt with a strike of 310 •~3 40 • . The strata exposed in the study area mainly include the Changcheng system, Cambrian system, Ordovician system, Quaternary system and locally exposed Yanshanian diorite intrusions ( Figure 2). The lithology of the exposed strata is mainly quartzite sandstone, mudstone, limestone and dolomite. average slope of more than 40°. The terrain of the study area is generally high in the west and low in the east, with elevations of 1300~1600 m and cutting depths of 400~700 m. The study area has a continental monsoon climate with an average annual precipitation of 543.8 mm. The precipitation is mainly concentrated from June to October, accounting for 75% of the annual precipitation. The study area is located at the southern edge of the Taihang Mountains uplift belt in the south of the North China Plate, and the tectonics in the study area is dominated by the Mesozoic Yanshanian structures [42]. As shown in Figure 2, the faults in the study area mainly extend along the north-northeast direction, and a few faults extend along the north-northwest direction. The main fault in the eastern part of the study area is the Huangyadi-Yangjiachi-Bada normal fault with a strike of 10°~18°. The fault group in the north of the study area belongs to the Xianli tectonic magmatic belt with a strike of 310°~340°. The strata exposed in the study area mainly include the Changcheng system, Cambrian system, Ordovician system, Quaternary system and locally exposed Yanshanian diorite intrusions ( Figure 2). The lithology of the exposed strata is mainly quartzite sandstone, mudstone, limestone and dolomite.

Spatial Distribution of Rockfall Events
While incorporating a DEM with a resolution of 5 m, the remote sensing interpretation was performed by using 1:10,000 aerial photographs and Google Earth satellite imagery to identify the distribution characteristics of the typical rockfall source areas [19]. Then, these characteristics were used to identify other potential rockfall points in the Taihang GCSA in the ArcGIS platform [21]. In addition, these potential rockfall points have been validated by field investigations.
As shown in Figure 3, a total of 302 potential rockfall points have been identified, including 38 in Wangmang Gorge, 28 in Wuzhi Gorge, 106 in Hongdou Gorge, 96 in Baquan Gorge and 34 in Qinglong Gorge. Potential rockfall points are mainly distributed on both sides of gullies, rivers and roads. The rockfall hazard-prone strata are mainly the quartzite sandstones of the Dahe formation and Zhaojiazhuang formation (7.95%), mudstones of the Mantou formation (10.93%) and limestones of the Zhangxia formation and Gushan formation (60.26%). Moreover, it should be noted that the canyon region is

Spatial Distribution of Rockfall Events
While incorporating a DEM with a resolution of 5 m, the remote sensing interpretation was performed by using 1:10,000 aerial photographs and Google Earth satellite imagery to identify the distribution characteristics of the typical rockfall source areas [19]. Then, these characteristics were used to identify other potential rockfall points in the Taihang GCSA in the ArcGIS platform [21]. In addition, these potential rockfall points have been validated by field investigations.
As shown in Figure 3, a total of 302 potential rockfall points have been identified, including 38 in Wangmang Gorge, 28 in Wuzhi Gorge, 106 in Hongdou Gorge, 96 in Baquan Gorge and 34 in Qinglong Gorge. Potential rockfall points are mainly distributed on both sides of gullies, rivers and roads. The rockfall hazard-prone strata are mainly the quartzite sandstones of the Dahe formation and Zhaojiazhuang formation (7.95%), mudstones of the Mantou formation (10.93%) and limestones of the Zhangxia formation and Gushan formation (60.26%). Moreover, it should be noted that the canyon region is dominated Remote Sens. 2022, 14, 3021 5 of 28 by high and steep cliffs, which severely limits the remote sensing interpretation and field identification of high-locality potential rockfall points. dominated by high and steep cliffs, which severely limits the remote sensing interpretation and field identification of high-locality potential rockfall points.

Types and Characteristics of Rockfall Events
Bedrock slopes are subject to various degrees of weathering, which may lead to fracturing, opening of joints and fissures [22,43]. The strata in the study area are generally characterized by soft-hard alternating structures (Figure 4), and the differences in lithology lead to a difference in the weathering resistance capabilities of the rock masses. Specifically, the underlying weak layers, such as mudstone and shale, are easily eroded to form rock cavities, which provide space for the development of rockfalls ( Figure 5). Simultaneously, it also intensifies the propagation and opening of vertical joint sets and fissures in hard strata. Finally, a rockfall occurs under the action of the predisposing factors.

Types and Characteristics of Rockfall Events
Bedrock slopes are subject to various degrees of weathering, which may lead to fracturing, opening of joints and fissures [22,43]. The strata in the study area are generally characterized by soft-hard alternating structures (Figure 4), and the differences in lithology lead to a difference in the weathering resistance capabilities of the rock masses. Specifically, the underlying weak layers, such as mudstone and shale, are easily eroded to form rock cavities, which provide space for the development of rockfalls ( Figure 5). Simultaneously, it also intensifies the propagation and opening of vertical joint sets and fissures in hard strata. Finally, a rockfall occurs under the action of the predisposing factors. dominated by high and steep cliffs, which severely limits the remote sensing interpretation and field identification of high-locality potential rockfall points.

Types and Characteristics of Rockfall Events
Bedrock slopes are subject to various degrees of weathering, which may lead to fracturing, opening of joints and fissures [22,43]. The strata in the study area are generally characterized by soft-hard alternating structures (Figure 4), and the differences in lithology lead to a difference in the weathering resistance capabilities of the rock masses. Specifically, the underlying weak layers, such as mudstone and shale, are easily eroded to form rock cavities, which provide space for the development of rockfalls ( Figure 5). Simultaneously, it also intensifies the propagation and opening of vertical joint sets and fissures in hard strata. Finally, a rockfall occurs under the action of the predisposing factors.   The exposed strata in Wangmang Gorge, Wuzhi Gorge, Hongdou Gorge in the w and Baquan Gorge in the middle of the study area are characterized by an upper hard a lower soft structure (Figures 4a and 5a). The bottom part is soft mudstone and shale, a the upper part is hard limestone and dolomite. The exposed thickness of the limest and dolomite in the upper part is large, which constitutes the main part of the cany region. Moreover, two sets of vertical joints are usually developed, so tensile and unlo ing failure along the joints can easily occur. Therefore, the rock mass structures are c sidered the predisposing factors of rockfalls. As shown in Figure 5a-d, the failure mo of the rockfalls in this area are dominated by the falling type, toppling type and stagge type [44].
The exposed strata of the Qinglong Gorge in the east of the study area are charac ized by a soft-hard alternating structure (Figures 4b and 5e). The bottom part is h quartzite sandstone of the Dahe formation with two sets of vertical joints, and the fail modes of rockfalls are mainly of the falling type and tension type (Figure 5f). The mid part is the stratum of the Zhaojiazhuang formation with soft-hard alternating sandw structures, and the rockfall failure mode is mainly of the tensile type and falling type (F ure 5g). The upper part is hard quartzitic sandstone of the Changzhougou formation, a the failure mode is mainly of the sliding type and falling type (Figure 5h). Due to presence of underlying weak layers and differential weathering, the upper part of the derlying weak layers in this region is prone to the development of large-scale, unsta tower-shaped rock masses [43].
In addition, the unloading disturbance to the rock slopes caused by human engine ing activities, such as road cutting, excavation and blasting, further intensifies the risk rockfall events in the Taihang GCSA, and the failure mode of these rockfalls is mainly the local small falling type and tension type ( Figure 6). The exposed strata in Wangmang Gorge, Wuzhi Gorge, Hongdou Gorge in the west and Baquan Gorge in the middle of the study area are characterized by an upper hard and lower soft structure (Figures 4a and 5a). The bottom part is soft mudstone and shale, and the upper part is hard limestone and dolomite. The exposed thickness of the limestone and dolomite in the upper part is large, which constitutes the main part of the canyon region. Moreover, two sets of vertical joints are usually developed, so tensile and unloading failure along the joints can easily occur. Therefore, the rock mass structures are considered the predisposing factors of rockfalls. As shown in Figure 5a-d, the failure modes of the rockfalls in this area are dominated by the falling type, toppling type and staggered type [44].
The exposed strata of the Qinglong Gorge in the east of the study area are characterized by a soft-hard alternating structure (Figures 4b and 5e). The bottom part is hard quartzite sandstone of the Dahe formation with two sets of vertical joints, and the failure modes of rockfalls are mainly of the falling type and tension type (Figure 5f). The middle part is the stratum of the Zhaojiazhuang formation with soft-hard alternating sandwich structures, and the rockfall failure mode is mainly of the tensile type and falling type ( Figure 5g). The upper part is hard quartzitic sandstone of the Changzhougou formation, and the failure mode is mainly of the sliding type and falling type ( Figure 5h). Due to the presence of underlying weak layers and differential weathering, the upper part of the underlying weak layers in this region is prone to the development of large-scale, unstable, tower-shaped rock masses [43].
In addition, the unloading disturbance to the rock slopes caused by human engineering activities, such as road cutting, excavation and blasting, further intensifies the risk of rockfall events in the Taihang GCSA, and the failure mode of these rockfalls is mainly of the local small falling type and tension type ( Figure 6).

Identification of Potential Rockfall Sources at the Regional Scale
To assess the hazard of rockfall events and formulate effective prevention an measures, the primary task is to accurately identify the potential rockfall sou Since the study area is a typical alpine canyon landform with high and steep clif 1), traditional optical remote sensing image interpretation and field investigation limited. To solve this issue, this study proposed an identification scheme for id potential rockfall sources at the regional scale suitable for alpine canyon region 7). The scheme includes two parts: static identification and activity detection o source areas. Specifically, based on the high-precision topographic informatio otechnical properties, the SAD method is used for static identification of potentia sources to preliminarily determine the possible location of rockfalls at the regio In addition, it is assumed that the occurrence of rockfall events is related to t mation of the rock surface [32,45]. Then, the small baseline subset InSAR (SBA time series technique is used to detect the active deformation zone at hard-to-re tions in the study area and determine the activity of the potential rockfall sourc

Identification of Potential Rockfall Sources at the Regional Scale
To assess the hazard of rockfall events and formulate effective prevention and control measures, the primary task is to accurately identify the potential rockfall source areas. Since the study area is a typical alpine canyon landform with high and steep cliffs (Figure 1), traditional optical remote sensing image interpretation and field investigation are often limited. To solve this issue, this study proposed an identification scheme for identifying potential rockfall sources at the regional scale suitable for alpine canyon regions ( Figure 7). The scheme includes two parts: static identification and activity detection of rockfall source areas. Specifically, based on the high-precision topographic information and geotechnical properties, the SAD method is used for static identification of potential rockfall sources to preliminarily determine the possible location of rockfalls at the regional scale. In addition, it is assumed that the occurrence of rockfall events is related to the deformation of the rock surface [32,45]. Then, the small baseline subset InSAR (SBAS-InSAR) time series technique is used to detect the active deformation zone at hard-to-reach locations in the study area and determine the activity of the potential rockfall source areas.

Static Identification of Potential Rockfall Sources at the Regional Scale
The core idea of the SAD method is to divide the geomorphologic units using the slope angle thresholds. Strahler (1950) [46] proposed that in an area with uniform lithology, soil, vegetation, climate and morphotectonic history, the SAD of each geomorphologic unit basically conforms to a Gaussian distribution with a low dispersion degree. Moreover, the mean slope angle is determined by a combination of the drainage density, topography and slope profile curvature.
Using the above assumptions, the study area is first divided into several homogeneous morphological areas (HMAs) following identical characteristics. Then, the SAD of each HMA is decomposed into the sum of several Gaussian distributions, in which each Gaussian distribution represents a specific morphological unit (MU). Since the SAD of each MU conforms to a Gaussian distribution, it is also called a GDMU. Finally, the slope angle range of the potential rockfall source area can be determined by analyzing the slope angle threshold of each GDMU [27]. The classification of MUs generally includes the following four types: A technical framework for rockfall hazard assessment incorporating regional-scale identification of potential rockfall source areas.

Static Identification of Potential Rockfall Sources at the Regional Scale
The core idea of the SAD method is to divide the geomorphologic units using the slope angle thresholds. Strahler (1950) [46] proposed that in an area with uniform lithology, soil, vegetation, climate and morphotectonic history, the SAD of each geomorphologic unit basically conforms to a Gaussian distribution with a low dispersion degree. Moreover, the mean slope angle is determined by a combination of the drainage density, topography and slope profile curvature.
Using the above assumptions, the study area is first divided into several homogeneous morphological areas (HMAs) following identical characteristics. Then, the SAD of each HMA is decomposed into the sum of several Gaussian distributions, in which each Gaussian distribution represents a specific morphological unit (MU). Since the SAD of each MU conforms to a Gaussian distribution, it is also called a GDMU. Finally, the slope angle range of the potential rockfall source area can be determined by analyzing the slope angle threshold of each GDMU [27]. The classification of MUs generally includes the following four types: (1) Plain: an area of low slope angles corresponding to the fluvial and fluvio-glacial deposits; (2) Footslope: an area of gentle slope angles at the lower part of a hillslope characterized by alluvial fans, debris flows and landslide deposits; (3) Steep slope: an area of steep slopes in till deposits and rocky outcrops covered with vegetation; (4) Cliff: an area of very steep slopes containing rocky outcrops.
As shown in Figure 8, a valley profile consists of two MUs (footslope and cliff) whose corresponding SAD varies randomly around their mean slope angles (20° for the footslope and 60° for the cliff). This indicates that the SAD of each MU conforms to a Gaussian distribution. Figure 7. A technical framework for rockfall hazard assessment incorporating regional-scale identification of potential rockfall source areas.
As shown in Figure 8, a valley profile consists of two MUs (footslope and cliff) whose corresponding SAD varies randomly around their mean slope angles (20 • for the footslope and 60 • for the cliff). This indicates that the SAD of each MU conforms to a Gaussian distribution. The SAD method mainly utilizes the high-resolution DEM (HRDEM). Since the effective surface area of each raster cell in the DEM depends on its slope angle, the actual terrain area represented by raster cells is determined by weighting each raster cell according to slope angle β. Specifically, the frequency of the slope angle distribution (wβ) is determined by the following equation: where Ahβ is the total horizontal area of the DEM cells with a slope angle β and AHMA is the The SAD method mainly utilizes the high-resolution DEM (HRDEM). Since the effective surface area of each raster cell in the DEM depends on its slope angle, the actual terrain area represented by raster cells is determined by weighting each raster cell according to slope angle β. Specifically, the frequency of the slope angle distribution (w β ) is determined by the following equation: where A hβ is the total horizontal area of the DEM cells with a slope angle β and A HMA is the total area of the HMA considered. The weighted slope angle frequency obtained was used in the subsequent SAD decomposition. The slope angles of each HMA were extracted from the 1:10,000 topographic map of the study area using the ArcGIS 10.1 software (ESRI Inc., Redlands, CA, USA), and the slope angles were reclassified at intervals of 1 • . Then, the frequency of the slope angles was corrected based on Equation (1) and normalized. In this study, the decomposition of GDMUs was performed using the Excel-based application Histofit [27]. Histofit requires manual adjustment of the mean and standard deviation of each GDMU and then iteratively calculates the most likely Gaussian curves, whose sum is close to the SAD of the HMA.
The key of the SAD method for identifying potential rockfall source areas is the selection of slope angle thresholds. For the SAD decomposition obtained through the above process (Figure 9), the slope angle range of each MU is determined by the intersection point between each GDMU on the graph. Specifically, the rules for determining the potential rockfall source areas in this study are as follows ( Figure 9) [47,48]: (1) The intersection between the GDMU cliffs and the GDMU steep slopes is defined as the threshold slope angle (noted as A). All areas with slope angles greater than this threshold are considered potential rockfall source areas; (2) In some cases, the SAD decomposition does not contain the GDMU cliffs. At this time, the highest GDMU is specified as the steep slope MU; (3) In very steep terrain conditions, the SAD decomposition may contain two GDMU cliffs. In this case, the lower GDMU cliff is selected in step 1 to obtain the threshold value A; (4) Rockfalls may occur not only in GDMU cliffs but also in less steep terrain. In addition, if all of the GDMU steep slopes are lumped into the potential rockfall source area, especially when the slope is covered by vegetation or quaternary sediments, the potential source area will be overestimated. Thus, the potential rockfall sources are defined as zones above the mode of the GDMU steep slopes (noted as B) with bare surfaces.
Remote Sens. 2022, 14, x FOR PEER REVIEW 10 of 30 In summary, the potential rockfall source area is identified by all DEM cells above the slope angle (point A), defined by the intersection of the GDMU steep slopes and GDMU cliffs and the rocky exposed areas steeper than point B. Notably, the thresholds of points A and B are directly controlled by the topographic and geomorphological conditions of the study area [27]. The literature review results indicate that the threshold ranges In summary, the potential rockfall source area is identified by all DEM cells above the slope angle (point A), defined by the intersection of the GDMU steep slopes and GDMU cliffs and the rocky exposed areas steeper than point B. Notably, the thresholds of points A and B are directly controlled by the topographic and geomorphological conditions of the study area [27]. The literature review results indicate that the threshold ranges of points A and B are (25 • , 61 • ) and (18 • , 44 • ), respectively [19,[47][48][49]. As the SAD method only uses the topographic information embedded in the HRDEM to identify the potential rockfall source areas, the above method is referred to as static identification of the rockfall source areas, and the results obtained are the static rockfall source areas.
3.1.2. Activity Identification of Potential Rockfall Source Areas at the Regional Scale To obtain more reliable deformation results, Berardino et al. [50] proposed the SBAS-InSAR technique, which can effectively reduce the influence of spatiotemporal decoherence and atmospheric delay and increase the temporal sampling frequency. The basic principle of the SBAS-InSAR technique is to use image pairs with short spatiotemporal baselines to generate interferograms to improve interferometry and then perform the phase-unwrapping of the SAR interferograms [51]. The specific technical flowchart of the SBAS-InSAR technique is shown in Figure 10. Assuming that N + 1 phase SAR data were obtained in chronological order (t0, … tN), one of them was selected as the main image, and the other N phase images were r istered with the main image. By setting reasonable thresholds for the spatial-tempo baselines, M differential interferograms with good coherence were formed, and M sa fied the following inequality (N is odd): For the interferogram j obtained from two SAR images at tA and tB (tA < tB), the int ference phase of the pixel with the azimuth and range coordinates (x, r) can be express as where ϕ(tB, x, r) and ϕ(tA, x, r) represent the phase value of the SAR image at moment and tA, respectively, ∆ϕdisp is the deformation phase along the line-of-sight (LOS) direct between tA and tB, ∆ϕtopo is the residual topographic phase caused by the inaccuracy of Assuming that N + 1 phase SAR data were obtained in chronological order (t 0 , . . . , t N ), one of them was selected as the main image, and the other N phase images were registered with the main image. By setting reasonable thresholds for the spatial-temporal baselines, M differential interferograms with good coherence were formed, and M satisfied the following inequality (N is odd): For the interferogram j obtained from two SAR images at t A and t B (t A < t B ), the interference phase of the pixel with the azimuth and range coordinates (x, r) can be expressed as where φ(t B , x, r) and φ(t A , x, r) represent the phase value of the SAR image at moments t B and t A , respectively, ∆φ disp is the deformation phase along the line-of-sight (LOS) direction between t A and t B , ∆φ topo is the residual topographic phase caused by the inaccuracy of the reference DEM, ∆φ atm is atmospheric phase error, ∆φ orb is the phase error caused by the baseline orbit and ∆φ noise is the random phase noise error. Therefore, the LOS deformation ∆φ disp can be obtained by removing the residual component in the interferometric phase. In addition, positive LOS displacement indicates the ground moved toward the satellite, while a negative value indicates moving away from the satellite.
For M differential interferograms, a set of M equations with N unknowns can be defined as where A is a matrix of dimensions M × N and φ = (φ(t 1 ), φ(t 2 ), . . . , φ(t N )) T is the unknown deformation phase N × 1 vector based on the measurement point. Therefore, the vector of the unwrapped phase values associated with differential interferograms can be expressed as δφ = (δφ(t 1 ), δφ(t 2 ), . . . , δφ(t N )) T . To obtain the deformation velocity, Equation (4) can be rewritten as follows: where B is a matrix of dimensions M × N and v is the mean phase velocity, which can be expressed as If the matrix B is a column rank-full matrix (i.e., M ≥ N), the deformation velocity can be obtained using Equation (5) by applying the least square method. In addition, when the rank of matrix B is deficient (i.e., M ≤ N), the deformation velocity can be estimated using the singular value decomposition (SVD) method. After solving for the deformation velocity, the velocity of each period can be integrated to estimate the shape variable of each period.
The InSAR data used in this study were from the European Space Agency observation satellite Sentinel-1A launched on 3 April 2014. The time span of the data used was from 4 January 2021 to 30 December 2021, with a total of 31 views (see Table 1 for detailed parameters). In this study, the 30-m resolution DEM from the NASA Shuttle Radar Topography Mission was used for removing the topographic phase from interferograms. Considering the undulating terrain and extensive vegetation coverage in the alpine canyon areas, the Delaunay minimum cost flow method [52] was employed for phase unwrapping to reduce the influence of vegetation and other factors. Then, multi-view and Goldstein filters were applied to remove noise, where the multi-view ratio was 1:6 [53]. Considering that the reliability of SBAS-InSAR results can be effectively improved by using stable ground control points (GCPs) [54], this study extracted permanent scatterers [55] as GCPs and introduced them into the processing flow of the SBAS-InSAR method ( Figure 10). The GCPs were used for estimating and removing the residual constant phases and phase ramps that still existed after unwrapping [56].
To estimate the deformation rate and residual topography, a linear model of the deformation coefficient and elevation coefficient for each coherent point was established and solved by the SVD method. Finally, the atmospheric phase at each time point was estimated and removed using time domain high-pass filtering and low-pass filtering techniques to obtain the final deformation results for the time series.

Rockfall Hazard Assessment in Key Potential Hazardous Areas
After locating the potential rockfall source area at the regional scale, a high-precision 3D real scene model can be built for the high and steep cliffs in the canyon region to simulate the rockfall movement process and obtain the disaster-causing factors, such as rockfall trajectory, frequency, bounce height, speed and kinematic energy. Then, the rockfall hazard assessment scheme can be established, which can lay the foundation for the subsequent refined prevention and control of rockfall events.

Construction of a 3D Real Scene Model in the Canyon
Rockfall hazard assessment in key potentially hazardous areas in the canyon region often relies on high-precision georeferenced data, such as the DEM [23]. To solve the problem of constructing 3D models of key potential rockfall hazardous areas in the canyon region, this study integrated UAV oblique photography and nap-of-the-object photogrammetry technology to construct a high-precision 3D real scene model in the canyon region ( Figure 11). The survey used the DJI Phantom 4 RTK, which was equipped with a GPS/GLONASS dual-mode satellite positioning system and inertial measurement unit. The images collected by this UAV contained high-precision coordinate information. Rockfall hazard assessment in key potentially hazardous areas in the canyon region often relies on high-precision georeferenced data, such as the DEM [23]. To solve the problem of constructing 3D models of key potential rockfall hazardous areas in the canyon region, this study integrated UAV oblique photography and nap-of-the-object photogrammetry technology to construct a high-precision 3D real scene model in the canyon region ( Figure 11). The survey used the DJI Phantom 4 RTK, which was equipped with a GPS/GLONASS dual-mode satellite positioning system and inertial measurement unit. The images collected by this UAV contained high-precision coordinate information. Considering the actual situation that the UVA flight is restricted in the narrow canyons, the UAV image acquisition adopted the combination of route flight and manual operation. The level flight route was arranged over the key potentially hazardous areas, and efficient UAV low-altitude oblique photogrammetry was used to realize the data acquisition of large-scale 3D scenes. Specifically, image data were collected from the vertical and four oblique views by the UAV. Considering the undulating terrain and inconspicuous features in the alpine canyon regions, the flight route was set with a heading overlap rate of 85% and a side overlap rate of 80% to ensure the accuracy of the data post-processing. However, the image resolution of the cliffs in the canyon is often low due to the inappropriate shooting distance or angle. Therefore, this study further adopted nap-ofthe-object photogrammetry technology to carry out image acquisition for the steep cliffs in the canyon. Due to the poor GPS signal and limited flight space in the canyon, the UAV was manipulated to fly vertically close to the cliff to obtain image data on the cliff walls.
The UAV images require basically reflecting the detailed contour and texture information of the ground objects, which can provide a data basis for the subsequent construction of 3D real scene models. To ensure the accuracy of the post-processing results, the Considering the actual situation that the UVA flight is restricted in the narrow canyons, the UAV image acquisition adopted the combination of route flight and manual operation. The level flight route was arranged over the key potentially hazardous areas, and efficient UAV low-altitude oblique photogrammetry was used to realize the data acquisition of large-scale 3D scenes. Specifically, image data were collected from the vertical and four oblique views by the UAV. Considering the undulating terrain and inconspicuous features in the alpine canyon regions, the flight route was set with a heading overlap rate of 85% and a side overlap rate of 80% to ensure the accuracy of the data post-processing. However, the image resolution of the cliffs in the canyon is often low due to the inappropriate shooting distance or angle. Therefore, this study further adopted nap-of-the-object photogrammetry technology to carry out image acquisition for the steep cliffs in the canyon. Due to the poor GPS signal and limited flight space in the canyon, the UAV was manipulated to fly vertically close to the cliff to obtain image data on the cliff walls.
The UAV images require basically reflecting the detailed contour and texture information of the ground objects, which can provide a data basis for the subsequent construction of 3D real scene models. To ensure the accuracy of the post-processing results, the original image data was first subjected to distortion correction and color matching. Then, the processed images and their corresponding coordinate information were imported into ContextCapture 10.16 software (Bentley Systems Inc., Exton, PA, USA), and ground control points were added. Finally, the images were encrypted by aerial triangulation to realize the construction of 3D real scene models.

Modeling of Rockfall Trajectories
In this study, RA 1.1 software based on the ArcGIS platform was used to simulate the rockfall movement characteristics [38]. RA 1.1 software can realize direct interaction with ArcGIS functions and can simulate 3D rockfall trajectories. In addition, the geostatistical analysis function in ArcGIS can be coupled with the 3D rockfall trajectories to create a distributed raster model of rockfalls. Compared with other rockfall motion simulation software, RA 1.1 software considers the offset of the motion space when simulating the rockfall trajectory, which is consistent with the actual uneven terrain.
The geometrical and geostructural characteristics of rock masses determine the kinematic conditions of rockfall events, and the slope morphology determines the rockfall motion mode [2,44]. In this study, a 1:2000 topographic map was used to create HRDEM data to extract high-precision information on the elevation, aspect and slope angle. The identification of potential rockfall source areas is the prerequisite for simulating the initiation mode and trajectory of a rockfall. Considering that the potential rockfall source areas in the study area are concentrated on the nearly vertical steep cliff walls, the rockfall source areas with the highest threat at the top of the cliff were selected to simulate the rockfall movement characteristics. In addition, the morphology and properties of the underlying slope surface have an important influence on the barrier effect produced by the falling rocks after collision, which directly determines the characteristics of a rockfall's trajectory, bounce height, kinematic energy and accumulation position. Considering the unique geomorphic characteristics of the canyon area, the supervised classification on the 3D real scene model were used to classify the ground surface material. Moreover, the geotechnical parameters were calibrated for each surface material through back analysis between the case histories and rockfall trajectories.
The RA 1.1 software performed the modeling of rockfall trajectories based on the lumped mass method, which considers the mass of the rock block to be concentrated in one point [57,58]. Rockfall sources (seeders) can be defined as points or polyline or polygon features which are stored in shape files. Notably, the rockfall source area drawn by polylines and polygons is actually represented by uniformly distributed points in space. In this study, the simulated rockfall source area was input into RA in the form of a polyline, and the interval of the point seeders was set to 5 m. In addition, 10 point seeders were positioned separately at each location with a direction interval of 6 • . The weight of a rockfall is determined by the lithology of the stratum where the rockfall source is located. After determining the parameters required for calculation, the RA 1.1 software could be used to simulate the rockfall frequency, bounce height, speed and energy and create the corresponding raster model.

Rockfall Hazard Assessment Process
Using the three raster models of rockfall frequency, bounce height and kinematic energy obtained by the RA simulation, the rockfall hazard assessment mapping could be achieved by the following equation: where the hazard degree D i of each raster cell was obtained by the weighted superposition of the rockfall frequency Rf i , bounce height Rh i and kinematic energy Re i in the cell, and their corresponding weights are a, b and c, respectively. In this study, the analytic hierarchy process (AHP) was used to determine the weights of the dominating factors. The AHP is a multi-objective evaluation decision-making approach that combines qualitative and quantitative methods [59]. The AHP involves a matrix-based pair-wise comparison of the contributions of different factors in rockfalls, which can help decision makers filter a preferred one from a set of alternatives. Specifically, the AHP includes the following five steps: (1) decompose a complex decision problem into several component factors (i.e., frequency, flying or bouncing height and kinematic energy in this study), (2) establish the AHP model according to the affiliation relationship, (3) assign numerical values to determine the relative importance of each factor based on subjective correlation ( Table 2), (4) construct a judgment matrix and (5) compute the eigenvector value of the judgment matrix to obtain the weights of each factor. Table 2. Scale of preference between two parameters in AHP [60]. To avoid logical errors in constructing the judgment matrix, the consistency test needs to be satisfied. The consistency ratio (CR) value was used to determine the random generation probability of the judgment matrix in this study [60]:

Degree of Preference Explanation
where RI is the random consistency index (Table 3) and CI is the consistency index of the judgment matrix, defined as follows: where n is the order of the judgment matrix and λ max is the maximum eigenvalue of the judgment matrix. In this study, it is considered that the judgment matrix has satisfactory consistency when CR < 0.1. Finally, combined with the natural breakpoint method, the final result of the hazard assessment was reclassified into five levels [62]: very low hazard area, low hazard area, medium hazard area, high hazard area and very high hazard area.

Results of Static Identification
Since the local geomorphology of Taihang GCSA is mainly controlled by the geological structure and rock mass properties, it was necessary to consider different combinations of lithology and tectonic units in dividing the HMAs. Based on the 1:50,000 geological map ( Figure 2) and the uniaxial compressive strength of the rocks obtained from field point load tests, the study area was divided into seven main HMAs (Table 4 and Figure 12). Figure 9 shows the results of the SAD decomposition of HMA A.   The specific thresholds of the seven HMAs in the Taihang GCSA were selected as shown in Table 5, and the final static identification results of the potential rockfall source areas are shown in Figure 13. The results show that the potential rockfall source areas are mainly concentrated in the canyon region, accounting for 19.91% of the total area of the study area ( Table 5). The HMAs with the largest potential rockfall source areas are zone E and zone D, whose lithology is dominated by hard limestone and dolomite. Due to the high and steep cliffs and developed unloading joints in HMAs E and D, rockfall events occur frequently in this region. The HMAs with the smallest potential source areas are zone F and zone G. The potential source area is mainly distributed in the artificial outcrops The specific thresholds of the seven HMAs in the Taihang GCSA were selected as shown in Table 5, and the final static identification results of the potential rockfall source areas are shown in Figure 13. The results show that the potential rockfall source areas are mainly concentrated in the canyon region, accounting for 19.91% of the total area of the study area ( Table 5). The HMAs with the largest potential rockfall source areas are zone E and zone D, whose lithology is dominated by hard limestone and dolomite. Due to the high and steep cliffs and developed unloading joints in HMAs E and D, rockfall events occur frequently in this region. The HMAs with the smallest potential source areas are zone F and zone G. The potential source area is mainly distributed in the artificial outcrops in the quaternary alluvium near the valley terraces.   To verify the accuracy of the static identification results, the static source areas of rockfalls identified by the SAD method were compared with the distribution of potential rockfall points determined by field investigation (Figure 14a). In addition, the 3D colored digital terrain model generated using 1:10,000 contours and aerial photographs was used to visually check the accuracy of the identification results. For the key potentially hazardous areas in the canyon area, visual verification was performed on the 3D real scene model constructed by UAV photogrammetry (5-cm ground sampling distance). The validation results show that the potential rockfall source areas identified by the SAD method had good inclusiveness, which was basically consistent with the actual field investigation. To verify the accuracy of the static identification results, the static source areas of rockfalls identified by the SAD method were compared with the distribution of potential rockfall points determined by field investigation (Figure 14a). In addition, the 3D colored digital terrain model generated using 1:10,000 contours and aerial photographs was used to visually check the accuracy of the identification results. For the key potentially hazardous areas in the canyon area, visual verification was performed on the 3D real scene model constructed by UAV photogrammetry (5-cm ground sampling distance). The validation results show that the potential rockfall source areas identified by the SAD method had good inclusiveness, which was basically consistent with the actual field investigation.

Results of Activity Identification
In this study, the ENVI SARscape module was used for interferometric processing of the SAR image data. The thresholds of the spatial baseline and temporal baseline were set to be 45% and 40d, respectively. A total of 87 interference pairs were eventually generated, and the interferograms were connected as shown in Figure 15, with the yellow dot (9 February 2021) image as the super main image. Using the SBAS-InSAR method introduced in Section 3.1.2, the final deformation results for the time series could be obtained ( Figure  16).

Results of Activity Identification
In this study, the ENVI SARscape module was used for interferometric processing of the SAR image data. The thresholds of the spatial baseline and temporal baseline were set to be 45% and 40d, respectively. A total of 87 interference pairs were eventually generated, and the interferograms were connected as shown in Figure 15, with the yellow dot (9 February 2021) image as the super main image. Using the SBAS-InSAR method introduced in Section 3.1.2, the final deformation results for the time series could be obtained ( Figure 16).

Results of Activity Identification
In this study, the ENVI SARscape module was used for interferometric processing of the SAR image data. The thresholds of the spatial baseline and temporal baseline were set to be 45% and 40d, respectively. A total of 87 interference pairs were eventually generated, and the interferograms were connected as shown in Figure 15, with the yellow dot (9 February 2021) image as the super main image. Using the SBAS-InSAR method introduced in Section 3.1.2, the final deformation results for the time series could be obtained ( Figure  16). Referring to the technical specification of InSAR monitoring for geohazards (T/CAGHP 013-2018) in China [63] and the monitoring results of rockfall pre-failure deformation rates [32,[64][65][66], an absolute value of the annual LOS deformation rate more than 10 mm/year was used as a pre-criterion to identify the active source areas of rockfalls. Then, we manually used the aerial and satellite image data to further exclude some raster cells that obviously did not have the conditions for rockfall development or have abnormalities [3,22]. For example, the raster cells distributed in loose deposits, plains, river channels and engineering construction areas were excluded to obtain the interpretation results of active rockfall sources in the study area ( Figure 17). Figure 17 shows that the active source areas of rockfalls are concentrated in Wangmang Gorge in the northwest and Baquan Gorge in the middle, with a total area of 2.96 km 2 , accounting for 1.9% of the total area of the study area. Most of the active rockfall source areas in the Taihang GCSA are located on the mountaintop slope, and a small part are located in the canyon region. As listed in Table 6, the LOS deformation rate of active rockfall source areas was mainly −20~−10 mm/year (59.39%), and the local deformation rate could reach −70~−60 mm/year (0.05%). Referring to the technical specification of InSAR monitoring for geohazards (T/CAGHP 013-2018) in China [63] and the monitoring results of rockfall pre-failure deformation rates [32,[64][65][66], an absolute value of the annual LOS deformation rate more than 10 mm/year was used as a pre-criterion to identify the active source areas of rockfalls. Then, we manually used the aerial and satellite image data to further exclude some raster cells that obviously did not have the conditions for rockfall development or have abnormalities [3,22]. For example, the raster cells distributed in loose deposits, plains, river channels and engineering construction areas were excluded to obtain the interpretation results of active rockfall sources in the study area ( Figure 17). Figure 17 shows that the active source areas of rockfalls are concentrated in Wangmang Gorge in the northwest and Baquan Gorge in the middle, with a total area of 2.96 km 2 , accounting for 1.9% of the total area of the study area. Most of the active rockfall source areas in the Taihang GCSA are located on the mountaintop slope, and a small part are located in the canyon region. As listed in Table 6, the LOS deformation rate of active rockfall source areas was mainly −20~−10 mm/year (59.39%), and the local deformation rate could reach −70~−60 mm/year (0.05%).  To verify the accuracy of the active identification results, the active rockfall source areas identified by SBAS-InSAR technology were compared and analyzed with the potential rockfall points with high hazard areas determined by field investigation (Figure 17).  To verify the accuracy of the active identification results, the active rockfall source areas identified by SBAS-InSAR technology were compared and analyzed with the potential rockfall points with high hazard areas determined by field investigation (Figure 17). In addition, the identified active source areas were visually checked again using the highprecision 3D real scene model obtained by UAV photogrammetry. The validation results show that the active rockfall source areas identified by SBAS-InSAR technology were basically consistent with the field situation.

Identification Results of Potential Rockfall Sources at the Regional Scale
The distribution map of potential rockfall sources in the study area was obtained by summarizing the static source areas determined by the SAD method and the active source areas drawn by the SBAS-InSAR technique. As shown in Figure 18, the total area of the potential rockfall source areas finally determined in the study area was 33.47 km 2 , accounting for 21.47% of the study area. Among them, the static source area was 31.04 km 2 , accounting for 92.74% of the identified source area, and the active source area was 2.96 km 2 , accounting for 8.84% of the identified source area. In this study, we identified the overlapped area of the static source area and active source area as the key potentially hazardous area, with an area of 0.53 km 2 , accounting for 17.91% of the active rockfall source area.

Rockfall Hazard Assessment in the Sky City Scenic Spot
The Sky City scenic spot is a must-see attraction in the Baquan Gorge (Figure 1), with 10,000s tourist per day during the peak season. As shown in Figure 19, the cliff wall where the outdoor sightseeing elevator is located is nearly upright, with a vertical height difference of 208 m. The stratum exposed on the cliff wall is huge, thick oolitic limestone with a weakly weathered surface. Since the bedding orientation of the rock layers in this region is nearly horizontal, and the vertical joints are developed, the rock cavities formed by falling rock blocks of different sizes can be seen everywhere on the steep cliff wall ( Figure  19c). Historical observation data since 2015 indicate that the Sky City scenic spot is a highly prone area for rockfall events.

Rockfall Hazard Assessment in the Sky City Scenic Spot
The Sky City scenic spot is a must-see attraction in the Baquan Gorge (Figure 1), with 10,000s tourist per day during the peak season. As shown in Figure 19, the cliff wall where the outdoor sightseeing elevator is located is nearly upright, with a vertical height difference of 208 m. The stratum exposed on the cliff wall is huge, thick oolitic limestone with a weakly weathered surface. Since the bedding orientation of the rock layers in this region is nearly horizontal, and the vertical joints are developed, the rock cavities formed by falling rock blocks of different sizes can be seen everywhere on the steep cliff wall (Figure 19c). Historical observation data since 2015 indicate that the Sky City scenic spot is a highly prone area for rockfall events. the outdoor sightseeing elevator is located is nearly upright, with a vertical height difference of 208 m. The stratum exposed on the cliff wall is huge, thick oolitic limestone with a weakly weathered surface. Since the bedding orientation of the rock layers in this region is nearly horizontal, and the vertical joints are developed, the rock cavities formed by falling rock blocks of different sizes can be seen everywhere on the steep cliff wall (Figure 19c). Historical observation data since 2015 indicate that the Sky City scenic spot is a highly prone area for rockfall events.

Construction of the 3D Real Scene Model
In this study, a total of 2758 images were collected through the combination of oblique photography and nap-of-the-object photogrammetry technology from the scenic spot of Sky City. The final 3D real scene model constructed by ContextCapture 10.16 software is shown in Figure 11b. Combined with the 3D model, five underlying surface materials were finally identified, and their corresponding average values of the coefficients of restitution and friction are listed in Table 7. The statistical data show that the falling rocks at the Sky City scenic spot are dominated by small rock blocks with a mass of 5-10 kg, with many fresh smash marks visible at the elevator entrance steps and ticket gates (Figure 19d,e). Therefore, the average mass of the falling rocks was set as 10 kg in the simulation. The rockfall trajectory obtained by the RA 1.1 software is shown in Figure 20a. Figure 20a shows that the rockfall trajectories endangered the tourist distribution center, footway, scenic road and parking lot. Due to the obstruction of the southern mountain, the falling rock blocks could move farthest to the foot of the front mountain. The statistical analysis of the 920 simulated rockfall trajectories in this study shows that the numbers for the falling rock blocks stopped and accumulated in the rivers, roads, distribution centers and shrubs were 35 (3.8%), 72 (7.8%), 375 (40.8%) and 438 (47.6%), respectively. This shows that the falling rock blocks deposited in the range of crowd activities reached more than half of the total number of simulated falling rock blocks.

Simulation Results of Rockfall Trajectories
The statistical data show that the falling rocks at the Sky City scenic spot are dominated by small rock blocks with a mass of 5-10 kg, with many fresh smash marks visible at the elevator entrance steps and ticket gates (Figure 19d,e). Therefore, the average mass of the falling rocks was set as 10 kg in the simulation. The rockfall trajectory obtained by the RA 1.1 software is shown in Figure 20a.   Figure 21 shows a typical rockfall trajectory, with a horizontal movement distance of 323.9 m and a falling height of 285.4 m. After the start of the falling rock block, due to its small initial horizontal velocity, the rock block mainly rolls down under its own gravity. In the falling process, the gravitational potential energy is constantly converted into kinematic energy. When suddenly encountering the platform or obstacle, the falling rock block collides and bounces. Finally, the kinematic energy of the falling rock block continues to decrease under the barrier effect of the underlying slope surface, and the movement mode is converted from bouncing to rolling until it stops.  Figure 20a shows that the rockfall trajectories endangered the tourist distribution center, footway, scenic road and parking lot. Due to the obstruction of the southern mountain, the falling rock blocks could move farthest to the foot of the front mountain. The statistical analysis of the 920 simulated rockfall trajectories in this study shows that the numbers for the falling rock blocks stopped and accumulated in the rivers, roads, distribution centers and shrubs were 35 (3.8%), 72 (7.8%), 375 (40.8%) and 438 (47.6%), respectively. This shows that the falling rock blocks deposited in the range of crowd activities reached more than half of the total number of simulated falling rock blocks. Figure 21 shows a typical rockfall trajectory, with a horizontal movement distance of 323.9 m and a falling height of 285.4 m. After the start of the falling rock block, due to its small initial horizontal velocity, the rock block mainly rolls down under its own gravity. In the falling process, the gravitational potential energy is constantly converted into kinematic energy. When suddenly encountering the platform or obstacle, the falling rock block collides and bounces. Finally, the kinematic energy of the falling rock block continues to decrease under the barrier effect of the underlying slope surface, and the movement mode is converted from bouncing to rolling until it stops. In the RA 1.1 software, the rockfall frequency was the number of rockfall trajectories passing through the same raster cell. Since the cliff wall at the scenic spot of Sky City is nearly upright, the falling rock blocks will produce a lateral offset after colliding with the underlying slope surface. This results in the dispersed distribution of rockfall trajectories and the frequency of rockfalls gradually decreasing along the trajectory (Figure 20b). Fig-Figure 21. Comparison of typical rockfall trajectories for the scenic spot of Sky City before and after placing passive protection nets (the location of the typical rockfall trajectory see Figure 20a).
In the RA 1.1 software, the rockfall frequency was the number of rockfall trajectories passing through the same raster cell. Since the cliff wall at the scenic spot of Sky City is nearly upright, the falling rock blocks will produce a lateral offset after colliding with the underlying slope surface. This results in the dispersed distribution of rockfall trajectories and the frequency of rockfalls gradually decreasing along the trajectory (Figure 20b). Figure 20c shows that the high-energy areas are mainly distributed in the front part of the rockfall trajectory. This is mainly due to the steep slope and large height difference of the steep cliff, and thus the potential energy of falling rock blocks is rapidly converted into kinematic energy during free fall. The kinematic energy of a falling block reaches its maximum near the foot of the cliff. The bounce height of a rockfall refers to the relative height of the rock trajectory above the slope surface. Figure 20d shows that the area with the largest bounce height at the scenic spot of Sky City is concentrated near the foot of the cliff. In addition, with the continued occurrence of rockfall collisions, the energy of a rockfall is gradually lost, resulting in the subsequent rockfall motion being dominated by rolling and sliding.

Results of Hazard Assessment
The rockfall hazard assessment for the scenic spot of Sky City used the method proposed in Section 3.2 of this study. The results of the AHP method show that the weights of the rockfall frequency, kinematic energy and bounce height were 0.5400, 0.2971 and 0.1629, respectively. Moreover, the judgment matrix passed the consistency test (CR = 0.0048 < 0.1) ( Table 8). As shown in Figure 22, the RA simulation results show that the high hazard areas and the very high hazard areas are mainly located in the tourist distribution center at the bottom of the cliff, which directly threatens the safety of tourists. The area of the threat area is 3199 m 2 , and the frequency, energy and bounce height of the most dangerous area are 12, 17.6 kJ and 34 m, respectively. Although the parking lots and scenic roads are low hazard areas, attention is still needed to ensure the safe operation of scenic spots. The cliff wall at the scenic spot in Sky City is a natural geological relic with great ornamental value. Large unstable rock blocks have been removed from the region near the cliff elevator. However, to preserve the view of the cliff face, only localized areas are protected with shotcrete ( Figure 19), which does not provide sufficient protection for the tourist activity area. To preserve the ornamental value of the cliff, it is recommended to set up passive protective nets to intercept the falling rock blocks above the gentle slope at the bottom of the cliff. For example, the RA simulation results show that the falling rock blocks can be effectively blocked by placing a 3-m-high passive protection net at the top of the rockfall deposit ( Figure 21). In addition, the monitoring and early warning equipment is suggested to be arranged for continuous observation.  The cliff wall at the scenic spot in Sky City is a natural geological relic with great ornamental value. Large unstable rock blocks have been removed from the region near the cliff elevator. However, to preserve the view of the cliff face, only localized areas are protected with shotcrete ( Figure 19), which does not provide sufficient protection for the tourist activity area. To preserve the ornamental value of the cliff, it is recommended to set up passive protective nets to intercept the falling rock blocks above the gentle slope at

Early Identification of Potential Rockfall Source Areas at the Regional Scale
For the early identification of potential rockfall source areas in the canyon region, this study used the SAD method to extract the angle threshold of each MU to achieve the static identification of potential rockfall source areas and also applied the SBAS-InSAR technique to calculate the small deformation in the study area and identify the activity of potential rockfall source areas, both of which achieved good results. The SAD method is a fast, efficient and systematic global method for identifying potential rockfall sources at a regional scale without considering the influence of the rock mass structure [27]. Moreover, the results of this study also show its excellent application value in evaluating the rockfall susceptibility and hazard. The SBAS-InSAR technique has the advantage of acquiring small deformations and long time series of slow surface displacements [67,68]. Hence, applying it to the activity monitoring of rockfall source areas can effectively identify dangerous rock blocks. In this study, the combination of the above two methods can greatly improve the efficiency and accuracy of the identification of potential rockfall sources, which can provide a reference for the early identification of rockfall events. However, the above framework for identifying potential rockfall sources mainly provides an easy-to-use solution for the rapid identification and delineation of potential rockfall hazard areas at the regional scale [49]. The proposed framework can be used to identify key potentially hazardous areas and clarify the priorities for subsequent more detailed hazard assessment. It should be noted that the SAD method identifies potential rockfall source areas on a morphological basis, thus leading to a rough identification of source area boundaries. In addition, the SBAS-InSAR technique is affected by the high vegetation coverage and steep slope in the alpines canyon region, resulting in reduced accuracy of the identification results, especially for the weak signals in the canyon regions. Moreover, the deformation of the rock surface does not always correspond to the occurrence of rockfall events [3]. Although good results have been obtained by using a threshold of 10 mm/year to determine the active rockfall source areas, the selection of the threshold still needs to be further discussed and verified in future studies. Therefore, the early identification framework of rockfall events proposed in this study cannot completely replace field investigation.

The Influence of the DEM Precision
The previous studies indicate that the identification of static source areas is greatly affected by the precision of the DEM [27]. Since the DEM's surface is a simplification of the real topography, some morphological features may be missing. Specifically, the lower the DEM resolution, the lower the slope angle of steep slopes, which leads to small slopes in the study area often being ignored. Therefore, a larger-scale topographic map can be used in subsequent studies to further improve the accuracy of the identification of potential rockfall sources at the regional scale. In addition, the precision of the DEM also affects the simulation results of the RA 1.1 software [40]. The trajectory and velocity of rockfalls simulated by high-resolution DEMs are more complex and variable, and the falling rock blocks will have a larger lateral offset after collision. In contrast, due to some microgeomorphic features being absent in the low-resolution DEM, the obtained kinematic energy and bounce height of rockfalls are too large, which leads to the predicted rockfall threat range not being able to be applied in practice. The results of this study show that the high-precision topographic data constructed by UAV photogrammetry technology can basically meet the needs of rockfall movement simulation. However, the influence of the DEM's precision on the simulation of rockfall movements should be further studied in combination with the monitoring of data for rockfall events.

Response Relationship between Rainfall and Rockfall Events
Many scholars have shown the influence of rainfall on the rockfall frequency [44]. Field surveys and historical statistics found a significant increase in rockfall events in the study area in 2021, especially during the extreme rainfall event in July [69]. As shown in Figure 23, cumulative displacement curves of the three active rockfall source areas in the study area were constructed and compared with the histogram of monthly cumulative rainfall in 2021. Figure 23 shows that the rainfall in the study area increased sharply in July 2021. In addition, the signs of accelerated deformation growth in the three active rockfall source areas can be clearly observed. Therefore, the good correspondence between the deformation trend and rainfall further illustrates the accuracy of the activity identification results. Among them, the deformation at point C in Baquan Gorge increased most significantly. Notably, the occurrence of rockfall events is controlled by a combination of internal and external factors [3]. The field survey found that the topographic and geostructural characteristics of point C made it more prone to rockfalls under rainfall conditions. Therefore, the deformation of point C continuously increases during the rainy season from July to September. DEM's precision on the simulation of rockfall movements should be further studied in combination with the monitoring of data for rockfall events.

Response Relationship between Rainfall and Rockfall Events
Many scholars have shown the influence of rainfall on the rockfall frequency [44]. Field surveys and historical statistics found a significant increase in rockfall events in the study area in 2021, especially during the extreme rainfall event in July [69]. As shown in Figure 23, cumulative displacement curves of the three active rockfall source areas in the study area were constructed and compared with the histogram of monthly cumulative rainfall in 2021. Figure 23 shows that the rainfall in the study area increased sharply in July 2021. In addition, the signs of accelerated deformation growth in the three active rockfall source areas can be clearly observed. Therefore, the good correspondence between the deformation trend and rainfall further illustrates the accuracy of the activity identification results. Among them, the deformation at point C in Baquan Gorge increased most significantly. Notably, the occurrence of rockfall events is controlled by a combination of internal and external factors [3]. The field survey found that the topographic and geostructural characteristics of point C made it more prone to rockfalls under rainfall conditions. Therefore, the deformation of point C continuously increases during the rainy season from July to September.  Figure 17).
The influence of the rainfall intensity on rockfalls can be explained by different processes, such as water pressure, strength weakening, accelerated weathering and the erosion of supporting material [3]. Although many scholars have shown the influence of rainfall on the rockfall frequency [22], rainfall is not the only causal factor. The geometrical properties and geostructural characteristics of the rock mass conditions also directly con-  Figure 17).
The influence of the rainfall intensity on rockfalls can be explained by different processes, such as water pressure, strength weakening, accelerated weathering and the erosion of supporting material [3]. Although many scholars have shown the influence of rainfall on the rockfall frequency [22], rainfall is not the only causal factor. The geometrical properties and geostructural characteristics of the rock mass conditions also directly control the frequency of rockfall events. Therefore, the correlation between rainfall and rockfall frequency revealed by the SBAS-InSAR technique in this study is only a preliminary result, and further analysis is needed in future studies.

Conclusions
Taking the Taihang GCSA as an example, this study applied the SAD method to achieve the static identification of potential rockfall source areas at a regional scale and then identified the activity of potential rockfall source areas with the help of SBAS-InSAR technology. On this basis, a typical potential rockfall point, named Sky City, was selected to analyze the rockfall movement characteristics by using RA 1.1 software. Finally, a rockfall hazard assessment in the canyon region was demonstrated. The main conclusions obtained are as follows: (1) For the early identification of potential rockfall source areas in alpine canyon regions, a region-scale identification framework combining static identification and dynamic detection was proposed in this study. The proposed framework was proven to be efficient and accurate. (2) The total area of potential rockfall sources in the Taihang GCSA is 33.47 km 2 , accounting for 21.47% of the study area, which is mainly distributed in strips on the cliffs on both sides of the canyon. Among them, the active source area determined by SBAS-InSAR technology is 2.96 km 2 , and their deformation rates are mainly −20~−10 mm/year. The field verification shows that the identification results were basically consistent with the field situation. (3) To meet the needs of high-locality rockfall hazard assessment in the key potential hazardous area of the canyon region, this study proposed a scheme to construct a high-precision 3D real scene model of the canyon region by integrating UAV oblique photogrammetry and nap-of-the-object photogrammetry technology. (4) Raster models of the rockfall frequency, bounce height and kinematic energy were obtained by simulating the rockfall movement process, and the weight of each disastercausing factor was determined by introducing the AHP method. The results of the hazard assessment indicated that the high-energy area of rockfalls at the scenic spots of Sky City was distributed at the foot of the steep cliff and primarily threatened the tourist distribution center below. (5) This study proposed an easy-to-use integrated solution from region-scale early identification to individual evaluation of rockfall events, which can provide a technical reference for rockfall prevention and control in similar alpine canyon regions. In future research, the applicability and accuracy of the InSAR technique in identifying the activity of potential rockfall source areas need to be further improved.