4D monitoring of mountain areas using the UAV-PPK workflow

The mapping and three-dimensional modeling of mountain areas and the study of land cover change dynamics are current tasks in preserving and maintaining protected natural parks and forests. In this context, recent developments in digital photogrammetry using the SfM-MVS method to process captured imagery and the development of unmanned aerial systems (UAS) allow for reducing the costs, time, and the use of human resources and obtaining and repeatable 3D topographic data for moun-tainous regions. We will call this acquired 3D high-resolution topographic data (HRTD) 4D data in the context of an additional temporal component. The main objective of this study is to evaluate the applicability of PPK (Post-Processing Kinematic) direct georeferencing of images captured by UAVs and processed through the SfM-MVS method to obtain HRTDs for 4D land cover analysis. We analyze a 3D HRTD with an acquisition interval of two years for a mountain test area in Plana Mountain near Sofia. The test area has a diverse vegetation cover, including coniferous forest, grassland, hay mead-ows, shrubs, and single deciduous trees. We conducted multiple surveys of the test area with a budget PPK-UAV configuration (DJI Phantom 4 Pro with a single-frequency PPK-GNSS kit installed) from March 2020 to October 2022. Two autumn surveys from September 2020 and October 2022 were selected, which possess the most -good performance on numerical data accuracy. We performed 3D data analysis on 1) Assessment of the accuracy of PPK-SfM-MVS photogrammetry generated topographic data (3D clouds and DSM); 2) Investigation of the errors in the individual specific surfaces (for the individual isolated sections) using the M3C2 tool for comparing and evaluating dense point clouds; 3) Determining land cover changes in the demarcated areas using a surface of differences (DoD).


