Accuracy Assessment of UAV-Photogrammetric-Derived Products Using PPK and GCPs in Challenging Terrains: In Search of Optimized Rockfall Mapping

Unmanned aerial photogrammetric surveys are increasingly being used for mapping and studying natural hazards, such as rockfalls. Surveys using unmanned aerial vehicles (UAVs) can be performed in remote, hardly accessible, and dangerous areas, while the photogrammetric-derived products, with high spatial and temporal accuracy, can provide us with detailed information about phenomena under consideration. However, as photogrammetry commonly uses indirect georeferencing through bundle block adjustment (BBA) with ground control points (GCPs), data acquisition in the field is not only time-consuming and labor-intensive, but also extremely dangerous. Therefore, the main goal of this study was to investigate how accurate photogrammetric products can be produced by using BBA without GCPs and auxiliary data, namely using the coordinates X0, Y0 and Z0 of the camera perspective centers computed with PPK (Post-Processing Kinematic). To this end, orthomosaics and digital surface models (DSMs) were produced for three rockfall sites by using images acquired with a DJI Phantom 4 RTK and the two different BBA methods mentioned above (hereafter referred to as BBA_traditional and BBA_PPK). The accuracy of the products, in terms of the Root Mean Square Error (RMSE), was computed by using verification points (VPs). The accuracy of both BBA methods was also assessed. To test the differences between the georeferencing methods, two statistical test were used, namely a paired Student’s t-test, and a non-parametric Wilcoxon signed-rank. The results show that the accuracy of the BBA_PPK is inferior to that of BBA_traditional, with the total RMSE values for the three sites being 0.056, 0.066, and 0.305 m, respectively, compared to 0.019, 0.036 and 0.014 m obtained with BBA_traditional. The accuracies of the BBA methods are reflected in the accuracy of the orthomosaics, whose values for the BBA_PPK are 0.039, 0.043 and 0.157 m, respectively, against 0.029, 0.036 and 0.020 m obtained with the BBA_traditional. Concerning the DSM, those produced with the BBA_PPK method present accuracy values of 0.065, 0.072 and 0.261 m, respectively, against 0.038, 0.060 and 0.030 m obtained with the BBA_traditional. Even though that there are statistically significant differences between the georeferencing methods, one can state that the BBA_PPK presents a viable solution to map dangerous and exposed areas, such as rockfall transit and deposit areas, especially for applications at a regional level.


Introduction
The use of unmanned aerial vehicles (UAVs) for remote sensing has become more common in the study of natural hazards and geomorphological processes, especially due to the development of technologies and mapping systems at larger scales, leading to a more comprehensive understanding of natural processes. UAVs are especially suitable for studying and monitoring rockfalls [1][2][3][4], due to their smaller spatial extent and the fact that the potential rockfall areas are often located on steep slopes above infrastructure and settlements, where the terrain is hardly accessible and dangerous for classical field Kinematics), and PPK (Post-Processing Kinematics) solutions, based on a multi-frequency multi-constellation of GNSS receivers. By using GNSS RTK, positional data are acquired from satellites and a base station providing real-time high-accuracy positions, such that corrections in real time are applied to the data onsite. Although the PPK method gathers positional data by using virtual reference station (VRS) or continuously operating reference stations (CORS) in a similar manner, the corrections are applied in the post-processing stage. With the accuracies of those methods being comparable to that of indirect georeferencing, the GCPs would not be needed, thus decreasing the time required for the acquisition and processing of the data [29,36,37].

Previous Work on UAV Images Georeferencing without GCPs
Several studies have considered using image georeferencing without GCPs in the production of UAV-photogrammetric-derived products. In the following paragraphs, we present a short overview of the results provided only by UAV imagery georeferencing, using BBA without GCPs (BBA_PPK). A complete overview of the reviewed studies is available in Appendix A.
The accuracy of photogrammetric-derived products can be assessed by comparing the photogrammetric models against a reference, or by using verification points (VPs) surveyed with traditional topographical methods, and identified within the BBA and orthomosaics/DSM [51]. In this review, we only included the results for the latter case, as we used the same approach, and focused on horizontal and vertical accuracies that were either compared for BBA, orthomosaics, and DSM. Observing the results for BBA, studies have reported a horizontal accuracy between 0.020 and 0.060 m, and vertical from 0.040 up to 0.300 m. Some studies have also reported lower accuracies: the specifics of these UAV flights were higher flying altitudes of the UAV, flying in only one direction, surveying vegetated areas and areas that are more diverse with extreme elevation changes, and with the CORS station at larger distances (especially affecting the vertical accuracy) [28,31,39,41,48,49]. Horizontal accuracies of orthomosaics have revealed similar results, with values of 0.020-0.050 m, except for Mian et al. (2016) [38] and Tomaštik et al. (2019) [28], where the accuracies resulting from different flight setups and vegetation conditions varied between 0.044 and 0.087 m. Vertical accuracy estimated from the DSM is, in the best case 0.030 m, reaching up to 0.200 m, with some exceptions exceeding these values [37,38].

Aim of the Study
A majority of the presented studies in Section 1.2. were carried out in flat areas with homogeneous terrain conditions (except for the studies by Hugenholtz et al. [36], Tomaštík et al. [28], and Tufarolo et al. [48]), where there is less satellite obstruction, better signal reception, and finally conditions for measuring GCPs and VPs. In our research, we wanted to test the accuracies of the photogrammetric products in the case of rockfalls, as the conditions for a survey of high-accuracy data (e.g., GCP positions) are highly demanding, due to the higher terrain roughness and larger elevation changes. It is crucial for surveyors of these environments to know whether it is possible to achieve results comparable to those acquired by using georeferencing, without requiring any GCP assistance. Therefore, we question which georeferencing method is more appropriate and suitable for providing highaccuracy photogrammetric-derived products, taking into consideration all of described environmental challenges and risks of field operations in the event of rockfall hazards. Either it is permissible to use the georeferencing with PPK (BBA_PPK) method, or the use of indirect georeferencing (BBA_traditional) is still preferable and needed, especially in the case of planning technical protection measures where an accuracy of the order of centimeters is required. As we wish to minimize the amount of time that field operators spend on an active rockfall slope, the main goal of the research was to analyze which BBA method (BBA_traditional or BBA_PPK) provides more accurate photogrammetric-derived products (e.g., orthomosaics and DSMs). In the research, we compared the photogrammetric-derived products that have been produced strictly by georeferencing, using BBA_PPK and indirect georeferencing (using BBA_traditional only), to determine the possible applications in the event of rockfall studies.
The reason that we considered the option of not using BBA_PPK in combination with the BBA_ traditional method is that the surveying of rockfalls can sometimes be too risky to place any GCP below, for example, rock walls in the transit corridors of the falling rocks. However, rockfalls located in remote and hardly accessible areas are still valuable and need to be included in past rockfall databases. In our study, we are focused only on rockfall transit and deposit areas, where we are collecting data on past rockfall events needed for calibration and validation of rockfall models. For rockfall modeling, we require data on the volume, location, shape of individual rock deposits, and surface roughness of the slope, among others, and know that the accuracy of photogrammetric-derived products is crucial for rockfall modeling at different scales.
The remainder of the paper is structured as follows: In the Materials and Methods section, we describe the test site, the experimental setup for comparing the results of georeferencing with two different methods (i.e., AT with solely GCPs and using auxiliary data provided by PPK and no GCPs), the UAV feature and flight plan, data processing, measuring and obtaining the positions of the VPs, and statistical evaluation of the results. In the Results section, we present the differences in the accuracy, using VPs, of the BBA (X, Y, and Z), orthomosaics (X and Y), and DSM (Z), as obtained by the two georeferencing methods mentioned above. In the Discussion section, we explain the significance of the results, draw the main outcomes, and discuss possible future investigations. Due to different applications (e.g., comparison of photogrammetric points clouds and laser scanner point clouds; spatial planning-orthomosaics; and so on) [51], our comparison of the results is carried out separately for BBA, for orthomosaics, and DSM.

