Rockfall feature investigation and kinematic simulation based on nap-of-the-object photogrammetry and GIS spatial modeling

On-site investigation and numerical simulation are the main methods for assessing the risk of rockfall disasters. Unmanned aerial vehicles (UAVs) have been widely used because they can quickly obtain geospatial data. However, in a complex terrain environment, it is difficult to ensure reliable positioning accuracy and resolution of geospatial data collected or constructed at the same flight altitude. Therefore, it is difficult to complete the detailed investigation of rock mass characteristics, and the reliability of the simulation results of rockfall motion characteristics will also be affected. This article summarizes a new image acquisition method and three-dimensional (3D) modeling ideas. Furthermore, in this study a risk assessment is conducted of a slope where a rockfall disaster has occurred. Through this method, the real-world 3D model with a positioning accuracy of less than 5 cm and digital surface model (DSM) data with sub-centimeter resolution of the research area are obtained. Using these data, the rock mass characteristics and the mechanism of disaster formation are analyzed, and the motion characteristics of potential rockfall are numerically simulated in a geographic information system (GIS). After comparing and analyzing the simulated superior trajectory with the on-site terrain, it is found that the trajectory simulated using high-resolution terrain data can better represent the actual movement of falling rocks. The method and process summarized in this paper can provide technical reference for the investigation and evaluation of the risk of rockfall disasters under complex terrain conditions.


