A New Rockfall Simulation Tool based on Unity 3D

Objectives: Rockfall occur frequently in the mountainous areas of Southwest China, causing huge casualties and property losses every year. A rapid simulation method of the collapse process is urgently needed. The commonly used collapse simulation software still has obvious shortcomings, such as low terrain accuracy, no consideration of the structural characteristics of rock mass, as well as block collision and fragmentation. The logic operation based on CPU still limits the calculation speed. In this study, a new method for rapid simulation of the three-dimensional motion process of collapse is proposed. Methods: The Unmanned aerial vehicle (UAV) aerial photography modelling combined with field investigation is used to obtain the slope surface model and determine the characteristic parameters of rockfall. The simulation software of a large-scale collapse movement process is developed by Unity 3D platform which integrates PhysX physics engine and CPU-GPU parallel computing capability. Results: The software can simulate the whole process of collapse, impact, fragmentation, and deposition. It can output the three-dimensional trajectory, velocity, energy and jumping height, which can be used for the engineering design. Finally, Xiejiayan rockfall in Nayong, Guizhou Province has been selected as an application case, to verify the feasibility and suitability in practice.


Introduction
Rockfall is a common type of geological disaster in mountainous areas, which has obvious concealment, sudden and catastrophic characteristics, and is easy to cause casualties and property losses [1,2] . With global warming, there were more and more collapse disasters in southwest China. For example, a largescale collapse of the mountain in Longchang Town, Kaili, Guizhou Province, on February 18, 2013, with a volume of 3×10 4 m 3 , caused 5 people dead and relocated 21 households with 79 persons [3] . On August 28, 2017, another mountain collapse occurred in Pusa Village, Zhangjiawan Town, Nayong County, Guizhou Province, with a volume of about 6×10 5 m 3 , which killed 35 people [4] . On June 24, 2017, a high-level landslide occurred in Xinmo Village, Diexi Town, Maoxian County, Sichuan Province, the Xinmo landslide with a volume of about 4.5×10 6 m 3 , caused 64 farmhouses to be buried, 83 people dead or missing, and blocked more than 1,000 meters of river channel [5] . It can be seen that rockfalls have a great role in mountainous constructions and lives.
The research on the mechanical mechanism and movement characteristics of the collapse movement process include field tests, model tests and numerical simulation methods [6][7][8][9][10] . Numerical simulation methods have the advantages of low cost, simplicity and portability. Then, it is commonly used for rockfall evaluation. Compared to the main software (Table 1), it can be seen these methods still have obvious shortcomings and did not consider the fragmentation process of collapsed bodies. E.g., Rockfall software is often used in the calculation of distance, bounce height and kinetic energy. But it needs to specify the terrain profile manually, which cannot truly reflect the rockfall motion trajectory [11] . Rockfall Analyst is a secondary development component of ArcGIS, which can realize the automated and semiautomated analysis of collapse motion trajectory simulation and risk mapping, but ignores rock fragmentation as well [12] . CRSP-3D is a collapse simulation program designated by the U.S. Federal Highway Administration related to engineering projects, but it does not consider the structural characteristics of the rock mass, and cannot set the location of the dangerous rock [13] . It has certain limitations in simulating large-scale collapse movement [14] . In addition, key parameters such as movement distance, jumping height, and impact force are required in the quantitative risk assessment of collapse. Therefore, further study on the rockfall movement simulation is to be important in scientific and social value. The Unity3D engine is a set of integrated development environments including graphics, sound, physics and other functions developed by the Danish company Unity Technologies. It can create interactive 3D scenes. The operating system is highly compatible. It is widely used in system simulation and virtual reality [15] . The engine integrates the PhysX physics engine, provides simultaneous detection of multiple object continuous collisions and supports GPU parallel computing functions [16] . Based on the theory of rigid body dynamics, it can simulate the behaviour of rigid body bouncing, rolling and sliding, solve the motion speed and energy, and obtain the motion trajectory. We think it might be a new method for solving the three-dimensional simulation of the large-scale collapse movement process. Chen Shushu (2018) has already used the Unity3D platform to develop the movement scene of the debris flow [17] . Chiu (2015) used Unity to create a landslide motion scene and uploaded it to the Web and mobile platforms for an application in science education on geological disasters [18] . However, the above applications have not yet introduced the physics engine technology. As a result, the simulation results lack physical support and the effect is poor.
Given this, according to a large number of collapse prototype cases investigation, this paper uses a physics engine to ensure the relationship between collapse dynamics theory and computer simulation. A rockfall simulation tool has been developed based on Unity3D to realize a physical process of fragmentation in the movement process. The Xiejiayan rockfall in Guizhou Province was selected as a case application to verify its reliability and applicability, which can provide a new simulation method and evaluation method for the prevention and mitigation of rockfall disasters.