Study Sites
In the study, we used three rockfall deposit sites, all located in the Julian Alps in the NW part of Slovenia ( Figure 1). Two of the rockfalls are located in glacial alpine valleys-Trenta (rockfall Kekec) and Krnica (rockfall Krnica)-while the third (rockfall Mangart) is located on the mountain pass, below the mountain Mangart. This part of the Alps is geotectonically part of the Southeastern Alps, the thrust unit of the Julian Alps that consists of sediment rocks (mainly from Lower Triassic and Cretaceous, and predominantly from Upper Triassic): limestone and dolomite, with limestone layers commonly passing into dolomite in vertical and lateral directions [52,53]. Faults in this area have transverse Dinaric direction (NE-SW), while significant faults are also present in the Dinaric direction (NW-SE).
The rockfall deposits at the Kekec site consist of massive and bedded dolomite and limestone (Carnian), at the Krnica site of bedded dolomite and limestone (Anisian-Ladinian; Triassic), while at Mangart Pass, they consist of platy and bedded micritic and Remote Sens. 2021, 13, 3812 5 of 31 crinoid limestone with chert (Malm) from the Jurassic period. The specialty of this area (i.e., Mangart Pass) is that these beds are of restricted dimensions in the Julian Alps and represent tectonically confined rests of unconformably overlying limestones from the Upper Triassic period (massive limestone) [54].
All release areas are characterized by steep, vertical slopes and walls from where the rockfall material was deposited directly below the walls, steeper slopes, or lower down in the valleys. The terrain is rough with great elevation and slope values variations which can be observed in Figures 2 and 3. The rockfall runout area was the smallest for the Kekec rockfall (4800 m 2 ), followed by Mangart (13,000 m 2 ) and Krnica (35,000 m 2 ). At all three rockfalls rock deposits, volumes vary from the size of the rocks on the scree slopes of up to 100-200 m 3 at the Krnica rockfall, up to 30 m 3 at the Mangart rockfall, and up to 70 m 3 at the Kekec rockfall. The conditions of surveying all three rockfalls are, therefore, extremely exposed, as the rocks deposited in the rockfall runout area can still be transported down the slope, and as there is always a possibility of new rockfall events occurring. Even though larger rockfall events in these sites have occurred in the past, the activity of falling rocks continues which can be observed by the colors of the rocks (brighter color of rocks indicate more recent events in comparison to darker deposits that are past events), and based on the vegetation (e.g., rocks covered by moss and grass; pioneer tree and shrub species).
Remote Sens. 2021, 13, x FOR PEER REVIEW 5 of 29 dolomite in vertical and lateral directions [52,53]. Faults in this area have transverse Dinaric direction (NE-SW), while significant faults are also present in the Dinaric direction (NW-SE).
The rockfall deposits at the Kekec site consist of massive and bedded dolomite and limestone (Carnian), at the Krnica site of bedded dolomite and limestone (Anisian-Ladinian; Triassic), while at Mangart Pass, they consist of platy and bedded micritic and crinoid limestone with chert (Malm) from the Jurassic period. The specialty of this area (i.e., Mangart Pass) is that these beds are of restricted dimensions in the Julian Alps and represent tectonically confined rests of unconformably overlying limestones from the Upper Triassic period (massive limestone) [54]. All release areas are characterized by steep, vertical slopes and walls from where the rockfall material was deposited directly below the walls, steeper slopes, or lower down in the valleys. The terrain is rough with great elevation and slope values variations which can be observed in Figures 2 and 3. The rockfall runout area was the smallest for the Kekec rockfall (4800 m 2 ), followed by Mangart (13,000 m 2 ) and Krnica (35,000 m 2 ). At all three rockfalls rock deposits, volumes vary from the size of the rocks on the scree slopes of up to 100-200 m 3 at the Krnica rockfall, up to 30 m 3 at the Mangart rockfall, and up to 70 m 3 at the Kekec rockfall. The conditions of surveying all three rockfalls are, therefore, extremely exposed, as the rocks deposited in the rockfall runout area can still be transported down the slope, and as there is always a possibility of new rockfall events occurring. Even though larger rockfall events in these sites have occurred in the past, the activity of falling rocks continues which can be observed by the colors of the rocks (brighter color of rocks indicate more recent events in comparison to darker deposits that are past events), and based on the vegetation (e.g., rocks covered by moss and grass; pioneer tree and shrub species).

Setting Up Ground Control and Verification Points
At each study site we set out VPs, which were later on used for accuracy comparison. Besides the VPs, we also set GCPs, which were used exclusively for indirect georeferencing and were not used as VPs. Before setting up VPs in the field, we have defined the area of interest at each rockfall by digitalizing the known rockfall runout area, using orthomosaics from previous UAV surveys. With the use of the Fishnet function in ArcGIS Pro 2.7.3   VP and GCP targets were squares with a 2 × 2 chess pattern (two white and two black squares; Figure 4), as recommended by the Pix4Dmapper software [61], attached to the ground by a nail through the center of the target. Where that was not possible (as, on the rockfall runout area, most of the terrain was represented by rocks), the targets were fixed with rocks at their edges. At the Krnica rockfall, targets with dimensions 21 cm × 29.7 cm were used, while the targets at Kekec and Mangart targets had dimensions of 29.7 cm × 42 cm. The targets at Krnica rockfall were smaller than at other two locations, as the terrain was too uneven to position larger targets. As all rockfalls are located within the Triglav National Park, any invasive procedures for placing or making different targets were not allowed.
All points were then surveyed, using a Leica Viva TS12 total station, while the orientation points and backside locations were measured with a Leica Zeno20 with GG04 Smart Antenna with accuracy of 0.010-0.030 m. Measurements of these points were carried out in two different parts of the same day (average location of coordinates was later calculated), with each coordinate being measured for an interval of 5 min (with a 6-h of difference). The differential correction data were acquired by using a real-time connection with the Slovenian real-time positioning service SIGNAL. The closest station to all rockfalls was station Bovec, which is located 13 km from Mangart rockfall, 15 km from Kekec rockfall, and 21 km from Krnica rockfall. The X and Y measurements were transformed to the