Introduction
Changes in the landscape and land cover of mountainous regions are caused by several factors: a) natural phenomena (erosion, avalanches, floods, extreme climatic events), b) human intervention (logging, afforestation, construction of roads and facilities, animal husbandry) c) natural change in plant cover (growth and aging of forest massifs and ecosystems).Today's rapid development of remote sensing through unmanned aerial systems (UAS) has increased the quality of high spatial resolution topographic data (HRTD).UAV-captured images are processed by applying a structure-from-motion (SfM) -Multi-View Stereo (MVS) ) photogrammetry workflow to obtain high-resolution topographic data (i.e., point clouds, orthophoto mosaics, digital surface models -DSM or digital elevation models -DTM).This process allows us to obtain a high-quality 3D reconstruction of an object (surface), starting from a series of two-dimensional images captured from different viewpoints (Westoby et al. 2012, Eltner et al. 2016).
HRTD generated by images from unmanned aerial systems (UAS) with an accuracy better than 20 cm opens new possibilities in studying and monitoring natural landscapes.Typically, the accuracy of such topographic data depends on the number and distribution of ground control points (GCPs) by which the digital model is georeferenced.Rapid and efficient mapping of hazardous or inaccessible areas is hampered by the need to place and measure multiple GCPs to georeference the UAV aerial imagery data.Direct georeferencing of UAV imagery using the PPK-GNSS method can overcome these limitations by providing centimeter-level HRTD accuracy for specialized land surface surveys without the need for ground-based GCPs (Zhang et al. 2019).
With the availability of good quality multitemporal HRTDs acquired over a sufficient vegetation interval, there is an opportunity to investigate land surface changes for areas of interest.Several authors in recent years have studied the topographical changes of riverbeds (Wheaton et al. 2010, Westaway et al. 2003), observations and monitoring of glaciers (Dall' Asta et al. 2017), (Ryan et al. 2015), snow cover monitoring (Nolan et al. 2015), landslides and agriculture (Clapuyt et al. 2017), (Mauri et al. 2021), coastal erosion (Brunier et al. 2016), as well as change in forest structural and biomass (Zahawi et al. 2015), (Messinger et al. 2016), (Mlambo et al. 2017), (Gobbi et al. 2022).The application of the direct method for georeferencing UAV images in the research of topographic changes is discussed in the publications of (Turner et al. 2016, Gonzalo 2018, Tomaštík 2019, and Zhang et al. 2019).
This study proposes an integrated multi-temporal land cover change analysis, including a conifer forest plot, a single trees study, uncultivated green areas, a hay meadow, and a low bushes area.High-resolution multi-temporal 3D digital land surface models (DSM) were used to assess land cover dynamics for a test mountain area.Two UAV-PPK surveys were selected for analysis at the test site with an interval of two years.Thus, UAV-based 4D-monitoring of land cover dynamics (e.g., tree canopy) over time represents an additional novel aspect of our work.Finally, the proposed study fills a gap in scientific knowledge regarding the possibility of efficiently analyzing land cover evolution and landform change in mountainous regions using a UAV-PPK-SfM workflow.This workflow starts from unmanned photogrammetric capturing with direct georeferencing of the imagery, proceeds through SfM-MVS photogrammetry for HRTD acquisition, and concludes with land surface change analysis.

Materials and methods
Figure 1 shows the complete UAV-PPK-SfM workflow (modified from Cucchiaro, 2018a) through which the observations, surveys, and processing of the results for 4D monitoring of mountainous areas were carried out.In the following sections, we give full details about: Data collection and PPK georeferencing (sections 1,2,3 in Fig. 1), data processing (sections 4,5,6,7 in Fig. 1), data and error analysis (sections 8 and 9 in Fig. 1), and finally, the evaluation of the topographic DEM difference and visualization of multitemporal profiles (DoD; sections 10 and 11 in Fig. 1).

Study site
The research experiment was conducted in mountainous terrain near the village of Plana, close to the Geodetic Observatory "PLANA" at the Bulgarian Academy of Sciences.The site is near the capital city and offers convenient year-round UAV survey access.The area is hilly terrain with a slope of about 8.3% from north to south, with an average altitude of 1200 m.For our survey, various land cover is  available to test for accurate reproducibility and changes in the surface from multitemporal imagery.Grass vegetation is dominant in the area, there are single trees and plots with shrubs, and there is a coniferous forest at the periphery of the test site (Fig. 2).In the northern part of the test site, the herbaceous vegetation was mowed in 2022 and can be examined for change from the 2020 baseline survey when the grassland was not cropped.

UAVs platforms and PPK-module
For the unmanned photogrammetric survey, a consumer-grade UAV was used -DJI Phantom 4 Pro (P4P) -Figure 3c.The P4P has a DJI FC6310 camera with a nominal focal length of 8.8 mm and a 1" CMOS 20-megapixel sensor with a 2.41 x 2.41 µm nominal pixel size.To implement the PPK-GNSS experiment, we used the platforms equipped with a single-band compact GNSS-RTK receiver (Reach M +, Emlid Ltd) with RTK / PPK capabilities (Fig. 3a) and a Single-band RTK GNSS receiver (Fig. 3b) with centimeter precision (Reach RS+, Emlid Ltd) as a base station.The base station is mounted on a tripod located in the outline of the test area for aerial capturing.The base receiver is configured to log the raw data to a RINEX file at 5 Hz using GPS and GLONASS satellites.The Phantom 4 Pro's front LEDs flash whenever a photo is taken.A phototransistor registers the blink, and the timestamp is written to the Reach m+ log file.Timestamps are extracted and combined with the photos as precise geotags (Dinkov 2020).

Data Acquisition
Two UAV-SfM-PPK surveys were conducted on 19.09.2020 and 05.10.2022.In addition, GCPs stabilized in the test area were used to evaluate the digital products and additional validation points chosen during the acquisition (Fig. 4).
The photogrammetric missions were predefined with the UgCS flight planning software (Fig. 4c).Summary metrics for the flight missions conducted are reported in Table 1

PPK-SfM-MVS workflow with accuracy assessment
The PPK-SfM-MVS workflow was applied to process the two sets of aerial imagery and produce high-resolution digital topographic data (Dinkov 2023), which we also used to assess the accuracy of the final products.The PPK-SfM-MVS process starts with initial data collection -photogrammetrically capturing the area and recording GNSS raw data to apply PPK precision georeferencing to the imagery.The next step is calculating the exact image coordinates after the rover (UAV) and Base Station joint processing of the GNSS records.Next, the SFM method calculates the camera parameters and generates an initial sparse point cloud.This step is followed by the internal accuracy evaluation of the initial model (via independent validation points) and iterative calibration of the camera parameters.If we use GCP for precise model orientation, now is the time to input their coordinates and use them for georeferencing.Once this initial computational stage is complete, processing continues through the MVS workflow to generate High-Resolution Topographic Data (HRTD) 1) Dense Point Cloud, 2) Digital Surface Model (DSM), and 3) Orthophoto Mosaic.An assessment of HRTD accuracy follows.In a previous study (Dinkov, 2022), the workflow for estimating the accuracy of point clouds and numerical surface model is discussed in detail.In the present study, PPK-direct georeferencing of aerial photographs was applied.In this case, GCP and validation points (VP) were used to estimate the spatial accuracy of HRTD without the points being involved in the georeferencing of the digital data.The most reliable estimation that can be made is based on the use of validation control points by measuring the differences between the locations reported in the generated HRTD (dense cloud points, orthophoto maps, and DSM) and the spatial coordinates measured in the field with a GNSS-RTK receiver (ASPRS 2015).The accuracy assessment of point cloud, orthophoto mosaic, and DSM was conducted using 24 validation points for Study-1, 2020, and 19 validation points for Study-2, 2022.For the internal point cloud evaluation (SfM Photogrammetry Workflow), point positions are validated directly in the 3D point cloud in Pix4D software using the built-in RayCloud 3D Editor.The following errors are calculated for the VP in the photogrammetric project: • Mean (μ): The average error in each direction (X, Y, Z) by formula (1) • Sigma (σ): The standard deviation of the error in each direction (X, Y, Z) by formula (2) • RMSEe: The Root Mean Square Error in each direction (X, Y, Z) by formula (3) Model (DSM); and 3. Orthophoto Mosaic.Next, the accuracy of the HRTD is evaluated.Our previous study (Dinkov 2023) discusses in detail the workflow for estimating the accuracy of point clouds and numerical surface models.In the present study, PPK-direct georeferencing of aerial photographs was applied, and in this case, GCP and VP were used to estimate the spatial accuracy of HRTD without the points being involved in the georeferencing of the digital data.The most reliable estimation that can be made is based on the use of validation control points by measuring the differences between the locations reported in the generated HRTD (dense cloud points, orthophoto maps, and DSM) and the spatial coordinates measured in the field with a GNSS-RTK receiver (ASPRS 2015) .The accuracy assessment of point cloud, orthophoto mosaic, and DSM was conducted using 24 validation points for Study-1, 2020, and 19 validation points for Study-2, 2022.For the internal point cloud evaluation (SfM Photogrammetry Workflow), the validation point positions are determined directly in the 3D point cloud in Pix4D software using the built-in RayCloud 3D Editor.The following errors are calculated for the validation points (VP) in the photogrammetric project: • Mean (μ) : The average error in each direction (X,Y,Z) -formula (1) • Sigma (σ): The standard deviation of the error in each direction (X,Y,Z) -formula (2) • RMSE e : The Root Mean Square error in each direction (X,Y,Z).-formula (3) (2) Where Δe is the error of each point for the given direction.It represents the dif-In these formulas, Δe i is the error of each point for the given direction.It represents the differences between the accurately measured coordinates (X, Y, or Z) and the reported points in the dense point cloud.N is the number of tie points.
The Mean Z error helps to recognize systematic errors due to bad GCP acquisition.
The Δe i differences provide a measure of the uncertainty of each point cloud.The mean of the differences μ indicates the accuracy of the registration process; the standard deviation σ of the differences indicates the accuracy level (Cucchiaro et al. 2018a).In addition, the root mean squared error (RMSE 3D ) was calculated in the x, y, and z directions to further check for potential systematic errors in the points (Remondino et al. 2017) using the formula: (4)