Theory and Simulation Method of Rockfall Dynamics
The simulation of the rockfall movement process is mainly to deal with the interaction between rockfall and slope, which may involve multiple collisions, rebounds, rolling, sliding and impact fragmentation processes [19] . When falling rocks hit the slope surface, if the average slope angle of the collision point is less than 45 degrees, it is considered that the rock starts to roll, and if it is greater than 45 degrees, we are more likely to see a rebound. In the numerical simulation of the movement process, the rock block and the slope are regarded as rigid bodies. The energy loss in the collision between rigid bodies is usually processed by a contact function including friction and recovery coefficients [20] . The friction and recovery coefficients use the parameters of the physical material assigned to the object for effective control in the Unity3D engine [21] . Physical materials allow users to control the interaction between objects through the dynamic friction coefficient, static friction coefficient, and bouncing coefficient used by the physics engine PhysX. Depending on the state of motion before the collision and the geometry of the slope and falling rock objects, each collision may contain multiple contact points ( Figure 1). The linear velocity and angular velocity of the block that bounces off the slope surface after the collision depends on the product of the iterative solution of the velocities of all contact points in each collision under the constraints of contact and friction.  Fig.1 Conceptual diagram of the multiple collisions, and contact points generated during a single rock-slope interaction. v and ω variables represent incident (i) and rebound (r) translational and rotational velocities.
The linear velocity and angular velocity of each contact depend on the friction and recovery equations. During the collision, the coefficient of friction (μ) is used to determine the ratio of the normal force (N) between the simulated rockfall surface and the slope surface. The ratio is used to calculate the friction force when the block moves. The acting director of the friction force (Ff) is opposite to the tangential component of the input velocity vector. Solve the friction force in each contact and use it to calculate the linear velocity and angular velocity after the contact. Under the standard definition of friction angle, the friction coefficient (μ) is usually equal to the tangent of the friction angle (φ). In the Unity environment, the surface collision is broken down into two contact points, causing the applied friction to double ones [22] . Therefore, the effective friction angle is equal to the arctangent angle twice the unit friction coefficient (Eq. 1). The formula for the friction force applied by each contact point (Eq. 2) and the conversion between the friction coefficient and the friction angle are shown below. Unlike the coefficient of dynamic friction applied to the tangential component of velocity, the coefficient of restitution RV only works in the normal direction. This coefficient is equal to the ratio of the relative speed of the object before and after each collision (Eq. 3).
After analyzing the friction and recovery coefficient of the engine, to simulate the rockfall fragmentation in the movement process, a component called a "fixed joint" was used in the engine ( Figure 2) to connect the rock object. The two rock objects connected by a fixed joint are restricted by constraints and cannot move relative to each other in any direction, nor can they rotate. When the external force does not reach the threshold of the fracture force of the fixed joint, it will keep moving relatively. If the external force of the rock block exceeds the threshold, the fixed joint will be disconnected, and the two rock objects will no longer be constrained to separate from each other. It can simulate the fracture of the structural plane caused by the external force of the rock block exceeding the strength of the structural plane.