Setting Up Ground Control and Verification Points
At each study site we set out VPs, which were later on used for accuracy comparison. Besides the VPs, we also set GCPs, which were used exclusively for indirect georeferencing and were not used as VPs. Before setting up VPs in the field, we have defined the area of interest at each rockfall by digitalizing the known rockfall runout area, using orthomosaics from previous UAV surveys. With the use of the Fishnet function in ArcGIS Pro 2.7.3 [55], we created a VP sample network (locations) for each rockfall. The final numbers of VPs were 42 at Kekec, 51 at Krnica, and 48 at Mangart rockfall. GCP locations were also preplanned and set up such that they were not at the same locations as VPs. The number and locations of the GCPs were set up according to the recommendations in the literature; the main condition was that the points were, as much as possible, regularly spaced within the area of interest and its edges. The horizontal distance between the points was approximately 30 m at Kekec, and 65 m at Krnica and Mangart [56][57][58][59][60]. We used 14 GCPs at Kekec, and 16 GCPs at Krnica and Mangart.
The known VP and GCP locations were then imported as shapefiles to the Leica Zeno20 with GG04 Smart Antenna. The points in the field were, thus, set up by using this antenna and the known locations of the VPs and GCPs. As the automatic generation of the VP sample network and predicted locations of GCPs cannot account for the conditions in the field, such that the points cannot be set due to particular terrain features (e.g., a more or less flat area, vegetation, unstable terrain, etc.), the final locations were adjusted when setting the points, based on the predicted locations, while also considering the outline of the rockfall runout area (Figure 3). VP and GCP targets were squares with a 2 × 2 chess pattern (two white and two black squares; Figure 4), as recommended by the Pix4Dmapper software [61], attached to the ground by a nail through the center of the target. Where that was not possible (as, on the rockfall runout area, most of the terrain was represented by rocks), the targets were fixed with rocks at their edges. At the Krnica rockfall, targets with dimensions 21 cm × 29.7 cm were used, while the targets at Kekec and Mangart targets had dimensions of 29.7 cm × 42 cm. The targets at Krnica rockfall were smaller than at other two locations, as the terrain was too uneven to position larger targets. As all rockfalls are located within the Triglav National Park, any invasive procedures for placing or making different targets were not allowed. horizontal geodetic datum "Slovenia 1996/Slovene National grid with Slovenia Geodetic Datum 1996", whilst the Z coordinate was converted to an orthometric height, using the geopotential model EGM96.

Flight Planning and Acquisition of Data with UAV
A DJI Phantom 4 RTK [34] drone with GNSS RTK/PPK technology was used to survey the study sites. The UAV weighs 1391 g, has a diagonal distance of 350 mm, and maximum flight time of 30 min. A camera that is mounted on the UAV has a 1-inch, 20-megapixel CMOS sensor with lens FOV of 84° 8.8 mm/24 mm (35 mm format equivalent) and f/2.8-f/11 auto-focus at 1-∞ m. The ISO range for photos is 100-3200 for auto mode, and 100-12,800 for manual. The maximum image size is 5472 × 3648 (3:2), with a pixel size of approximately 2.4 μm. Multi-Frequency Multi-System High-Precision RTK GNSS enables GPS (L1/L2), GLONASS (L1/L2), and Galileo (E1/E5a) tracking. The manufacturer states that it enables a positional accuracy (RMSE) of 0.015 m (vertical) and 0.010 m (horizontal), both with 1 ppm, which means that the error has a 1 mm increase for each 1 km of movement by the UAV.
The UAV surveys were performed on three separate days for each location: Krnica rockfall was surveyed on 29 August 2019, Mangart rockfall on 20 September 2019, and Kekec rockfall on 7 November 2019. UAV flights were planned such that the flying height was kept constant, using LiDAR-derived digital terrain model (DTM) with a spatial resolution of 1 m. This feature allows for the acquisition of images with the same scale, despite the terrain being very rugged with changing elevation (Figure 3). The flight pattern was perpendicular to the surface slope. The flying height was set to 80 m at all study sites, to achieve a ground sampling distance of 2-3 cm. The image overlaps were 80%. Pictures were taken at 45° and 90° angles, as the combination of oblique and nadir images contributes to a more accurate representation of terrain [62]. For planning the UAV mission, we used the DJI GS RTK application.
The surveyed area for each site is provided in Table 1. We surveyed all sites in PPK mode, such that GNSS data were recorded during the flight and post-processed by using correctional data from the base station that was located within or next (<1 km) to the surveyed area. The location of the base station was calculated with RTK VRS method based on the Slovenian SIGNAL service (interval observations were performed with time spacing of more than one hour and a half). The cadence of observations was one observation All points were then surveyed, using a Leica Viva TS12 total station, while the orientation points and backside locations were measured with a Leica Zeno20 with GG04 Smart Antenna with accuracy of 0.010-0.030 m. Measurements of these points were carried out in two different parts of the same day (average location of coordinates was later calculated), with each coordinate being measured for an interval of 5 min (with a 6-h of difference). The differential correction data were acquired by using a real-time connection with the Slovenian real-time positioning service SIGNAL. The closest station to all rockfalls was station Bovec, which is located 13 km from Mangart rockfall, 15 km from Kekec rockfall, and 21 km from Krnica rockfall. The X and Y measurements were transformed to the horizontal geodetic datum "Slovenia 1996/Slovene National grid with Slovenia Geodetic Datum 1996", whilst the Z coordinate was converted to an orthometric height, using the geopotential model EGM96.