Methods for measuring the distance between 3D HRTDs (surfaces, point clouds)
Below we present the main approaches (Lague et al. 2013) used to measure the distance between two surfaces and point clouds in the context of 4D -Land Cover Change Detection.

DEM of Difference (DoD)
DEM of Difference is the most common method for comparing two surfaces in the earth sciences when the large-scale geometry of the scene is planar (i.e., riverbed, cliff).The two surfaces (DEM or DSM) are compared pixel by pixel, which is equivalent to a vertical distance measurement.This technique is very fast and now includes explicitly calculating the uncertainties associated with point cloud registration, data quality, and point cloud roughness (Wheaton et al. 2010).However, the DoD technique suffers from two significant drawbacks (Lague et al. 2013): A) It cannot work correctly in 3D environments with highly dissected terrain, as DEM cannot handle overhanging parts (cliffs and drop-offs of banks, large blocks ), B) The fixed resolution of DEM imposes a limitation on the level of detail retained by the raw data, which can be a substantial limitation for surfaces showing many different characteristic scales.

A minimum level of change detection in DoD
The most commonly adopted procedure for managing DoD uncertainty involves setting a minimum detection threshold (minLoD) to distinguish real surface changes from intrinsic noise (Brasington et al. 2003, Wheaton et al. 2010).
The minimum detection level (minLoD) is calculated as in (Brasington et al. 2003): minLoD=t.√(δzold) 2 + (δz new) 2 ( 5) (2) Where Δe i is the error of each point for the given direction.It represents the differences between the accurately measured coordinates (X,Y or Z) and the reported points in the dense cloud.N is the number of tie points.
The Mean Z error helps to recognize systematic errors due to bad GCP acquisition.The Δe i differences provide a measure of the uncertainty of each point cloud.The mean of the differences μ indicates the accuracy of the registration process; the standard deviation σ of the differences indicates the accuracy level (Cucchiaro et al. 2018a).In addition, the root mean squared error (RMSE 3D ) was calculated in the x, y, and z directions to further check for potential systematic errors in the points (Remondino et al. 2017): In this formula, minLoD is the distributed error in DoD, which represents the minimum level of detection threshold (LoDmin), and δzold and δznew are the individual error in the reference DSM and comparative DSM, respectively, t is the critical t-value at the chosen confidence level (for a conservative approach, a value of t = 1.96 was used, corresponding to a confidential interval of 0.95).