Three-dimensional simulation of the rockfall movement process
The three-dimensional rockfall movement process simulation software is composed of three modules: slope surface model build, dangerous rock mass fine model and rockfall movement process simulation. The technical route process is shown in Figure 3.

Slope surface model build
There are many ways to build a slope surface model, including reading the grey map to generate the terrain model and using the Sketch Up software to make the 3D model based on a terrain map. This paper proposes a simple and accurate method for building a slope surface model. First, Unmanned Aerial Vehicle is used to obtain digital images, perform uniform colour distortion correction, and import Context Capture for aerial triangulation. The new reconstruction project builds a triangulated mesh model in FBX format and a JPG format texture map, which is imported into Unity3D. Then, the coordinate should be adjusted and collision components should be added to complete the slope surface model.

Rockfall failure mass model build
The rockfall failure mass is the main body of rockfall movement simulation, and various parameters need to be set to the model beforehand. We developed a structural plane measurement program on the slope surface model, which was developed in C++ and OSG engine. Using the OSGB format model loading plug-in provided by OSG, the model can be loaded directly into the scene. In terms of user interaction, QTGUI is used as the bottom-level dependency library for GUI development, and the software interface is designed as shown in Figure 4. After the user loads the model and selects three points that are not collinear on a flat surface, the software can read the three-point space coordinates and use the three-point method [23] to calculate the dip direction and dip angle.  Fig.4 Interface of a structural plane measurement program

Rockfall movement simulation
Rockfall movement simulation is a core module of the software. It simulates the whole movement process of the rockfall failure mass from falling from instability to stopping at the bottom of the slope and obtaining its movement characteristic parameters. The specific steps are: use C# to write a script program to obtain the three-dimensional coordinates, speed, distance from the underlying surface and other motion characteristic parameters of the object, and add the script program to each dangerous rock block, output the data classification to a table file, and finally analyze and process the data to obtain information such as kinetic energy, the movement path of the block stone, and the ground elevation under the movement path of the block stone.

Geological setting
The Xiejiayan rockfall is located at a high mountainous landform with strong weathering and unloading. The elevation of the slope foot is about 1740m, the top is about 1884m. The elevation of the source area of rockfall failure mass is about 1840 meters, and the altitude difference is 140-150m. The rocks are limestone, milestone and mudstone, with a dip inclination is 189-193°, and the slope inclination is about 228°. The slope is a dip-bedding slope with a dip angle of 8-18°. The slope is steep with an angle of 70-85°, forming a cliff with few vegetations and lots of rockfall failure masses. The rockfall deposit is 15~25°, which is covered by low shrubs. The bottom of the slope is scattered with collapsed rocks before, and a passive protection net has been built to protect the residential houses.  Here, rockfalls often happen and to be scattered on the bottom of the slope surface. Based on field investigation and UAV detection, the deposition is mainly distributed close to shrubs, particular some large-size rocks. The overall view of the source area and deposition area of rockfall is shown in Figure  6. It can be seen, most of the rock volume on the slope is less than 1m³, some large-size rocks can be up to 5.0m 3 , the largest volume can be up to 24m 3 . To understand the mechanism of rockfall, the discontinuities in rock mass have been checked by a self-developed structural joints measurement program compared with a field survey. Three sets of structural joints (C, J1 and J2) in the source rockfall failure masses are obtained, (C: 191∠14°, J1: 239∠73°, J2: 96∠85°), which can be drawn by a polar stereographic projection of the upper hemisphere (Fig. 7). It can be seen, each set of the dominating discontinuity is ensured by measurement 100 times. The values of dip tendency and dip angle are almost the same with an average error of 3.7° and 4.5°, respectively. According to the combination of rock joints (C, J1 and J2), the boundary, scale and movement direction of potential rockfall can be determined. The spacing of rock joints and volume of rockfall failure masses is listed in Table 2. It can be seen, the volume of rockfall blocks is at the range of 0.5～ 4m³, and the volume of historical rockfall events at the bottom of the slope is concentrated at the range of 0.5～5m³. There is a strong correlation between them.