Introduction
A rockfall is common geological disaster and a natural phenomenon in which some dangerous rock masses on a slope roll down after being subjected to a force. It is usually characterized by unpredictability and high-speed motion and can seriously threaten villages, communication lines, and traffic. Even if a small boulder rolls off a slope from a height of tens of meters, it can cause casualties (Chau et al. 1998;Choi et al. 2009;Vanneschi et al. 2019; Extended author information available on the last page of the article Youssef et al. 2015). The mechanisms of inducing rockfall disasters are usually diverse and include freeze-thaw erosion, earthquakes, vegetation growth, and artificial excavation. These factors lead to differential weathering of rock masses, joint openings, pore pressure increases, and stress imbalances, which lead to rockfall disasters (Abebe et al. 2010;Chen et al. 1994;Dorren 2003;McCarroll et al. 1998;Tong et al. 2022aTong et al. , 2022bValagussa et al. 2014;Wang et al. 2022a, b, c;Wasowski and Del Gaudio 2000). Once a rockfall occurs, its trajectory and energy are controlled by multiple complex factors, including the size of the rock (Haas et al. 2012), shape of the rock (Buzzi et al. 2012;Vijayakumar et al. 2011Vijayakumar et al. , 2012, and surface properties of the slope (roughness, material properties, and covering) (Basharat et al. 2018;Wang and Lee 2010). The magnitudes and frequencies of different rockfall hazards vary in space and time, and the associated hazards often cannot be eliminated (Lan et al. 2010). Therefore, it is extremely important to study rockfall trajectories and energies through kinematic simulations to assess the risks of a rockfall disaster.
The methods and processes used for analyzing the spatiotemporal motions of rockfalls primarily include field investigations, theoretical analyses, and numerical simulations. First, through field investigations, the characteristics of dangerous rock areas, historical rockfall trajectories, and slope terrain characteristics have been identified, and rock mass stability and failure modes have been analyzed. Then, the trajectories, kinetic energies, and threat ranges of rockfalls were obtained by two-or three-dimensional (2D or 3D) numerical simulation methods (Bozzolo et al. 1988;Lan et al. 2007Lan et al. , 2010Li and Lan 2015;Radtke et al. 2014;San et al. 2020;Sun et al. 2017;Wang et al. 2022a, b, c;Wei et al. 2014); these results provide decision-making bases for subsequent disaster prevention and control (Choi et al. 2009;Fanos and Pradhan 2019;Vanneschi et al. 2019). However, describing rock mass characteristics using traditional field investigation methods becomes challenging when faced with difficult-to-approach slopes. In rockfall hazard assessments, the accuracy of the rockfall source location and the accuracy of the geographic data are two important factors that affect the reliability of the simulation results. The lack of geospatial slope data severely limits the application of numerical simulation methods (Lan et al. 2010). Sturzenegger (2009) proposed a combination of field surveys and ground remote sensing technology to investigate and analyze rock mass characteristics in 2009, and there have been many applications of this method since then. However, for a rock mass on a high and steep slope, since the collection point of ground remote sensing data is at the bottom of the slope, when the location distribution of the rock mass is relatively high, the geometric information, spatial distribution characteristics, and location of the rockfall source cannot be accurately determined. High-resolution geographic data for slopes were also unavailable. The research and application of certain technologies, such as satellite remote sensing and interferometric synthetic aperture radar, have become relatively mature for use with geological disasters (Liu et al. 2022;Sun et al. 2016;Wang et al. 2021). However, these two methods are more suitable for early identification and monitoring of large-scale disasters, and for small-area single geological disasters, these methods have the obvious limitations of redundant data collection, processing difficulty, and high cost (Lian et al. 2020).
With the rapid development of unmanned aerial vehicle (UAV) technology, near-Earth aerial remote sensing technology (with flying altitudes within 1000 m) overcomes the limitations of ground remote sensing. Oblique photogrammetry technology based on lowaltitude remote UAV sensing is convenient and flexible and is not affected by terrain during operation. This method can produce a 3D reality model that can reflect the real texture, color, and geometric shape of the modeling subject (Giordan et al. 2015;Ming et al. 2019;Vanneschi et al. 2019;Youssef et al. 2015). Additionally, high-resolution digital surface models (DSMs) (Mora et al. 2019) and digital orthophoto maps (DOMs) (Manfreda et al. 2018(Manfreda et al. , 2019 can be acquired using this method. However, in Tibet and other remote alpine and valley areas, the GPS positioning accuracy of UAVs is low. The accuracy is also affected by the terrain environment and image collection methods, making it difficult to guarantee the accuracy of image data collected in the field. When capturing images, due to the large vertical height difference of the slope, when collecting images by the conventional method of fixed flight height, the farther the area is from the camera position, the lower the image resolution. When building a 3D reality model, due to the various advantages and characteristics of different modeling methods, it is still difficult to obtain a 3D reality model that can satisfy both position accuracy and clear texture if only one method is used to process images. When investigating the characteristics of rock mass, the low resolution of the model makes it impossible to identify fine joints and fractures. Low positioning accuracy will lead to large errors in the interpretation of structural surface information, and even the statistics of the length and width of cracks may not match the actual ones. Completing the modeling work with high positioning accuracy and high resolution in complex environments requires a professional background in the field of surveying and mapping science, which is difficult for general geological researchers. Therefore, the comprehensive investigation and research of UAV technology in complex environments such as disaster characteristic investigation, formation mechanism analysis, and risk assessment has not leveraged its greatest advantage. Therefore, this paper combines nap-of-the-object photogrammetry (He 2019;Wang et al. 2022a, b, c;Wang et al. 2022a, b, c) and numerical simulation methods to summarize a set of detailed investigation processes suitable for rockfall disasters in complex terrain environments. According to this set of process methods, a hazard investigation and study were completed on a collapsed rockfall in Chaya County, Tibet. Unlike the oblique photogrammetry method, nap-of-the-object photogrammetry achieves 3D capture of the target by planning the route in the vertical direction. During the capture process, the drone flies along the surface, and the capture angle and distance can be adjusted according to the elevation of the slope. The higher-resolution geospatial data obtained ensure the fineness of the identification of rock mass characteristic information and the reliability of the simulation results of rockfall motion characteristics.

Study area
Changdu is located in the eastern part of the Qinghai-Tibet Plateau, in the high mountains and deep valleys east of the Hengduan Mountains. The altitude ranges from 3000 to 6378 m, and the terrain is high in the north and low in the south. The tectonic location belongs to the Changdu-Lanping bidirectional back-arc foreland basin of the Tibet-Sanjiang orogenic Sanjiang Arc Basin, and the Changdu Basin is sandwiched between the Jinsha River junction and the Lancang River junction (Fig. 1). Since 1950, there have been five earthquakes with magnitudes of Ms ≥ 5.0 in Changdu and its surrounding areas. The earthquake with the most complete record and the largest degree of damage is the earthquake of Ms = 6.0 that occurred in 1979. As of 2018, there were a total of 2192 geological disasters and hidden dangers in the Changdu area, with a total of 647 collapses, accounting for 29.52% of the total disasters, which included a total of 403 small rockfall disasters.
Chaya County is located approximately 50 km southeast of the Tibet Autonomous Region. The branch fault of the Lancangjiang fault passes through the town of Chaya (NW-SE direction), the upper wall of the fault exposes the Upper Triassic Duogela Formation (T3d), and the lower wall is the Lower Jurassic Wangbu Formation (J1w) (Fig. 1b). At present, there are a large number of hidden dangers of geological disasters around the county. The north side is dominated by collapses and falling rocks, and the south side presents a large number of large-super-large landslides distributed in stripes (Fig. 1b).
The study area for this investigation was located on a high and steep slope on the northwest side of the county seat. The elevation at the foot of the slope is 3100 m and the elevation at the top of the slope is 3270 m, with an average slope of approximately 45°. The G349 National Road passes through at the foot of the slope. Many rockfall disasters have occurred in this area before, accumulation areas have formed at the foot of the slope, and there are many sources of danger. The complex terrain environment causes significant difficulty for disaster investigation, assessment, and prevention and seriously threatens the safety of passing vehicles and surrounding facilities.

Nap-of-the-object photogrammetry
Nap-of-the-object photogrammetry, a photogrammetry method, uses drones to capture high-definition images by capturing images close to a slope surface. During the measurement process, the capture angle and distance of the drone can be adjusted with the elevation of the slope, ensuring that every image captured by the drone maintains high resolution. Compared with conventional fixed flight height photogrammetry technology, the images captured by the close photogrammetry technology have higher spatial resolution and multi-angle observation advantages. The disaster investigation process based on UAV oblique photogrammetry was primarily composed of four parts: data preparation in the survey area, field photography, internal image data processing, and field verification. Among these, field image capture and internal 3D modeling are critical steps.

Field image acquisition
Nap-of-the-object photogrammetry work follows the "coarse to fine" work idea. First, images were taken according to the measurement method of fixed flight height, and the initial 3D point cloud and digital surface model (DSM) of the study area were obtained. Then, a 3D route based on these initial terrain data was planned, and automatic close flight was completed. In the high-altitude areas of Tibet, the Continuously Operating Reference Stations (CORSs) of the Global Navigation Satellite System (GNSS) are not fully covered. Therefore, the real-time kinematic (RTK) carried out by a UAV cannot obtain reliable position accuracy. To obtain a 3D real scene model that meets both high resolution and high positioning accuracy, image control points must be set up before close-to-flight capture, and the points should be arranged evenly and cover the entire area as much as possible. When collecting the coordinates of the image control points, one GPS device was first set as the base station mode, and then another GPS device was connected as the collecting receiver, which ensured that the accuracy of the collected image control points reached the centimeter level. For high and steep slopes, it is impossible to arrange the image control points in the center of the study area; thus, a highprecision total station was used to encrypt the control points by means of prism-free measurement. The flowchart of image acquisition is shown in Fig. 2.

Construction of 3D model
The sparse reconstruction was based on the structure from motion (SFM) algorithm of computer vision, and the dense reconstruction was based on the multi-view stereo (MVS) algorithm (Ferreira et al. 2017;Meinen and Robinson 2020;Piermattei et al. 2016). The reconstruction process was mainly based on the scale-invariant feature transformation (SIFT) (Lakshmi and Vaithiyanathan 2017;Wang et al. 2013) algorithm, which extracts the feature points of the photographs. Then, the structure was matched using the nearest-neighbor algorithm (NN) (Gbash and Saleh 2017) to calculate the correspondence between feature points of different images. Then, the bundle adjustment (BA) (Huang et al. 2021) method was used to iterate step-by-step to continuously reduce the re-projection error between a projected point and an observed image point and obtain the best camera orientation and 3D point cloud coordinates for the scene.
Currently commonly used modeling software includes ContextCapture, PhotoScan, and Pix4DMapper. It was impossible to obtain a 3D model with reliable position accuracy by only processing images with one program because the advantages of different software under complex terrain conditions are different. For example, the 3D model built by Con-textCapture has high resolution and can truly reflect texture. PhotoScan enables higherprecision aerial triangulation before building a 3D model. To complete the detailed investigation of rock masses, both model resolution and position accuracy must be guaranteed. Therefore, the advantages of various software packages should be considered when modeling, and a reliable 3D model could be obtained by using software that complements each other. The image processing flow is shown in Fig. 3.

Model accuracy assessment
The constructed 3D model should be evaluated for positioning accuracy first, and only when the accuracy is reliable can the rock mass feature information be further interpreted. The accuracy assessment includes the inherent accuracy of the model itself and the actual positioning accuracy of the model in the Earth coordinate system. Model intrinsic accuracy can be checked using image processing software. The actual positioning accuracy evaluation must calculate the mean-square error between the coordinates of the image control points extracted from the model and not involved in the modeling and the coordinates measured by GPS. The calculation formula is shown in Eq. (1): where δ represents the mean-square error, Δs is the difference between the coordinates interpreted in the model and those measured by GPS, and n denotes the number of control points.

Modeling approach for 3D rockfall trajectory in GIS
Disaster simulation analyses and risk assessments are very important for disaster prevention and control, and the commonly used analysis and assessment method is the physical-mechanical model analysis method (Agliardi et al. 2009;Li and Lan 2015;Liu and Lan 2012). Based on Newton's law of motion and collision theory, the numerical simulation mechanics model summarizes the characteristic parameters affecting rockfall through many experimental analyses and uses kinematics and dynamics to simulate the trajectories and velocities of rolling stones (Rammer et al. 2010). In this study, the ArcGIS-based 3D numerical simulation software Rockfall Analyst was used to simulate the motion characteristics of rockfalls in the study area. The model uses the lumped mass method, and the interaction between a rock and a slope is specifically described by two coefficients of restitution (R N and R T ) (Ji et al. 2019; Siming et al. 2009;Wenli andRong 2015Wong et al. 2000) and a coefficient of friction. Rockfall dynamics are not affected by the masses and shapes of the boulders, and only the inherent randomness of the source location and initial orientation is considered in the analysis (Lan et al. 2007).
The software can simulate the three primary motion processes of rockfall flight (collision, bouncing, and rolling), and rolling, and the motion path, height, and speed are simulated by parabolic equations. The displacement and velocity expression are shown in Eqs. (2) and (3): In Eqs. (1) and (2), g represents the acceleration of gravity (9.8 m/s 2 ), (X 0 , Y 0 , Z 0 ) is the initial position of a rolling stone in 3D space, and (V x0 , V y0 , V z0 ) is the initial velocity of the rolling stone in the x, y, and z directions.
During the collision and rebound phase, the first step is to determine the impact point at the end of the rockfall flight, and then, the bounce velocity vector is represented by the particle collision and rebound model (Fig. 4). The bounce velocity vector is calculated according to the restitution coefficients (R N , and R T ). The mathematical expression and bounce velocity vector are defined by Eqs. (4) and (5): (2) The slope surface parameters of the study area should be acquired before the numerical simulation. The accumulation location of rockfall disasters that have occurred is investigated and analyzed, and the slope parameters of the study area are obtained through reverse simulation based on the existing literature and empirical data. The specific process is shown in Fig. 5.

UAV results
The UAV model used in this research was DJI (Da-Jiang Innovations) Phantom4 RTK, the image sensor was a 1-inch CMOS, and the number of camera pixels was 20 million. Because the geometric information for the joint surface must be extracted from the 3D model, not only is the resolution of the 3D model required to be high, but the positional accuracy must also be reliable. First, the coordinates of 12 image control points were collected by combining GPS and high-precision total station, and then, the 3D flight route was planned based on the rough geospatial data obtained at a flight altitude of 150 m. Finally, Numerical simulation process the 3D route was imported into the flight control system of the UAV to complete the close photogrammetry. The flying height of this close photogrammetry was 10 m, and the original image resolution was 3 mm. The photographs were imported into PhotoScan, and eight image control points were selected for aerial triangulation. After the accuracy was qualified, the photographs were imported into ContextCapture to generate a 3D model and obtain a digital surface model (DSM) and a digital orthophoto (DOM) simultaneously. The internal accuracy evaluation table of the model derived from PhotoScan is shown in Fig. 6a. The maximum plane error of the eight image control points involved in the modeling was 3.5 cm, and the maximum elevation error was 6 cm. Four image control points not involved in modeling were used to evaluate the actual positioning accuracy of the model in the Earth coordinate system. Calculated from formula (1), the maximum mean-square error of the plane was 4 cm and that of the elevation was 6 cm. The evaluation results are shown in Fig. 6b. It is obvious that the plane error was within 5 cm whether it was the internal accuracy of the model or the absolute positioning accuracy.

Rock mass characteristics
Dangerous rock mass refers to the isolated rock mass with relatively large elevation difference or the concave steep slope with large free surface in front of the slope body, with fracture development in the slope body, in which the rock mass structure is not complete potential collapse rock mass (Wang et al. 2020a, b;Shi et al. 2020;Wang et al. 2020). The potential collapse rock mass in the study area was mainly siltstone, which was distributed in a "band" shape and was divided into three areas, A, B, and C (Fig. 7). The size of area A was approximately 745 m 3 , the rock mass was cut into blocks, the outcrops are relatively new, and there were obvious signs of overall collapse (Fig. 7a); the size of area B was approximately 1850 m 3 , the surface of the rock mass was weathered into clastic shape, the "X-shaped" shear joints were extremely developed, and many tensile fractures could be seen, with a length of 0.7-2.7 m and a maximum width of approximately 8 cm (Fig. 7b and c); the size of area C was approximately 500 m 3 , the erosion was serious, and many empty surfaces had formed. Signs of caving could be seen from the rock outcrop (Fig. 7d). The height differences of areas A, B, and C and the highway were 22, 48, and 70 m, respectively. As the height of the provenance distribution position increased, the spacing between structural planes gradually decreased, the Fig. 6 Model accuracy evaluation. a Internal relative accuracy of model obtained after eight image control points participate in modeling; x, y, and z represent the three directions of coordinates. b Absolute positioning accuracy of model in eEarth coordinate system based on four image control points not involved in modeling more complex were the secondary structural planes, the more fragmented was the rock mass, and the higher was the degree of weathering. The coordinates of the rock mass surface were extracted from the 3D model, and the directions and inclinations of the three sets of main control structure planes were calculated. The distance between the structural planes was directly measured, as shown in Table 1.
The direction and dip for the three groups of main control structural planes, J1, J2, and J3, as well as the slope aspect and slope (SL) of the slope were drawn on a polar stereographic projection map, as shown in Fig. 8. The intersection point M of the intersection lines of J1 and J2 was located on the outside of the SL, the inclination of the intersection line was smaller than the slope angle, and the inclination of the intersection line was the same as the slope aspect, so the rock mass was unstable; the main slip direction was approximately 138°. The combined intersections of J1, J3 and J2, J3 were opposite to the slope, and the rock mass was relatively  stable. In addition to the main control structural plane affecting the stability of the siltstone, some secondary structural planes could also lead to instability of the rock mass. From a combination of the comprehensive analysis of direction and dip for the structural planes and the on-site investigation, it was determined that the damage form of the siltstone was the falling type. Using the HOEK wedge model, the stability coefficient of the siltstone was calculated to be 0.96 for rainfall and water-filled conditions and 2.16 for dry conditions. These results further show that the rock mass was prone to instability under rainfall conditions. The reference basis for the stability coefficient is shown in Table 2.

Existing rockfalls
There are signs of rockfalls in both provenance areas A and C, and there was rock accumulation at the foot of the slope on the side of the G349 National Highway (Fig. 9). According to the statistics of the 3D model, there were 15 pieces with diameters larger than 1 m, and the maximum diameter was 2.3 m. An analysis of the slope topography and the distribution location for the siltstone found that the rocks falling from area C were more likely to accumulate on the side of the township road. By combining these results with a comprehensive analysis of the rock mass joint spacing and the rock diameters, it was determined that the rock accumulated at the foot of the side slope of the G349 National Highway came from area A. The accumulation at the foot of the slope was near the national highway,  and if rocks collapse from higher places, the scope of the threat requires further study and assessment.

Freeze-thaw erosion and differential weathering
A common rock mass failure mode is soft and hard rocks becoming differentially weathered to form an empty face, which leads to rock mass instability. The rock mass in the study area was silty mudstone and siltstone distributed in the form of "soft-hard-soft" interlayers. Under alternating cold and heat in typical plateau alpine regions, rock masses are easily damaged by freeze-thaw erosion. Since the weathering resistance of siltstone is stronger than that of silty mudstone, the weathering and denudation rate of silty mudstone is accelerated, and thus, an empty face is formed (Fig. 10). The siltstone loses its supporting force, and it is simultaneously cut into unstable blocks by the structural surface. When its own gravity and external forces act to a certain extent, it becomes unstable and falls.

Structural planes
In addition to the three main control structural planes, some secondary structural planes also developed in the rock mass, "X" joints developed in local areas, and tension fractures formed (Fig. 5c). Since the strength of a structural plane is lower than that of the rock block, the combination of different structural planes produces a variety of failure modes. In addition, water easily accumulates in the joints and fissures, and alternating freezing and thawing of ice and snow in the alpine zone expands the gaps between the structural

Changes in the slope shape
The average annual rainfall in Changdu is approximately 488 mm, and it primarily occurs from March to October. After being affected by freeze-thaw erosion, the shape of the slope becomes steeper under the action of rain erosion, which promotes the occurrence of rockfall. In addition, owing to road construction, many excavations at the foot of the slope, as well as changes to the shape of the slope, inevitably affect the stress distributions and will eventually lead to deformation and damage of the rock mass.

Faults
The Changdu Basin in Tibet is located between the Jinsha River junction belt and the Lancang River junction belt. The branch thrust fault of the Lancang River fault zone passes through the county seat (in the NW-SE direction), and the direction and dip angle of the fault plane are 223°-250°∠61°-73°. The "X" shear joint in the provenance area was extremely developed, and it was inferred to be affected by this fault.

Earthquakes
Changdu and the surrounding area are in an earthquake zone, and there have been many earthquakes in its history. Since 1950, there have been five earthquakes with magnitudes of Ms ≥ 5.0, and an earthquake with a magnitude of Ms ≥ 6.0 that occurred in 1979 is the most completely and destructive earthquake recorded since then.

Parameter acquisition
During rockfall simulation analyses performed by Rockfall Analyst, the accuracy of the geospatial data, the values of the coefficients of restitution (R N and R T ), the internal friction angle (Deg), and the location of the rockfall source are the primary factors that affect To ensure more reliable simulation results, the resolution of the 3D model and the DSM obtained in this study reached 5 cm. The location of the rockfall source was directly identified from the visual 3D model by analyzing the rock mass fissures, the air surface formed by weathering and denudation, and the fragmentation situation. The high-precision DSM ensured that the modeled terrain was as close to the actual terrain as possible. R N , R T , and Deg were determined by the inversion simulation method according to the accumulation trajectory of previous rockfalls. First, approximate parameter ranges were determined based on previous research results. Then, the parameter values were adjusted within a certain range to simulate the trajectory of the collapsed rock. Finally, a comparison was made with the rockfall accumulation location obtained from the on-site investigation. If the accumulation position of the simulated track was essentially consistent with the previous rockfall accumulation position, this parameter was used as the simulation parameter for the slope. By adjusting the parameters, the trajectory of the rockfall motion in area A was simulated, and the trajectory end point obtained was in good agreement with the rockfall accumulation position obtained from the field investigation. Therefore, the slope parameters that generated the trajectory were used as the simulation parameters of the slope for the entire study area, as shown in Table 3.

Motion feature analysis
According to the parameters determined by the simulation described above, the potentially caving siltstone blocks at 263 in area B and 188 in area C were simulated, and the motion trajectories, the flight changes, and the speed changes during the rolling processes were obtained, as shown in Figs. 11, 12, and 13. The simulated trajectories show that rockfalls in areas B and C threaten the G349 National Highway and the township road to varying degrees (Fig. 11). On this side of the G349 National Highway, the farthest rockfall track in area B crosses the road, the end points of multiple trajectories are on the road, and the trajectories that do not reach the road also threaten the road. The rockfall track in area C does not reach the foot of the slope, and the overall threat is not large. On the side of the township road, the end points of the rockfall trajectories for areas B and C are nearly all on the road, and some of the trajectories for area C cross the township road to reach the recreational area and the G349 National Highway.
To further study the motion status of the rockfall trajectories and the threat to the road using the simulation, a 3D map of the rockfall trajectories was constructed in this simulation (Fig. 12). The 3D motion trajectories show that the primary rockfall motion forms for the two areas can be divided into rolling, collision bouncing, and flight. On the side of the G349 National Highway, the rockfall that directly threatens the road surface bounces from the foot of the slope and then jumps to the road surface, and some rockfalls directly fly over the G349 National Highway and stop on the outside of the road. On the side of  the township road, the flying altitude during rockfall jumping is relatively high, and many rocks from area C fly directly over the township road, threatening the recreational area. A grid analysis was carried out on the bouncing height and flight speed of the rockfall to study the changes of the bouncing height and flight speed during the rockfall movement (Fig. 13). The maximum flying altitude for the rockfalls in area B is 10 m (Fig. 13a), and the maximum speed is 12 m/s (Fig. 13b). The maximum flying altitude for the rockfalls in area C is 11 m (Fig. 13c), and the maximum speed is 17 m/s (Fig. 13d). The peak flight altitudes of the rockfall trajectories in the two areas occur during the last bounce near the foot of the slope, and the flying speed greatest during the fall after this bounce. The changes in the flight height and velocity vectors corresponding to different nodes of the same rockfall trajectory are correlated.
To further verify the reliability of this numerical simulation, a comprehensive analysis was conducted that combined the simulation results with the actual terrain of the site. On the side of the G349 National Highway, the provenance distribution in area C is high, and the end points of the rockfall tracks are at the middle and upper positions of the slope, while the provenance distribution in area B is low, but the trajectories directly threaten the road. By combining these results with the 3D reality model, it was found that the slope in area C on the side of the national highway is relatively gentle, the weathering debris accumulation is thick, and there are small shrubs growing in the area, which hinder rockfall motion well. The slope is steeper on the side of the township road, and a cliff approximately 3 m high exists at the foot of the slope due to road construction and excavation. These conditions are conducive for faster speeds and higher flight altitudes during rockfall motion, which conclusions are consistent with the simulation results of rocks flying directly over the road and entering the recreational area. Overall, the simulated rockfall trajectory distribution aligns with the geospatial characteristics of the slope, and the results are reasonable.

Advantage trajectory analysis
Two dominant motion trajectories were selected from the trajectories simulated for areas B and C, as shown in Fig. 14a and b, respectively. The maximum bouncing height of the dominant rockfall trajectory in area B is near 6 m. From the perspective of the falling rock traveling distance, compared with the entire trajectory, it is obviously unreasonable to reach the maximum flying altitude at a position where the traveling distance is less than 20 m. By combining these results with the analysis of the slope profile, it was found that the maximum flight altitude node displayed by the trajectory is exactly equal to the cliff height formed by the rock caving in within area A of the provenance, so the flight altitude in this location is relatively high. It is precisely owing to the influence of the cliff that the gravitational potential energy of the falling rock increases during the subsequent rock motion, which further increases the speed. At 45 m along the trajectory, the rock hits the slope for the last time and bounces. The starting point of the bounce is approximately 10 m from the G349 National Highway, and the height is approximately 2 m. At this time, the motion speed reaches its maximum, nearly 9 m/s, and the trajectory finally ends on the G349 National Highway.
The dominant rockfall trajectory in area C is approximately 110 m in total length, the maximum flight altitude node occurs at 80 m along the trajectory, and the maximum height is approximately 9 m. By combining these results with the slope profile, it was found that the highest point occurs approximately 3 m from the rural road, and the falling rock directly leaps over the rural road and enters the recreational area. At this time, the falling rock reaches the ground with its maximum speed, nearly 13 m/s, and then bounces and rolls continuously until it finally stops. The variation trends of the bounce height trajectory and velocity trajectory have a strong correlation with the terrain profile of the motion path, indicating that the simulation results in the high-resolution terrain environment can express the characteristics of rockfall movement more finely and can better represent the actual movement situation.

Discussion
The potential caving siltstone in the study area was distributed at a distance of more than 20 m from the road at the foot of the slope, the overall angle of the slope was approximately 45°, and there were many blind spots from various viewing angles. Traditional survey methods, such as using hand-held compasses and rangefinders, can hardly complete the survey. In the actual investigation process, it was necessary to obtain the coordinate data for the rock mass surface using a 3D model to obtain the geometric information for the rock mass, so the accuracy of the model determined the reliability of the investigation results. During numerical simulation, a higher DSM resolution indicates a better fit of the actual terrain and simulation results that are more aligned with the actual situation. In aerial photogrammetry, the resolution of the photographs and the accuracy of the image control points are the two most important factors affecting the accuracy of the 3D model for terrain conditions that include mountains and valleys. Under such terrain conditions, the height differences in the flight area vary significantly. If the flight altitude is too high, it will inevitably affect the image resolution and thus the model accuracy. In this paper, the technical method nap-of-the-object photogrammetry was applied based on the basis of conventional fixed-height flight survey, which solves the limitation of low resolution and uneven images in the previous slope rock mass survey. The sub-centimeter-level image resolution is sufficient to identify rock fractures, fine joints, rock surface roughness, weathering degree, and other information. To accurately extract the information of joint spacing and even the length and width of small cracks on the basis of identifying rock mass characteristics, a combined modeling method is adopted, and the model accuracy is evaluated by introducing mean-square errors. There are some limitations when using the research method in this paper. The distance between the UAV and the slope was only 10 m when the images were captured. More flights were needed to cover the entire study area with images. For small UAVs, many flights will lead to insufficient endurance. Compared with the conventional method, the number of images collected was greater, so more time must be consumed in real scene 3D modeling, which imposes higher requirements on computer hardware equipment. At the same time, gray texture and voids still appear in the local domain of the model established in this paper, mainly in the water area and places with dense vegetation, which will be further improved in future research.
The shear joints in area B were extremely developed from the identification of the 3D model. Field investigation found that a thrust fault developed at the foot of the slope, with a direction of 220°-240° and a dip of 65°-75°. According to the early geological data, the branch faults of the Lancangjiang fault zone pass through Chaya County (in the NW-SE direction). It is inferred from the survey results that the branch faults of the Lancangjiang fault zone pass through the study area. At present, the stability of the slope is affected by freeze-thaw erosion, differential weathering, structural plane, fault, slope shape change, and earthquake. We consider that the differential weathering and faults of rock mass are the main hazard factors.
Rockfalls in areas B and C of the study area pose a serious threat to the roads and recreational areas on both sides, and preventive measures should be taken. On the side close to the G349 National Highway, the rockfall trajectories indicate bouncing approximately 10 m from the national highway and therefore threats to the national highway. If prevention and control are conducted in this range, the protection height is smallest, but the rockfall speed and the kinetic energy are largest here; therefore, the protective material wears easily. On the side near the township road, building a stone wall at the foot of the slope can prevent some falling rocks from endangering the township road, but it cannot prevent falling rocks from leaping into the recreational area. Based on the terrain profile and the flight altitude trajectory, using a passive net protection at a position 13 m away from the township road would be a more effective prevention and control measure.
The ArcGIS-based 3D simulation software Rockfall Analyst uses the concentrated mass method when simulating rockfall trajectories, so it is assumed that the rock does not disintegrate during the entire rolling process and that its mass remains unchanged. However, in an actual motion process, there must be disintegration when the boulder rolls down from a height of several tens of meters, which makes simulating the motion trajectory more complicated. The assessment of rockfall motion is based on a probabilistic model, and the factors affecting the results of numerical simulations are complex and diverse. To ensure that the location of the rockfall source point was accurate and that the geographic data represented the actual terrain to the greatest possible extent, the probability model and the deterministic model were combined to ensure the reliability of the simulation results in this study.

Conclusions
Through use of nap-of-the-object photogrammetry technology, a visual 3D real-world model and high-resolution DSM and DOM data of the research area were obtained. Based on these data, the rock mass characteristics were interpreted, the formation mechanism of rockfall disasters was analyzed, and the location and volume of potential rockfalls were determined. Using the 3D numerical software Rockfall Analyst, the potential rockfall movement process was numerically simulated, and the threat range of rockfall and the reliability of the simulation results were analyzed. Finally, based on the visual 3D rockfall trajectory obtained from the simulation, as well as the movement speed and flight altitude, a simple prevention and control suggestion was put forward. The research results presented in this paper can provide a certain data basis for the planning and design for disaster prevention and control of collapse and rockfall in Chaya County. At the same time, they provide technical reference for the investigation and research of collapse and rockfall disasters in southeastern Tibet.
Although the study area does not cover the Continuously Operating Reference Stations (CORSs) of the Global Navigation Satellite System (GNSS), high-precision 3D reality models and high-resolution DSMs can still be obtained using UAV technology. It is possible to directly complete the comprehensive research of fine rock mass information investigation, disaster mechanism analysis, and numerical simulation by these results.
Compared with oblique photogrammetry, the resolution of a single image collected based on the nap-of-the-object photogrammetry method is higher. The 3D model constructed in combination with the methods summarized herein guarantees positioning accuracy and resolution at the same time. Although this technical method may lead to redundancy of image data, it still has certain application potential because it can provide sub-centimeter-level 3D model resolution.
DSM data can be directly imported into GIS for simulation, and raster analysis of the results can be performed, due to the fact that Rockfall Analyst is based on the ArcGIS operating environment. Simultaneously, the trajectory of the rockfall can also be displayed in 3D space, which is conducive to visualization research on the movement characteristics of the rockfall.
There were three dangerous areas in the study area, with a total provenance volume of approximately 2690 m 3 . After the rockfall occurred, the trajectory showed a "tree branch"shaped diversion phenomenon affected by the terrain. Rockfalls in areas B and C both threaten the G349 national road and township road, and the township road will be affected more seriously when the rockfall occurs.