Direct cloud-to-cloud comparison with closest point technique (C2C)
This method is the simplest and fastest method for direct 3D comparison of point clouds, as it does not require gridding or data linking nor computation of surface normals (Lague et al. 2013).For each point in the second point cloud, the closest point in the first point cloud can be determined.This method was developed for fast change detection in very dense point clouds, not for accurate distance measurement (Girardeau-Montaut et al. 2005).

Cloud-to-mesh distance or cloud-to-model distance (C2M)
This approach is the most common technique in inspection software.The surface change is computed using the distance between a point cloud and a 3D grid reference or theoretical model (Olsen et al. 2010).This approach works well on flat surfaces as a mesh can be constructed based on the average positions of the points in the cloud.(Kazhdan et al. 2006)

Algorithm M3C2 (Multiscale Model to Model Cloud Comparison)
The M3C2 algorithm finds the most appropriate normal direction for each point and then calculates the distance between the two point clouds along a cylinder of a given radius projected in the direction of the normal (Lague et al. 2013).In addition to the calculated distance, M3C2 reports the number of points in the projection cylinder (a measure of point density) and the standard deviation of the points in the cylinder (a measure of local roughness).The standard deviations are combined with the registration error to provide error estimates and delineate significant change areas.

Quantification of land cover changes by DoD (digital elevation models to calculate differences)
Land cover changes affecting the test site were analyzed by comparing DSMs, explicitly using the DEM method for Difference (DoD) calculation.DoDs were calculated with the Geomorphic Change Detection (GCD) software (Wheaton et al. 2010).However, due to the morphological complexity of the study area, it was not possible to use uniform error and, therefore, we applied a minimum detection level (minLoD) for the DoD calculation (minLoD) for the DoD calculation (Mauri et al. 2021;Bossi et al. 2015), to achieve efficient results.
The method for estimating the individual level of uncertainty was performed depending on the different surface types in the study area for which the DoD was calculated (Fig. 5) The different land cover types for the sections are 1) Conifer Forest, 2) Rough grassland, 3) Mowed meadow, 4) Single trees, 5) Bushes and low trees.In this study, the uncertainty level of each model was estimated using CloudCompare's M3C2 tool, calculating the cloud-to-cloud distance between each distinct plot.Specifically, the standard deviation of this distance is used as an indicator of the uncertainty associated with the analyzed pair of point clouds, and we adopted it as the uncertainty ("DSM") affecting the corresponding surfaces (Cucchiaro et al. 2018b).
Finally, using the Geomorphic Change Detection (GCD) software (Wheaton et al. 2010), the changes and uncertainties in each DSM and their propagation to the Difference DEM were quantified.