Flight Planning and Acquisition of Data with UAV
A DJI Phantom 4 RTK [34] drone with GNSS RTK/PPK technology was used to survey the study sites. The UAV weighs 1391 g, has a diagonal distance of 350 mm, and maximum flight time of 30 min. A camera that is mounted on the UAV has a 1-inch, 20-megapixel CMOS sensor with lens FOV of 84 • 8.8 mm/24 mm (35 mm format equivalent) and f/2.8-f/11 auto-focus at 1-∞ m. The ISO range for photos is 100-3200 for auto mode, and 100-12,800 for manual. The maximum image size is 5472 × 3648 ( both with 1 ppm, which means that the error has a 1 mm increase for each 1 km of movement by the UAV. The UAV surveys were performed on three separate days for each location: Krnica rockfall was surveyed on 29 August 2019, Mangart rockfall on 20 September 2019, and Kekec rockfall on 7 November 2019. UAV flights were planned such that the flying height was kept constant, using LiDAR-derived digital terrain model (DTM) with a spatial resolution of 1 m. This feature allows for the acquisition of images with the same scale, despite the terrain being very rugged with changing elevation (Figure 3). The flight pattern was perpendicular to the surface slope. The flying height was set to 80 m at all study sites, to achieve a ground sampling distance of 2-3 cm. The image overlaps were 80%. Pictures were taken at 45 • and 90 • angles, as the combination of oblique and nadir images contributes to a more accurate representation of terrain [62]. For planning the UAV mission, we used the DJI GS RTK application.
The surveyed area for each site is provided in Table 1. We surveyed all sites in PPK mode, such that GNSS data were recorded during the flight and post-processed by using correctional data from the base station that was located within or next (<1 km) to the surveyed area. The location of the base station was calculated with RTK VRS method based on the Slovenian SIGNAL service (interval observations were performed with time spacing of more than one hour and a half). The cadence of observations was one observation per second, while the number of station satellites varied between 4 and 7. GNSS data from the flights were post-processed, using Leica GeoOffice and RTKlib.

Data Processing
The number of images used for each site are presented in Table 1. The images were processed with Pix4Dmapper (Version 4.6.4) by Pix4D [61]. For each site, we carried out two separate processing procedures, following the georeferencing methods: (a) processing by carrying out BBA_traditional and (b) processing by carrying out BBA_PPK. For the latter, a geolocation file was used to add the coordinates to the perspective centers of the images, while, for the former, the geolocation data of images were removed. The coordinate system used for the horizontal geodetic datum was Slovenia 1996/Slovene National grid with Slovenian Geodetic Datum 1996, while that for the vertical was the EGM96 Geoid model. The configuration of processing parameters was the same under all georeferencing methods: initial processing was performed in full mode; for point cloud densification, we used an image scale of 1/2 (i.e., half image size) and high point density; and the additional photogrammetric products were orthomosaic and DSM ( Figure 5).
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 29 method of inverse distance weighting. The results of photogrammetric processing are provided in Figure 5. The presented products were all produced by using the BBA_traditional method.

Accuracy Assessment of Bundle Block Adjustment, Orthomosaic and Digital Surface Model
The accuracies of the BBA_traditional and BBA_PPK, namely of BBA (X, Y and Z), orthomosaic (X and Y) and DSM (Z) were assessed with the help of the VPs. The coordinates of VPs, as computed by BBA_traditional and BBA_PPK, were retrieved through Pix4Dmapper, namely by using Checkpoints. They were marked in the same manner as GCPs, i.e., by measuring the center of the target (in this study, we marked each VP at 30 images). Then, the difference between the initial and computed position of the checkpoints (for X, Y and Z) was calculated and shown in the quality report [48]. The calculated differences were further used in comparative analyses.
Coordinates of the VPs of the orthomosaics (X and Y) and DSMs (Z) were retrieved by using GIS software, namely ArcGIS Pro 2.7.3 [55]. The VP targets were vectorized based on the orthomosaic. To accurately determine the center of the target, the procedure of vectorizing was carried out in three steps: (i) In the first step, a polygon with the actual size of target was made, orientated in the same direction as the target, and centralized according to the center of the target, such that it was covered by a polygon. (ii) In the second step, two polylines were drawn from each polygon's angle, such that they crossed in the middle of the target. Lastly, (iii) in the third step, the point that represents the center of the target was placed at the crossing of the polylines. The location of the point was used The precision values of the perspective centers of coordinates for BBA_PPK method are given in Table 2, while they were not needed for BBA_traditional, as all positional data of the images were removed before the processing.  Table 1). For the creation of the DSM, we used default options in Pix4DMapper, namely noise filtering to correct the altitude of points with the median altitude of neighboring points, as well as sharp surface smoothing. The generation of the point cloud can lead to noisy and erroneous points, and with the noise filtering, the altitude of these points was corrected with the median altitude of the neighboring points. Once the noise filter was applied, a surface was generated by using these points. However, the surface can still contain areas with erroneous small bumps, and with surface smoothing these areas are then corrected by flattening. Sharp surface smoothening is by default, and this type of smoothening tries to preserve the orientation of the surface, and keeps sharp objects. In this process, only quasi-planar areas were flattened. The DSM was generated in the raster GeoTIFF file format, using the method of inverse distance weighting. The results of photogrammetric processing are provided in Figure 5.

Accuracy Assessment of Bundle Block Adjustment, Orthomosaic and Digital Surface Model
The accuracies of the BBA_traditional and BBA_PPK, namely of BBA (X, Y and Z), orthomosaic (X and Y) and DSM (Z) were assessed with the help of the VPs. The coordinates of VPs, as computed by BBA_traditional and BBA_PPK, were retrieved through Pix4Dmapper, namely by using Checkpoints. They were marked in the same manner as GCPs, i.e., by measuring the center of the target (in this study, we marked each VP at 30 images). Then, the difference between the initial and computed position of the checkpoints (for X, Y and Z) was calculated and shown in the quality report [48]. The calculated differences were further used in comparative analyses.
Coordinates of the VPs of the orthomosaics (X and Y) and DSMs (Z) were retrieved by using GIS software, namely ArcGIS Pro 2.7.3 [55]. The VP targets were vectorized based on the orthomosaic. To accurately determine the center of the target, the procedure of vectorizing was carried out in three steps: (i) In the first step, a polygon with the actual size of target was made, orientated in the same direction as the target, and centralized according to the center of the target, such that it was covered by a polygon. (ii) In the second step, two polylines were drawn from each polygon's angle, such that they crossed in the middle of the target. Lastly, (iii) in the third step, the point that represents the center of the target was placed at the crossing of the polylines. The location of the point was used to calculate the X and Y location from the orthomosaic. The whole procedure was performed three times, separately [37], and, in the end, an average value of three coordinates was used as the final one, which was used in further analyses. The Z coordinate was extracted from the DSM layer, based on the planimetric (XY) locations of VPs as measured in the field, namely by using the function Extract Values to Points.
The coordinate differences, as mentioned above, were then used to compute several statistics. These concern the minimum, maximum, median, mean, standard deviation, and root mean square error (RMSE) of the differences between the coordinates measured on photogrammetric-derived products (X M , Y M and Z M ) and given (X VP , Y VP and Z VP ) coordinates of the VPs (Equations (1)-(3)) [63]: The RMSE, given by Equation (4) for the differences concerning each coordinate (C) separately (n is a total number of coordinates), was computed as an accuracy value in planimetry (Equation (5), [63]), and a total accuracy value (Equation (6), [63]) was also computed as an accuracy value in planimetry (Equations (4)-(6); [63]): The computed differences (diff X , diff Y and diff Z ) are shown in the form of histograms for each photogrammetric product and georeferencing method separately (Section 3.1, Section 3.2, Section 3.3). Additionally, we also plotted the differences spatially to study any correlation with the terrain topography; the differences were plotted for each photogrammetric product and georeferencing method (Section 3.4).
To compare the accuracy of both the AT, using the BBA_traditional and BBA_PPK (X, Y and Z), and the produced orthomosaics (X and Y) and DSM (Z) for each of the three study areas, two statistical tests were carried out. Whereas the Shapiro-Wilk test [64] indicated that the difference values were normally distributed (p-value > 0.05), a paired Student's t-test [65] was used, otherwise a non-parametric Wilcoxon signed-rank test [66] was used. All calculations and statistical analyses were carried out in RStudio [67]. Table 3 shows the several statistics computed to assess the accuracies of BBA_traditional and BBA_PPK. On the average, the BBA_traditional presented the smallest minimum value of the differences per coordinate. The differences range from −0.098 m (Krnica, diff X ) to 0.091 m (Krnica, diff Y ) in the case of BBA_traditional, while the differences with BBA_PPK were higher; they range from −0.171 m (Mangart, diff Y ) to 0.431 m (Mangart, diff Z ). The mean value of differences ranges from −0.006 m (Kekec, diff Y ) to 0.009 m (Kekec, diff Z ) in the case of BBA_traditional. The mean value of the differences with BBA_PPK were higher, and range from −0.109 m (Mangart, diff Y ) to 0.245 m (Mangart, diff Z ). The standard deviations of the differences were similar for both BBA methods for all coordinates when observing the Kekec (from 0.008 to 0.010 m) and Krnica site (from 0.014 to 0.031 m). At Mangart site, the BBA_PPK had larger standard deviations than BBA_traditional, namely reaching 0.028 for diff X , 0.047 m for diff Y , and 0.092 m for diff Z , while BBA_traditional had values between 0.007 and 0.008 m.  The overall accuracy of BBA (RMSE XYZ ) was higher in the case of BBA_traditional (range from 0.014 to 0.036 m); namely, in the case of the BBA_PPK method, the RMSE XYZ ranged from 0.056 to 0.305 m. Figure 6 shows the frequency distributions of the diff X , diff Y , and diff Z at the sites for both georeferencing methods. It can be observed that the differences related to BBA_traditional method are closely clustered around 0, having uni-modal distribution. The differences related to the BBA_traditional coincided with BBA_PPK, with the latter having a larger span of values, with the exception of the Mangart location, where the differences of BBA_traditional and BBA_PPK coincided in a small part or not at all.  A paired Student's t-test, for testing whether there was a statistically significance difference in the mean value of differences between the BBA_traditional and BBA_PPK, was used in three cases: diffX, diffY, and diffZ in Mangart. In all cases, there was a statistically significant difference between the mean values (p ≤ 0.05). For the rest of the pairs, a nonparametric Wilcoxon signed-rank test was used, which showed that, only in the case of diffY at Krnica location, there was not a statistically significant difference between the median values between the georeferencing methods (p > 0.05). A paired Student's t-test, for testing whether there was a statistically significance difference in the mean value of differences between the BBA_traditional and BBA_PPK, was used in three cases: diff X , diff Y , and diff Z in Mangart. In all cases, there was a statistically significant difference between the mean values (p ≤ 0.05). For the rest of the pairs, a non-parametric Wilcoxon signed-rank test was used, which showed that, only in the case of diff Y at Krnica location, there was not a statistically significant difference between the median values between the georeferencing methods (p > 0.05). Table 4 shows the statistics, as detailed in Section 2.5, computed by measuring the VP at orthomosaics produced for the three study areas, and produced by carrying out the AT with the BBA_traditional and BBA_PPK methods. Similarly, as for BBA accuracy, the BBA_traditional method achieved the smallest difference in average in each coordinate, compared to the VPs. The differences range from −0.104 m (Krnica, diff X ) to 0.078 m (Krnica, diff Y ) in the case of BBA_traditional, while the differences with BBA_PPK were higher; they range from −0.181 m (Mangart, diff Y ) to −0.016 m (Mangart, diff X ). The mean value of differences ranges from −0.017 m (Kekec, diff Y ) to 0.018 m (Kekec, diff X ) in the case of BBA_traditional. The mean value of differences with BBA_PPK were higher, and range from −0. In Figure 7, the frequency distributions of the diff X and diff Y for both georeferencing methods and all sites can be observed. Similar as with the assessment of BBA, the lowest values of differences were achieved by BBA_traditional, which had the densest distribution in all cases, and differences are clustered around 0. BBA_traditional's differences did have an overlap with BBA_PPK, but the differences of BBA_PPK experienced larger variability. The location that stood out is Mangart, where the differences of X and Y retrieved from orthomosaic produced by BBA_PPK were larger and do mostly did not match with those of BBA_traditional. The overall accuracy of orthomosaic (RMSE XY ) was higher in the case of BBA_traditional (range from 0.020 to 0.036 m), while the RMSE XY of orthomosaics produced by the BBA_PPK ranged from 0.039 to 0.157 m.

Assessing Accuracy of Orthomosaic
Remote Sens. 2021, 13, x FOR PEER REVIEW 15 of 29 the following cases: diffX and diffY in Kekec, and diffX in Mangart. At all locations and coordinates, there were significance differences between the mean values (p ≤ 0.05). With other pairs, a non-parametric Wilcoxon signed-rank test was used, which also showed that there were statistically significant differences between the median values between the georeferencing methods in all of the cases (p ≤ 0.05).

Figure 7.
Frequency distribution of the X and Y coordinate differences at all three rockfall sites for orthomosaic_traditional (produced by BBA_traditional) and orthomosaic_PPK (produced by BBA_PPK). Bar width is set to 0.01 m.

Assessing Accuracy of DSM
The accuracy of DSMs was assessed through the differences between the VP coordinate (Z) and coordinate extracted from the BBA_traditional and BBA_PPK methods. Table  5 shows the basic results for difference statistics for Z coordinate and the final RMSEZ values. The differences range from −0.079 m (Krnica) to 0.165 m (Krnica) in the case of BBA_traditional, while the differences with BBA_PPK were higher; they range from −0.129 m (Krnica) to 0.453 m (Mangart). The mean value of differences ranges from −0.022 Figure 7. Frequency distribution of the X and Y coordinate differences at all three rockfall sites for orthomosaic_traditional (produced by BBA_traditional) and orthomosaic_PPK (produced by BBA_PPK). Bar width is set to 0.01 m.

LEGEND
diff( X/Y )-MIN the minimum value of the differences in X/Y diff( X/Y )-MAX the maximum value of the differences in X/Y diff( X/Y )-MEAN the mean value of the differences in X/Y diff( X/Y )-MEDIAN the median value of the differences in X/Y diff( X/Y )-SD the standard deviation of the differences in X/Y RMSE( X/Y ) root mean square error of X/Y A paired Student's t-test, for testing if there was a statistically significance difference in the mean value of differences between the BBA_traditional and BBA_PPK, was used in the following cases: diff X and diff Y in Kekec, and diff X in Mangart. At all locations and coordinates, there were significance differences between the mean values (p ≤ 0.05). With other pairs, a non-parametric Wilcoxon signed-rank test was used, which also showed that there were statistically significant differences between the median values between the georeferencing methods in all of the cases (p ≤ 0.05).

Assessing Accuracy of DSM
The accuracy of DSMs was assessed through the differences between the VP coordinate (Z) and coordinate extracted from the BBA_traditional and BBA_PPK methods. Table 5 shows  LEGEND diff Z -MIN the minimum value of the differences in Z diff Z -MAX the maximum value of the differences in Z diff Z -MEAN the mean value of the differences in Z diff Z -MEDIAN the median value of the differences in Z diff Z -SD the standard deviation of the differences in Z RMSE Z root mean square error of Z RMSE Z was generally smaller in the case of the BBA_traditional method; however, in comparison to the BBA_PPK method for the Kekec and Krnica, the difference was only 0.027 m and 0.012 m, respectively. Larger differences in RMSE Z were presented for the Mangart rockfall, namely they reached 0.231 m. The frequency distribution of the Z coordinates ( Figure 8) indicates that DSM produced with BBA_traditional achieved the lowest differences, while the values had larger span, compared to the assessment of BBA alone. The differences of BBA_traditional and BBA_PPK coincided in a larger part of the frequency distribution, with the latter achieving larger differences. The results were, similar as with the assessment of BBA and orthomosaic, different at the Mangart location; where, for the DSM produced with BBA_PPK, the differences were larger and had greater span than those ones produced with BBA_traditional. Moreover, the differences did not coincide at all at the Mangart location.
To test for statistically significant differences in the mean value of differences between the BBA_traditional and BBA_PPK, a paired Student's t-test was used for diff Z at the Krnica and Mangart site. It showed that there were statistically significant differences in the mean To test for statistically significant differences in the mean value of differences between the BBA_traditional and BBA_PPK, a paired Student's t-test was used for diffZ at the Krnica and Mangart site. It showed that there were statistically significant differences in the mean value between the georeferencing methods values (p ≤ 0.05). A non-parametric Wilcoxon signed-rank test was used for diffZ at Kekec, and it showed that there was a statistically significant difference between the median values of georeferencing methods (p ≤ 0.05).

Spatial Distribution of X, Y and Z Differences
To observe the spatial pattern of the differences for individual coordinates, they were plotted for each VP and for all photogrammetric products produced (i.e., BBA, orthomosaics, and DSM). All plots (Figures 9 and 10) include the marked outlines of the rockfall areas. Figure 9 shows the spatial distribution of the difference values diffX, diffY and diffZ of respectively BBA_traditional and BBA_PPK methods. When observing the results of BBA_traditional, those differences appear to be randomly distributed through the rockfall areas, except for Krnica where larger diffX and diffY are concentrated in the central part

Spatial Distribution of X, Y and Z Differences
To observe the spatial pattern of the differences for individual coordinates, they were plotted for each VP and for all photogrammetric products produced (i.e., BBA, orthomosaics, and DSM). All plots (Figures 9 and 10) include the marked outlines of the rockfall areas.  The results for the orthomosaics produced by BBA_traditional (Figure 10, X and Y coordinates) showed similarity to the BBA accuracy assessment for X and Y coordinates. The spatial pattern and the difference degree were both similar for all rockfalls. For the orthomosaics produced by the BBA_PPK method, no spatial pattern was recognized for Kekec and Krnica, as they varied across the whole rockfall sites. At the Mangart location, larger differences can be recognized towards the northwestern part of the rockfall site.
As for the DSM result for BBA_traditional ( Figure 10, Z coordinate), Krnica rockfall was the only one that had larger differences in the central part of the rockfall while, for the other two, the spatial pattern of the differences could not be recognized. On the other Figure 9. Spatial distribution of the differences in X, Y and Z coordinates, calculated for accuracy assessment of BBA produced by the BBA_traditional and the BBA_PPK method. Figure 9 shows the spatial distribution of the difference values diffX, diffY and diffZ of respectively BBA_traditional and BBA_PPK methods. When observing the results of BBA_traditional, those differences appear to be randomly distributed through the rockfall areas, except for Krnica where larger diffX and diffY are concentrated in the central part of the rockfall area. BBA_PPK results show that at Krnica diffZ expressed larger differences in the southeastern part of the rockfall site, whereas, at the Mangart site, there was a systematic increase of diffX and diffY in towards the northwestern part and towards the southwestern part of the rockfall for the diffZ, which also had the largest values.
Remote Sens. 2021, 13, x FOR PEER REVIEW 18 of 29 BBA_PPK methods, when assessing the accuracy of DSMs showed that, at the Kekec site, the differences were larger in the southern direction; at the Krnica site, they were larger in the southeastern direction; and, at the Mangart site, they were larger in the southwestern direction. The largest differences were observed for the Mangart site. Figure 10. Spatial distribution of the differences in X and Y coordinates, calculated for accuracy assessment of orthomosaic produced by the BBA_traditional and the BBA_PPK method, and of the differences in Z coordinates, calculated for accuracy assessment of DSM produced by the BBA_traditional method and the BBA_PPK method.

Discussion
Georeferencing of UAV-acquired images has been studied by various authors, as already summarized in the Introduction section (see the references therein) and in Appendix A. The vertical and horizontal accuracies vary depending on several factors, including the BBA methods used, namely BBA with and without GCPs. The RMSE reported by the majority of studies varies between 0.020 and 0.060 m for horizontal accuracy, and, in the case of vertical accuracy, 0.020-0.090 m, while studies have also reported values of approximately 0.100 m up to three decimeters (in some cases, even higher).
The best accuracies in our study were achieved by the BBA_traditional method, when comparing the accuracies of BBA, orthomosaic, and DSM. The accuracies for the BBA_traditional method were, for either BBA, orthomosaic, or DSM, in the worst case, under 0.030 m for X and Y coordinates, and under 0.060 m for the Z coordinate. The BBA_PPK method generally achieved lower accuracy than the BBA_traditional method for all photogrammetric-derived products, but it must be highlighted that they were still under 0.035 m for X and Y coordinates and 0.075 m for the Z coordinate (except at the Mangart site), and, comparing the BBA_traditional with BBA_PPK, the results differed, in the worst case scenario, by 0.013 m for X and Y coordinate, and 0.036 m for the Z coordinate. The following results are in similar ranges as those reported in studies carried out by the authors mentioned in the Introduction [25,27,32,33,36,40,[43][44][45][46][47]. In the case of BBA_traditional, the smallest variations (diff-SD) are present at Kekec and Mangart, achieving maximum SD up to 0.010 m (X, Y, or Z coordinate) when observing assessment of the BBA. With orthomosaics the variations were larger for BBA_traditional and in some cases comparable with Figure 10. Spatial distribution of the differences in X and Y coordinates, calculated for accuracy assessment of orthomosaic produced by the BBA_traditional and the BBA_PPK method, and of the differences in Z coordinates, calculated for accuracy assessment of DSM produced by the BBA_traditional method and the BBA_PPK method.
The results for the orthomosaics produced by BBA_traditional (Figure 10, X and Y coordinates) showed similarity to the BBA accuracy assessment for X and Y coordinates. The spatial pattern and the difference degree were both similar for all rockfalls. For the orthomosaics produced by the BBA_PPK method, no spatial pattern was recognized for Kekec and Krnica, as they varied across the whole rockfall sites. At the Mangart location, larger differences can be recognized towards the northwestern part of the rockfall site.
As for the DSM result for BBA_traditional ( Figure 10, Z coordinate), Krnica rockfall was the only one that had larger differences in the central part of the rockfall while, for the other two, the spatial pattern of the differences could not be recognized. On the other hand, a systematic increase of differences for the DSM produced by BBA_PPK is present for Mangart, where the differences had larger values in the southwestern part of the rockfall site. The Kekec rockfall had larger differences in the central part of the rockfall, while differences are larger when moving in the southeastern direction at Krnica.
Observing the spatial distribution of the differences in X, Y and Z coordinate between BBA_traditional and BBA_PPK, when assessing the accuracy of BBA, it is possible to recognize that larger differences were present at the Mangart location. For the X and Y coordinates, the differences were similar across the rockfall area but, for the Z coordinate, they increase in the southwestern direction. Kekec and Krnica had homogeneous differences across the rockfall area with only Krnica having larger differences in Z coordinate in the southeastern part of the rockfall area. Similar conclusions can be drawn for the accuracy assessment of orthomosaics, when comparing the differences in X and Y coordinates of VPs as obtained by using the BBA_traditional and BBA_PPK methods. The difference in BBA accuracy assessment is with the Mangart site, where differences were increasing in northwestern direction for both with X and Y coordinates. No spatial pattern could be recognized for the Kekec and Krnica rockfalls. The spatial distribution of the differences between the Z coordinate of the VPs as obtained by using the BBA_traditional and BBA_PPK methods, when assessing the accuracy of DSMs showed that, at the Kekec site, the differences were larger in the southern direction; at the Krnica site, they were larger in the southeastern direction; and, at the Mangart site, they were larger in the southwestern direction. The largest differences were observed for the Mangart site.

Discussion
Georeferencing of UAV-acquired images has been studied by various authors, as already summarized in the Introduction section (see the references therein) and in Appendix A. The vertical and horizontal accuracies vary depending on several factors, including the BBA methods used, namely BBA with and without GCPs. The RMSE reported by the majority of studies varies between 0.020 and 0.060 m for horizontal accuracy, and, in the case of vertical accuracy, 0.020-0.090 m, while studies have also reported values of approximately 0.100 m up to three decimeters (in some cases, even higher).
The best accuracies in our study were achieved by the BBA_traditional method, when comparing the accuracies of BBA, orthomosaic, and DSM. The accuracies for the BBA_traditional method were, for either BBA, orthomosaic, or DSM, in the worst case, under 0.030 m for X and Y coordinates, and under 0.060 m for the Z coordinate. The BBA_PPK method generally achieved lower accuracy than the BBA_traditional method for all photogrammetric-derived products, but it must be highlighted that they were still under 0.035 m for X and Y coordinates and 0.075 m for the Z coordinate (except at the Mangart site), and, comparing the BBA_traditional with BBA_PPK, the results differed, in the worst case scenario, by 0.013 m for X and Y coordinate, and 0.036 m for the Z coordinate. The following results are in similar ranges as those reported in studies carried out by the authors mentioned in the Introduction [25,27,32,33,36,40,[43][44][45][46][47]. In the case of BBA_traditional, the smallest variations (diff-SD) are present at Kekec and Mangart, achieving maximum SD up to 0.010 m (X, Y, or Z coordinate) when observing assessment of the BBA. With orthomosaics the variations were larger for BBA_traditional and in some cases comparable with BBA_PPK (e.g., Kekec Y coordinate), while maximum differences were present at BBA_PPK. On the other hand, the Z coordinate always achieved the lowest diff-SD in the case of BBA_traditional, including the diff-MAX values, where the differences exceed 0.006-0.008 cm in the case of Kekec and Krnica, and 0.041 cm in the case of Mangart, when compared to BBA_PPK.
The Mangart study site was an exception compared to the other two locations, where RMSE Z for the BBA_PPK method was larger by 0.093 m (X coordinate) and 0.111 m (Y coordinate), and 0.254 m for the Z coordinate, when observing the accuracy of BBA. The results for the orthophoto and DSM accuracies were in a similar range: 0.077-0.114 m for the X and Y coordinate and 0.231 m for Z coordinate. The major cause of such a large difference in accuracy was the unfavorable satellite position configuration [8,27,29] on the day of the UAV survey (more or less three visible satellite positions in the time of survey), compared to that on the day when the GCP/VP survey was performed. Consequently, it is recommended to predict the satellite position quality [29] before the UAV flight, both from the perspective of the final accuracy of the photogrammetric products and for the safety of the flight. Due to the extent of the UAV survey (rapidly changing weather conditions; clouds and fog), and the GNSS survey of the points, both measurements could not be performed on the same day.
Even though our study was carried out in mountainous areas with larger elevation changes along the area of interest (rockfall deposit areas), the final accuracies of BBA, orthomosaics, and DSM were in line with the values reported in other studies. The reliability of the results is therefore high, as similar results have been reported for two rockfall locations that we included in the study; however, the third location highlights that the BBA_PPK method still may-particularly in remote areas-cause larger deviations, and GCPs should be used to validate the results.
The spatial distribution of the differences between the X, Y and Z coordinates of VPs, as obtained by using BBA_traditional and BBA_PPK methods, did not reveal major spatial patterns, while outliers were located randomly across the rockfall sites. The only location that showed a distinctive spatial distribution was Mangart, and that only for the photogrammetric-derived products produced by BBA_PPK. The differences increase towards the northwestern direction for the X and Y coordinates and towards the southwestern direction for the Z coordinate. The spatial pattern was similar with all photogrammetricderived products. The spatial distribution of the differences is a result of the lower accuracy of the BBA_PPK method, and it is not correlated with terrain topography.
The highest differences in coordinates (X, Y and Z) of the VPs when comparing the VP and GCP locations, were achieved at the Krnica rockfall, which was also the largest area that we surveyed. In our experience, measuring the locations of GCPs and VPs can be very challenging. Under the steep rock walls, there is no signal, meaning that the points cannot be measured with GNSS that has RTK; in such cases, the use of the total station is then needed. Additional challenge presented at this study site was moving the measurement equipment, and in the end measurement of the GCPs/VPs. Namely, not all GCPs/VPs could be measured with the total station from one standing position, since the area was too large. Consequently, the total station was positioned on several locations, and due to the large area that needed to be surveyed, not all measurements could be performed in one day. Moreover, due to the high surface roughness, there were many obstacles (rock deposits that are a few meters of high) that obstructed the measurements with total station, and the GCPs/VPs were in some cases located on scree slope (as there was no other option), meaning that they were not attached firmly to the ground. In latter case, targets were more unstable, and could be moved by simple slope processes or simply by walking past it. Consequently, all this factors could potentially influence the measurements, and result in higher errors at this rockfall site. The RMSE values of the BBA_traditional method were larger at the Krnica rockfall, compared to the other two locations, and comparable with the results of the BBA_PPK method. This was partially confirmed by a non-parametric Wilcoxon signed-rank test that showed that there were no statistically significant differences between the georeferencing methods for Y coordinate when assessing the accuracy of the BBA. However, in the assessment of other photogrammetric-derived products and coordinates, the differences were statistically significant between the georeferencing methods.
Taking into consideration this factor, BBA_PPK definitely represents a viable possibility and solution for application in such challenging terrains. Still, for generating products with high accuracy (if possible)-especially vertical accuracy-BBA_PPK should be used at least in combination with few GCPs and possible interferences or obstructions that can influence the final results. Benassi et al. [32] have reported that, by using at least one GCP point, the accuracy can be almost as good as with the use of GCPs for horizontal accuracy, and only slightly worse for vertical accuracy. The fact is that the accuracy can be significantly enhanced with the use of a UAV with PPK, especially when the GCPs cannot be properly distributed and measured in the hazardous areas [27,29,43,47,59,67]. To improve the vertical accuracy, Štroner et al. [47] proposed the use of duplicate flights (double-grid) and geometrically different combinations (different altitudes and camera angle of images); similar conclusions were drawn by Wiącek and Pyka [33]. Harwin et al. [68] concluded that including oblique images is advised, as they can improve self-calibration of the bundle block adjustment and increase horizontal and vertical accuracy. Similar findings were presented by Kyriou et al. [62]; they stated that the combination of nadir and oblique imagery can effectively be used for geomorphological mapping in areas with complex topography and steep slope (>60 • ), to increase the geometric accuracy of UAV data.
Accuracies of PPK solution up to 0.100 cm should provide a satisfactory photogrammetricderived product at a representational scale [28,31,33,39]. Teppati Losè et al. [31] mentioned the use of BBA approaches without GCPs at the scale of 1:500. Przybilla and Bäumker [69] proposed to use RTK/PPK, depending on the application intended for the results. For topography mapping, where the accuracy range is between 0.010 and 0.020 m, they proposed the use of RTK/PPK; with cadaster mapping having an accuracy range between 0.010 and 0.030 m, they proposed the use of RTK/PPK with the GCPs; and, for engineering surveying (accuracy range 0.010 m), the use of GCPs was recommended. Padró et al. [37] stated that the PPK method is consistent enough in generating a large-scale mapping with less effort than the GCP method, and that highly detailed maps (spatial resolution <0.05 m) can be achieved. Tomaštík et al. [28] concluded that the BBA_PPK method can provide an optimal solution for mapping in remote areas.
Based on the results of our study, we can come to similar conclusions. Regarding the use of BBA_traditional or BBA_PPK, the decision regarding the georeferencing method must follow the accuracy requirements of data application and the scale needed, and the conditions for a safely conducting of the field survey. In the event of extreme rockfall events, where rockfall activity is high and the conditions for field work are not safe, we would recommend using the BBA_PPK method only. When the results will be used at a regional level (decimeter accuracy), such as for rockfall mapping of release areas, or the modeling of potential rockfall propagation and deposit areas, the use of the BBA_PPK method could prove to be satisfactory (achieving horizontal accuracies around 0.040 m and vertical up to 0.080 m). On the other hand, when the results of the rockfall modeling are to be used for the planning of the technical protection measures, additional GCPs should be included, to increase the accuracy of the results (up to 0.010-0.020 m) and their reliability. To find an optimal ratio between the accuracy of the photogrammetric products, the work load related to georeferencing and, most importantly, the safety of the field operations, we will continue our research by finding an optimal combination of BBA_traditional and BBA_PPK methods. The following research will, therefore, investigate the number of GCPs needed besides the BBA_PPK method, and the importance of the location where they are situated, especially taking into consideration that they are not located on the exposed and hazardous areas.

Conclusions
The use of UAV platforms equipped with GNSS receivers that provide georeferencing without the use of GCPs is desirable when mapping large and dangerous areas, such as active rockfall areas, as there is no need to use GCPs, as with indirect georeferencing. In this study, we only tested the use of only the georeferencing of the bundle block adjustment (BBA) with PPK method, in comparison to the use of georeferencing of the BBA with GCPs, to see if it is possible to produce high-accuracy photogrammetric products in the case of rockfalls. The main motivation of the study was to improve the condition for field surveying of GCPs, which is labor-intensive, time-consuming, and, most importantly, dangerous. In the majority of cases, the accuracy of the X, Y and Z coordinates, in the case of the georeferencing of the BBA with the PPK method, did not exceed 0.050 m; in the case of orthomosaic, the accuracy of the X and Y coordinates did not exceed 0.034 m; and in case of DSM, the accuracy of the Z coordinate did not exceed 0.072 m. At one location, the results of georeferencing of the BBA with the PPK method were not comparable to those of georeferencing of the BBA with GCPs for all photogrammetric products, which was mostly due to the satellite configuration on the day of the UAV survey with the BBA_PPK method.
We can conclude that the georeferencing of the BBA with the PPK method can provide high-accuracy products, presenting a viable alternative to the georeferencing of the BBA with GCPs. The decision regarding the use of the georeferencing method must be in line with the purpose of the final results; for example, in the regional modeling of rockfalls, at decimeter accuracy, the BBA georeferencing with the PPK method may be satisfactory, while the planning of technical protection measures with rockfalls, at centimeter accuracy, may require the use of GCPs. The PPK method presents a safer option for mapping in mountainous areas. To improve the accuracy of the BBA georeferencing with the PPK