Model building
The model building includes two parts: the terrain model and the rockfall failure mass model (Fig.  8). The terrain model is obtained by UAV-based aerial photography, and the volume of rockfall failure masses can be obtained by measuring the terrain model. Each rockfall failure mass is formed by a group of rockfall block models arranged according to the characteristics of the structural plane. The rockfall block model includes block shape and size with parameters of density, fixed joint fracture force threshold, dynamic friction coefficient, static friction coefficient and bounce coefficient. The rockfall block model was set at the cube shape for the reason to be close to the shape of the rockfall deposits and the divided rock mass in the rockfall source area as well. As shown in Table 3, d1, d2, and d3 represent the size of the rockfall blocks, which is selected according to the structural plane spacing of rockfall failure mass and the size characteristics of the historical rockfall block mass. The density of limestone is selected as 2.26g/cm 3 based on the uniaxial tests and "Rock Mechanics Parameter Handbook" [24] . The fixed joint fracture force threshold is set as 1000N to be the rock mass shear strength. The dynamic friction coefficient is set to be 0.3, the static friction coefficient is 0.3, and the bounce coefficient is 0.6.

Rockfall movement simulation
Based on the developed rockfall movement simulation software, the Xiejiayan rockfall has been done. Taking the No. 75 rockfall block as an example (Fig. 8), it can be seen that the rockfall is to be unstable due to gravity, then it becomes a free-fall movement. While impacting on the slope surface for the first time, the force exceeds the shear strength then the rock fractured into three independent rock blocks. The No. 75 rockfall block bounced after it collided with the slope surface. Due to the steep slope, it collided with the slope surface four times and bounced to reach the gentle terrain at the bottom of the slope. The No. 75 rock block rolled a certain distance along the slope, and finally stopped under the action of friction, without reaching the passive net (Fig. 9).

Result
In terms of the proposed simulation method, the three-dimensional movement path of each rock can be obtained (Fig. 10). It can be seen, the simulated rockfall deposits are in good agreement with the results from field investigation. The rockfall movement simulation is confined by the physical principle, the kinetic energy of rockfall movement is transformed from the rock mass static energy. When rockfall collides the slope surface, the movement speed will be reduced under the effect of the coefficient of restitution, until the end of the movement. Usually, when the slope is steeper, the boulder has a larger bouncing. The kinetic energy often reaches the maximum when it reaches the middle of the slope, and the maximum bounce height as well (Fig. 11). The Xiejiayan rockfall frequent happens in recent decades of years. So far, the unstable rock mass in the source area still exists, which is prone to become rockfall events in the rainy season and earthquakes. The proposed rockfall simulation method might play a great role in rockfall prevention and mitigation. The potential movement route, the kinetic energy and the bounce height are very important during a support design or risk management for the residents.

Conclusions
This paper proposes a simulation method to obtain rockfall movement features based on the Unity3D engine. The Xiejiayan rockfall is selected as a case study to illustrate its suitability and feasibility to mitigate the rockfall losses. Three remarkable conclusions can be ensured as follows: (1)The terrain model and rockfall failure mass model can be determined according to the UAV-based photography, combined with a field survey. Unity3D and PhysX physics engine is used to develop the simulation software under CPU-GPU high-performance parallel calculation method.
(2)The presented software is used to simulate the movement process in Xiejiayan rockfall. Compared with traditional two-dimensional simulation, the three-dimensional route of the rockfall movement process is more reliable.
(3)The three-dimensional dynamic visual display of the rockfall movement process can be presented, which is much more realistic and intuitive. The simulation not only can help geological disaster management but also promote public risk perception.
For technical reasons, we haven't realized the function of dividing rock and soil with different properties and giving different surface parameters in the software. We will continue to study this problem in-depth, to improve the accuracy of simulation results.