Point clouds accuracy
Here we summarize the achieved point cloud accuracies from the two airborne missions considered, calculated by comparing validation point coordinates.Table-2 shows the number of points in the two computed dense point clouds, the accuracy of each point cloud (described by the absolute mean of the VP differences), and the precision of the point cloud (represented by the standard deviation of the coordinate differences of VP) for the 2020 and 2022 surveys and the overall root mean squared error (RMSE 3D ), in the position of the validation points.

DSMs generation and accuracy
The SfM-MVS image processing was performed in the Pix4Dmapper software, and the generated surfaces of DSM have a resolution of 1.9 cm for both the 2020 and 2022 models, respectively.
The accuracy of the DSM was assessed using the open-source GIS software QGIS, and details of the HRTD accuracy assessment can be found in a previous study by (Dinkov 2023).The planimetric X and Y coordinates of the validation points were output to an external point shapefile, where the points represent the centers of the landmarks recognized in the orthophoto mosaic.Afterward, Z coordinates for the vectorized points were added using an additional tool (plugin) to extract the Z values from the raster file of DSM.The values of the coordinates thus obtained were exported as a text file and compared with the reference positions (measured by GNSS-RTK) of the GCP and VP. Figure 6 shows a 3D visualization of the DSMs with the locations of validation points marked.For each DSM, a line chart with the reported errors is shown.Each curve (with nodal points) shows the deviation between the reference (measured with a high-precision GNSS receiver) and the DSM-reported values in altitude.Table 3 shows the mean squared height error for each DSM.

Comparing point clouds and identifying the minimum threshold for change
The dense SfM point clouds for the two epochs (2020 and 2022) measurements were compared to determine the minimum threshold of change for each test section.The following Figures 7,8,9,10,and 11 show the individual sectors and the results of calculations with the M3C2 plugin (CloudCompare) regarding the calculated mean distances between point clouds and their standard deviations.

Digital elevation models of Difference
The two DSMs obtained by SfM-MVS photogrammetry from Survey-1-2020 and Survey-2-2022 were compared to establish the level of land cover change.The study was done separately for each section with different land cover in the test area.The individual DoD models with corresponding estimates were obtained using Geomorphic Change Detection (GCD) software (Wheaton et al. 2010) to identify vertical and volumetric changes.We assumed that the standard deviation σ of each dataset from the cloud comparison using M3C2 represents the DSM error and minLoD.Therefore, those DoD cells with absolute values below minLoD are considered uncertain and are not used to calculate topographic changes (i.e., threshold DoD).
The topographic changes for the sections with different land cover are shown in Figures 11,12,13,14,and 4.

Discussion
For a general investigation of the feasibility of the UAV-PPK-SfM workflow, we performed a total of 11 studies of the test site "Plana." The experiments were realized with the equipment discussed in this paper for the period March 2020 -October 2022; more details in (Dinkov 2020(Dinkov , 2021)).For the land cover change analysis, we selected two surveys from different epochs for the same season (autumn) with the most robust performance on accuracy and precision of the final high-resolution topographic data.

