Accurate analysis is required due to the high hazard and risk status of rockfalls. From this aspect, modeling studies have great importance (Mary Vick et al. 2019). For rockfall modelling studies, geomorphometric approaches and modelling based on high-resolution DSM data provide more accurate and sensitive results (Loye et al. 2009, Zhang et al. 2019, Pérez-Rey et al. 2019, Francioni et al. 2020, Rodriguez et al. 2020, Akın et al. 2021). Within this scope, UAVs offer significant advantages in producing DSM and orthophoto data for large areas in a short duration and gain great importance for monitoring the distribution of natural disasters, especially, in spatial and temporal terms (Feng and Röshoff 2004, Abellán et al. 2006, Armesto et al. 2009, Alejano et al. 2013, Giordan et al. 2015). High-resolution DSM and orthophoto data are very effective for accurate and sensitive trajectory modeling and the detection of source zones in rockfall prone areas. For this reason, rockfall studies in recent years have moved beyond traditional methods (Colomina and Molina 2014, Boccardo et al. 2015, Gomez and Purdie 2016, Manconi et al. 2019) and begun to use high-resolution DSM and orthophoto images obtained with UAVs (Matasci et al. 2015, Török et al. 2017, Byrne 2018, Manconi et al. 2019).
3.1. UAV Survey and Photogrammetric processing
In this study, a DJI Phantom 4 Advanced UAV was used to produce high-resolution orthophoto and digital surface models for rockfall modeling. The general flowchart followed when producing this data is given in Fig. 3.
Preparations for flight height, scanning width and stereo image matching proportions were completed with Pix4dcapture software for the study area. With the aim of producing DSM and orthophoto data with good quality and sensitivity, 40 ground control points (GCP) were placed with D-GPS before the flight. Five parallel flights were made at 100 meters height, and 621 stereo aerial photographs were obtained with 70% overlap ratio (Fig. 4a-b). Processing of these images was performed using the Pix4d trial version. Eventually, point cloud data with a 3.03 cm resolution were obtained after the photogrammetric processing of orthophoto images with low (1.4 cm) RMSE error (Fig. 4c-d).
3.2. Evaluation of the quality of elevation data of DSMs
There are some errors when DSM data is obtained using different platforms. These errors are caused by the structure of the land surface and land use/cover features.
One of these errors is the vertical error rate (Ajayi et al. 2017, Coveney and Roberts 2017). Ground control points (GCPs) and checkpoints (CPs) are used to fix these errors during collecting stereo images by UAV. In order to understand the vertical accuracy error rate, CP data is also collected in addition to the GCP data. Determining the vertical accuracy rate in DSM data is important for the performance of the study and the precision of the results (Tamminga et al. 2015).
In this respect, using 35 checkpoints in the study area, the vertical accuracy rate of the DSM data root mean square error (RMSE, Formula 1), mean square error (MSE, Formula 2), mean absolute deviation (MAD, Formula3) and mean absolute tested with percentage error (MAPE, Formula 4) (Tamminga et al. 2015, Akturk and Altunel 2019). These formulas are expressed as,
The results obtained according to the specified formulas are given in Table 1. Accordingly, the error rates of the DSM data are MAD: 0.169, MSE: 0.066, RMSE: 0.258, MAPE: 0.01 cm. Accordingly, DSM data was obtained with a low vertical accuracy error rate.
Table 1
The vertical accuracy results of the DSM data based on different statistical parameters using RMSE, MSE, MAD, MAPE.
Check point
|
Actual Elevations
|
Forecast Elevations
|
Error
|
Absolute Value of Error
|
Square of Error
|
Absolute Values of Errors Divided by Actual Values.
|
|
At
|
Ft
|
At -Ft
|
| At -Ft|
|
( At -Ft)^2
|
| (At -Ft)/At|
|
1
|
1630.25
|
1630.10
|
-0.15
|
0.150
|
0.023
|
0.0001
|
2
|
1638.29
|
1638.13
|
-0.16
|
0.160
|
0.026
|
0.0001
|
3
|
1637.38
|
1638.00
|
0.62
|
0.620
|
0.384
|
0.0004
|
4
|
1628.94
|
1628.87
|
-0.07
|
0.070
|
0.005
|
0.0000
|
5
|
1619.62
|
1619.55
|
-0.07
|
0.070
|
0.005
|
0.0000
|
6
|
1633.70
|
1632.93
|
-0.77
|
0.770
|
0.593
|
0.0005
|
7
|
1596.07
|
1595.99
|
-0.08
|
0.080
|
0.006
|
0.0001
|
8
|
1590.10
|
1590.09
|
-0.01
|
0.010
|
0.000
|
0.0000
|
9
|
1593.89
|
1593.85
|
-0.04
|
0.040
|
0.002
|
0.0000
|
10
|
1589.94
|
1589.85
|
-0.09
|
0.090
|
0.008
|
0.0001
|
11
|
1587.47
|
1588.00
|
0.53
|
0.530
|
0.281
|
0.0003
|
12
|
1593.24
|
1593.13
|
-0.11
|
0.110
|
0.012
|
0.0001
|
13
|
1631.93
|
1631.87
|
-0.06
|
0.060
|
0.004
|
0.0000
|
14
|
1569.82
|
1569.78
|
-0.04
|
0.040
|
0.002
|
0.0000
|
15
|
1612.38
|
1612.30
|
-0.08
|
0.080
|
0.006
|
0.0000
|
16
|
1566.89
|
1566.78
|
-0.11
|
0.110
|
0.012
|
0.0001
|
17
|
1622.20
|
1623.00
|
0.80
|
0.800
|
0.640
|
0.0005
|
18
|
1595.68
|
1595.60
|
-0.08
|
0.080
|
0.006
|
0.0001
|
19
|
1562.14
|
1562.00
|
-0.14
|
0.140
|
0.020
|
0.0001
|
20
|
1574.44
|
1574.32
|
-0.12
|
0.120
|
0.014
|
0.0001
|
21
|
1571.70
|
1571.60
|
-0.10
|
0.100
|
0.010
|
0.0001
|
22
|
1589.67
|
1589.60
|
-0.07
|
0.070
|
0.005
|
0.0000
|
23
|
1566.96
|
1567.10
|
0.14
|
0.140
|
0.020
|
0.0001
|
24
|
1599.51
|
1599.40
|
-0.11
|
0.110
|
0.012
|
0.0001
|
25
|
1612.30
|
1612.20
|
-0.10
|
0.100
|
0.010
|
0.0001
|
26
|
1592.02
|
1591.99
|
-0.03
|
0.030
|
0.001
|
0.0000
|
27
|
1547.38
|
1547.70
|
0.32
|
0.320
|
0.102
|
0.0002
|
28
|
1556.64
|
1556.50
|
-0.14
|
0.140
|
0.020
|
0.0001
|
29
|
1586.68
|
1586.60
|
-0.08
|
0.080
|
0.006
|
0.0001
|
30
|
1576.71
|
1576.50
|
-0.21
|
0.210
|
0.044
|
0.0001
|
31
|
1583.12
|
1583.00
|
-0.12
|
0.120
|
0.014
|
0.0001
|
32
|
1553.27
|
1553.20
|
-0.07
|
0.070
|
0.005
|
0.0000
|
33
|
1548.69
|
1548.60
|
-0.09
|
0.090
|
0.008
|
0.0001
|
34
|
1555.24
|
1555.13
|
-0.11
|
0.110
|
0.012
|
0.0001
|
35
|
1552.49
|
1552.40
|
-0.09
|
0.090
|
0.008
|
0.0001
|
|
Totals
|
|
-1.090
|
5.910
|
2.326
|
0.004
|
|
N
|
35
|
|
|
|
|
|
MAD
|
MSE
|
RMSE
|
MAPE
|
|
|
|
0.169
|
0.066
|
0.258
|
0.01
|
|
|
3.3. Discontinuity analysis based on point cloud data
Discontinuity analysis, which is used in the identification of fracture and crack joints in rock surface, is important to understand behavior of the blocks. Discontinuities are explained by creating contour diagrams depending on the angle of the extension of the cracks and fractures in the rock units with the north (Zhang 2006, Riquelme et al. 2018). In recent years, LIDAR, TLS, and UAV platforms, which are actively used and highly preferred in the creation of digital surface data (DSM-DTM) and play a major role in the creation of discontinuities (Dewez et al. 2016, Riquelme et al. 2017, 2018, Valkaniotis et al. 2018). Discontinuities were performed in 3 part of the high slope source zone in this study (Fig. 5) and discontinuities were analyzed based on UAV-3D Point cloud data using qFacet plugin that allows to extract characteristics of the rock surface such as planar facets and their orientation and distance (Dewez et al. 2016, Valkaniotis et al. 2018). The mean dip angel and dip direction for the study area gives a similar result based on high-density point cloud datasets (Fig. 5. 1a-2b-3c). Besides, the discontinuity results obtained from 3 different locations gave very similar values. Accordingly, the obtained orientations generally refer to the areas corresponding to the vertical cooling cracks due to the lithology of the area. The general values obtained accordingly show that there are 3 main discontinuity surfaces represents Fig. 5b-c-d and the mean dip direction (°) of these surfaces is between 190–208°, on the other hand, the dip angle (°) values are between 33–35° (Fig. 5 and Table 2).
Table 2
The summary of the discontinuities in the study area based on qFacet plugins
Section
|
No. of planes
|
Mean
|
Dip (°)
|
Dip direction (°)
|
RMSE
|
1
|
7602
|
30
|
190
|
0.143
|
2
|
4269
|
35
|
203
|
0.123
|
3
|
2034
|
33
|
208
|
0.124
|
3.4. Block geometries
In the study area, which is at high risk of rockfall, the source zones of blocks generally correspond to cooling cracks located on steep slopes of cliffs (Figure. 2). To identify blocks that have potential falling risks, field studies were performed and 17 blocks were identified (Fig. 6). The height, width, and length were measured for each block and according to these values each rock is modelled in 3D RAMMS: ROCKFALL 1.6. 70 module software (Table 3). While 10 of the modeled rock blocks (R1-R2-R4-R6-R7-R10-R11-R15-R16-R17) have equant and real-equant geometry, 5 rocks (R3-R9-R12-R13-R14) have long and real-long geometry and R5-R8 blocks have flat geometry (Fig. 7). RAMMS software provides a great advantage for accurate and high-quality modeling of rocks with different shapes and volumes measured in field studies (Mary Vick et al. 2019) and takes into account the real size of blocks (Dorren and Kühne 2016, Torsello et al. 2021).
Table 3
Properties of the simulated rocks
No
|
X (length)
m
|
Y (width)
m
|
Z (height)
m
|
Volume
m³
|
Mass
kg
|
Rock shape
|
R1
|
3.2
|
2.5
|
1.5
|
7.96
|
21485.4
|
Equant
|
R2
|
2.0
|
1.8
|
1.0
|
1.545
|
4172.5
|
Equant
|
R3
|
1.5
|
1.6
|
1.2
|
0.719
|
1941.7
|
Real long
|
R4
|
1.42
|
1.23
|
1.0
|
0.95
|
2565.2
|
Real Equant
|
R5
|
1.2
|
1.4
|
0.8
|
0.925
|
2497.8
|
Real Flat
|
R6
|
1.5
|
1.25
|
1.0
|
1.124
|
3035.6
|
Equant
|
R7
|
1.5
|
1.25
|
1.0
|
1.124
|
3035.6
|
Equant
|
R8
|
1.22
|
1.17
|
1.0
|
0.687
|
1856
|
Real Flat
|
R9
|
0.96
|
1.16
|
0.98
|
0.464
|
1253.5
|
Real Long
|
R10
|
1.42
|
1.23
|
1.0
|
0.950
|
2565.2
|
Real Equant
|
R11
|
1.42
|
1.23
|
1.0
|
0.950
|
2565.2
|
Real Equant
|
R12
|
1.15
|
1.49
|
1.0
|
0.698
|
1.885
|
Real Long
|
R13
|
3.24
|
1.62
|
2.23
|
0.524
|
1414.1
|
Real Long
|
R14
|
1.56
|
2.01
|
1.35
|
0.719
|
1941.7
|
Real Long
|
R15
|
1.42
|
1.23
|
1.0
|
0.950
|
2565.2
|
Real Equant
|
R16
|
1.61
|
1.34
|
1.07
|
1.124
|
3035.6
|
Equant
|
R17
|
2.18
|
2.34
|
1.87
|
1.242
|
3352.4
|
Real Equant
|
3.5. Rockfall simulations
A range of rockfall models were developed in order to assess rockfalls in terms of dynamics, trajectories, kinetic energy, velocity and jump height and many other aspects (Azzoni et al. 1995, Dorren 2003, Volkwein et al. 2011, Frattini et al. 2012, Žabota et al. 2021, Liu et al. 2021). RAMMS, CONEFALL, STONE, Georock, Rockyfor3D, FlowR and Rotomap produce two-dimensional (2D) and 3D rockfall models based on different algorithms and GIS for spatial analysis (Guzzetti et al. 2002, Jaboyedoff and Labiouse 2011, Topal et al. 2012, Bartelt et al. 2016a). 3D rockfall simulations are performed on 3D surface/elevation models to determine rockfall trajectories and runout zones (Dorren 2003). Results obtained from this modeling provide advantages in producing hazard maps for rockfall susceptibility, protection, vulnerability, failure and distribution of rockfalls (Vo 2015; Bonneau et al. 2018; Sarro et al. 2018; Singh et al. 2018; Sazid 2019). For this reason, implementing the necessary protective systems (fences, ditches and forests etc.) as a result of modeling by identifying these blocks and performing the necessary measurements may prevent these hazards and risks (Piacentini and Soldati 2008).
In this study, after determining the geometric and volumetric properties of these blocks with very different shapes and volumes, 3D rockfall analyses were executed with the RAMMS: Rockfall software which is a GIS-based rockfall simulation program (Leine et al. 2014; Bartelt et al. 2016a). In the rockfall modeling process, a range of data inputs are required in addition to modeling obtained linked to DSM data. These inputs include features like plant cover (open, dense, etc.), terrain types (hard, extra hard, medium, soft, medium soft, etc.) and block types (flat, equitant, round, etc.) for blocks or source zone areas to be modelled by the RAMMS program (Vo 2015). For rockfall modeling in the RAMMS program, DSM data with resolution from 1–10 m is sufficient. For this reason, DSM data with cm-resolution was rescaled to 1-meter resolution to obtain the format for use in modeling. Considering the volcanic rocks in the field, the terrain type was selected as hard terrain, covered with very sparse plant cover (Fig. 4c). Besides, the friction parameters (Table 4) were determined based on terrain type that account for; “Rocks jump over ground, includes different size of blocks and absence of the vegetation” (Bartelt et al. 2016b).
Table 4
The friction parameters of the simulated rockfalls (Bartelt et al. 2016b).
Material strength
|
Material weakening
|
Rock ejection
|
Material behaviour
|
µ min
|
µ min
|
β
|
κ
|
v
|
0.55
|
2
|
185
|
3
|
0.4
|
3.6. Rockfall Hazard Index
Rockfall Hazard Index (RHI) method was employed for the hazard mapping of rockfalls (Crosta and Agliardi 2003). In this method, the kinetic energy, jump height and block count parameters are taken into consideration which was determined through the RAMMS program in this study. These basic parameters are reclassified within the scope of hazards. During the reclassification, the classification values of Crosta and Agliardi (2003) were used. Necessary classifications for kinetic energy, jump height and block count values are presented in Table 5. For calculation of block count value, the c/(5*n) formula was used and then reclassification was performed. As we can see in the formula, c represents rockfall count, and n represents the number of blocks from each cell value. The values obtained as a result of the modeling of the rock blocks were evaluated separately for each block based on RHI. The values obtained were separately calculated for reach rockfall model with the results of classification of the 3 parameters combined to identify high-moderate-low hazard areas.
Table 5
Classification of parameters according to the Rockfall Hazard Index method (Crosta and Agliardi, 2003).
Class
|
Block count
(local scale)
|
Kinetic energy
(kJ)
|
Jump height
(m)
|
1
|
˂ 0.01
|
≤ 700
|
≤ 4
|
2
|
0.01–0.1
|
700–2500
|
4–10
|
3
|
˃ 0.1
|
≥ 2500
|
≥ 10
|