UAV-PPK-SfM workflow with an evaluation of digital topographic data
Direct PPK georeferencing provided centimeter-level accuracy and precision to our topographic data.Applying the UAV-PPK-SfM workflow using L1-PPK equipment  provides close accuracy to traditional GCP georeferencing techniques established by a previous study (Dinkov 2023).As found in other studies (Zhang et al. 2019, Štroner et al. 2021, Cirillo et al. 2022), direct georeferencing with precise positioning can replace the conventional method using ground control points and allow obtaining of HRT data with centimeter accuracy.
In Study-1, performed in September 2020, there were offsets in the vertical position of PPK GNSS coordinates due to false solutions (e.g., false correction when resolving ambiguities).Implementing one GCP improved our survey results, and the centimeter accuracy of the digital data was achieved (Dinkov 2021).It should be noted that applying a GCP only moves the whole project to the approximate location without internal georeferencing, as shown by (Zhang et al. 2019).
Accuracy assessment of dense point clouds and terrain surfaces (DSM) was performed using validation points surveyed with high-precision GNSS-RTK receivers.This assessment indicates the centimeter accuracy of the digital data illustrated in Table 2.However, these validation points are fixed and marked on the terrain.The accuracy of the SfM point cloud depends strongly on the characteristics of the measured surface (Cook 2017).In mountainous areas, it is essential to understand the effects of different surface types on SfM results.The main parameters and errors in the calculation of SfM point clouds are presented in Table 2, highlighting the accuracy in centimeters of the overall accuracy of the SfM surveys described by the corresponding GCP and VP errors.

Surface error analyses
As noted above, due to the morphological heterogeneity and complexity of the test site, we consider separate specific areas (sections) with a distinct vegetation cover type.For the same reason, we did not apply a single land-cover error criterion.Accordingly, we calculated a specific threshold (minLoD) for detecting land-cover changes for each section.Then, applying the M3C2 tool in CloudCompare, we obtained the mean distances between separated multitemporal cloud points for individual sections with their corresponding standard deviations, as shown in Figs. 7,8,9,10,and 11 are summarized in Table 5.The results show that for the forest section (Section -1), we have the highest value for uncertainty (standard deviation σ = 0.748 m).Vegetation has long been recognized as a source of error in photogrammetry due to leaf clumping caused by wind, and illumination variation, which increases the complexity of the images (Messinger et al. 2016).For this reason, the use of a spatially explicit error threshold for topographic change detection is necessary to improve the registration reliability and quantification of changes over time.Not so high values but above 0.15 m are observed for "Section-4: Single Trees" and "Section-5: Bushes and Tree", while for the sections with herbaceous vegetation cover (Sections -2 and 3), the error is below 0.10 m.
The error analysis demonstrates the applicability of the UAV-PPK-SfM workflow in surveying mountainous areas with varying vegetation cover.However, the accuracy and precision of the SfM-generated point clouds depend strongly on the characteristics of the measured surface (Cook 2017).Therefore, all these errors must be considered in the workflow to obtain more realistic spatial representations of topographic changes (Cucchiaro et al. 2018a).The results of the DoDs allowed the analysis of land cover changes for each delineated plot.The GCD budget partitioning function permits us to evaluate in detail the areas where we report vegetation cover increases and decreases, vertical vegetation changes, and the volumes of vegetation decrease and increase in different areas, plus a representation of the net volume change.This is shown in Figures 12,13,14,15,and 16 by Diagrams of Topographic Change (2).Our data show that over the two years at the test site, there was an increase in the coniferous forest by an average of about 1.80 m and an increase in single deciduous trees by about 1.10 m.The increase and decrease in vegetation cover in the hay meadow, rough grassland, and scrub zones were correctly recorded.To better visualize the changes that have occurred in vegetation cover, we have also applied multitemporal profiles (across the 2020 and 2022 DSM surfaces) along a selected transect line "A-B" for each distinct section.These profiles are shown in Figures 12,13,14,15,and 16 by graphs (3).The profiles correctly show land cover changes that are also reported by DoDs.For example, in Fig. 11(2) for section-1, the bars reporting " Vertical Changes" show an average increase in conifer forest cover of about 1.80 m, which is also evident in the profile plot (Fig. 11(3)).In Fig. 13 for section-3, the total reduction in grass cover due to mowing in 2022 is accounted for by GCD.A similar pattern is also seen in Figure 15 (section-5), where the existing shrubs in 2020 are mowed in 2022.In Fig. 15(3), the flattened section (2022) of existing 2020 shrubs is clearly visible in the multitemporal profile, while the transect through nearby trees accounts for their canopy growth.

Conclusions
This study evaluates the application of a direct PPK method for georeferencing UAV images acquired over a large interval (2 years) to detect changes in vegetation cover in a mountainous region.We analyzed dense 3D point clouds and high-resolution digital surface models obtained from two photogrammetric surveys and processed them by the SfM method.Our study showed that applying PPK in direct georeferencing can provide centimeter-level accuracy and precision, significantly improving field survey efficiency.In addition, we assessed the positional accuracy and repeatability of the DSM by repeating the same flight plans.The PPK solution had similar accuracy to the traditional approach using GCP-based georeferencing.However, some flights (Survey-1, 2020) were characterized by vertical offset, resolved by including one GCP in the workflow.
The applied workflow for 4D UAV-PPK monitoring of an area with hilly terrain and varied land cover demonstrates the relevance and feasibility of the method for obtaining multi-temporal HRTD.The same workflow and hardware can be applied in difficult and inaccessible mountainous areas.This technique enhances field survey efficiency by avoiding the need for a GCP.A particular advantage is that it will increase the accuracy of topographical reconstructions of hard-to-reach areas such as active landslides, volcanic terrain, and glaciers (Clapuyt et al. 2017), where the accessibility of GCP placement is critically limited.
The complex and inhomogeneous structure of mountainous landscapes makes it difficult to adequately pre-assess the quality of digital data acquired by the SfM method (Smith et al. 2016).Therefore, to achieve quality monitoring data for mountainous areas applying the 4D UAV-PPK-SfM workflow, it is necessary to refine both the PPK-georeferencing of the captured images and the uncertainty assessment of the digital data.
In computing changes in vegetation cover in single (detached) trees with no homogeneous surface, the volumes are loaded with error due to forming a cylindrical DSM surface from the canopy to the ground surface.In such cases, more precise ground-based measurement methods (e.g., ground-based laser scanning (Brodu and Lague 2012)) and results processing (direct comparison of 3D point clouds (Lague et al. 2013, Williams et al. 2021) should be applied.
Finally, our test study on 4D land cover monitoring in mountainous areas provides additional value as a basis for further research to integrate new sensors, technologies, and methods to investigate processes concerning vegetation cover status and dynamics in the mountains.

Figure 1 .
Figure 1.General workflow chart presented in this paper is composed of 11 main interrelated sections (modified, Source: Cucchiaro, 2018a).

Figure 2 .
Figure 2. (A) Location of the study area within the borders of Bulgaria (42°28' N; 23°24'E).(B) Location of the test site near the Geodetic Observatory in the Plana mountain.

Figure 4 .
Figure 4. (A) Test area -control and validation points for 2020, (B) Test area -control and validation points for 2022; (C) Prepared flight plan with Phantom 4 Pro in UgCS software.

Figure 5 .
Figure 5. Test area with plotted test sections for comparison of surfaces for the two periods.

Figure 6 .
Figure 6.DSMs of the PLANA test site with signs in validation points: (a) for Survey-1 (2020); (b) for Survey-2 (2022); (c) Line chart of the reported errors in the validation points for Survey-1; (d) Line chart of the reported errors in the validation points for Survey-2 15.For each experimental sector, a numerical difference (DoD) model and land cover change plots are visualized using colored bars (red indicates an increase in vegetation cover, and green indicates a decrease).Diagrams illustrate areas with registered change (m 2 ), vertical change (m), and volumetric (m 3 ) changes.Note that the colored bars represent the mean values, and the error bars show the potential variation (±) associated with the propagated error.Each figure also shows a multitemporal profile along transect A-B on the DSM surfaces for 2020 and 2022.Details of the land cover change results for each section are shown in Table

Figure 8 .
Figure 8. Signed distances (m) between the point clouds of the Survey-1 (2020) and survey-2 (2022) for Section (2) "Rough Grassland": (a) Orthophoto mosaic ; (b) Distance map between dense point clouds (c) Gaussian distribution plot of the M3C2 distance between the two point clouds -mean μ = 0.052 m with a standard deviation σ = 0.054 m.

Figure 9 .
Figure 9. Signed distances (m) between the point clouds of the Survey-1 (2020) and survey-2 (2022) for Section (3) "Mowed Meadow": (a) Orthophoto mosaic ; (b) Distance map between dense point clouds (c) Gaussian distribution plot of the M3C2 distance between the two point clouds -mean μ = 0.053 m with a standard deviation σ = 0.068 m.

Figure 10 .
Figure 10.Signed distances (m) between the point clouds of the Survey-1 (2020) and survey-2 (2022) for Section (SA) "Single Trees": (a) Orthophoto mosaic ; (b) Distance map between dense point clouds (c) Gaussian distribution plot of the M3C2 distance between the two point clouds -mean μ = 0.070 m with a standard deviation σ = 0.173 m.

Figure 11 .
Figure 11.Signed distances (m) between the point clouds of the Survey-1 (2020) and survey-2 (2022) for Section (5) "Bushes and Trees": (a) Orthophoto mosaic ; (b) Distance map between dense point clouds (c) Gaussian distribution plot of the M3C2 distance between the two point clouds -mean μ = 0.001 m with a standard deviation σ = 0.209 m.

Figures
Figures 12, 13, 14, 15, and 16 -Pattern of vegetation cover changes (1) show the DoD (after applying minLoD) for each delineated reach based on surveys conducted in September 2020 and October 2022.The results of the DoDs allowed the analysis of land cover changes for each delineated plot.The GCD budget partitioning function permits us to evaluate in detail the areas where we report vegetation cover increases and decreases, vertical vegetation changes, and the volumes of vegetation decrease and increase in different areas, plus a representation of the net volume change.This is shown inFigures  12, 13, 14, 15, and 16  by Diagrams of Topographic Change (2).Our data show that over the two years at the test site, there was an increase in the coniferous forest by an average of about 1.80 m and an increase in single deciduous trees by about 1.10 m.The increase and decrease in vegetation cover in the hay meadow, rough grassland, and scrub zones were correctly recorded.To better visualize the changes that have occurred in vegetation cover, we have also applied multitemporal profiles (across the 2020 and 2022 DSM surfaces) along a selected transect line "A-B" for each distinct section.These profiles are shown inFigures 12, 13, 14, 15, and 16 by graphs (3).The profiles correctly show land cover changes that are also reported by DoDs.For example, in Fig.11(2) for section-1, the bars reporting " Vertical Changes" show an average increase in conifer forest cover of about 1.80 m, which is also evident in the profile plot (Fig.11(3)).In Fig.13for section-3, the total reduction in grass cover due to mowing in 2022 is accounted for by GCD.A similar pattern is also seen in Figure15(section-5), where the existing shrubs in 2020 are mowed in 2022.In Fig.

Table 1 .
Summary metrics for the flight missions

Table 2 .
Main errors in the UAV-SfM-PPK surveys point clouds.

Table 3 .
Vertical root means square errors (RMSEz) of validation points used in Pix4Dmapper software (in meters).

Table 4 .
Changes in vegetation cover volumes for individual sections.The table shows the DoD threshold and corresponding raw results for each sector (area).

Table 5 .
Cloud-to-cloud distance values in the Sections of the Plana study area