The Methodological Aspects of Constructing a High-Resolution DEM of Large Territories Using Low-Cost UAVs on the Example of the Sarycum Aeolian Complex, Dagestan, Russia

Unmanned aerial vehicles (UAV) have long been well established as a reliable way to construct highly accurate, up-to-date digital elevation models (DEM). However, the territories which were modeled by the results of UAV surveys can be characterized as very local. This paper presents the results of surveying the Sarycum area of the Dagestan Nature Reserve of Russia with an area of 15 sq. km using a DJI Phantom 4 UAV, as well as the methodological recommendations for conducting work on such a large territory. As a result of this work, a DEM with 0.5 m resolution as well as an ultrahigh resolution orthophotoplane were obtained for the first time for this territory, which make it possible to assess the dynamics of aeolian processes at a qualitatively different level.


Introduction
Digital elevation models (DEM) are the main source of information about relief, making it possible to perform geomorphological surveys in laboratory conditions. DEMs are necessary for performing calculations of such parameters as slope, profile and longitudinal curvature of the relief, terrain roughness index and length-slope factor (LS-factor), on the basis of which estimates are made according to erosion models, to predict the development of landslides and other exogenous processes. Different DEMs are used for different study scales: thus, for large-scale studies, researchers usually use DEMs obtained with topographic and geodetic tools, i.e., ground and airborne laser scanners [1][2][3]; for medium-and small-scale studies, global DEMs (GDEMs) of different detail are used [4,5]. Many works are focused on assessing the usefulness of GDEMs as a source of information. For example, in 2012, researchers from Indonesia compared the 2nd and 1st versions of the Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER) mission's GDEM and the Shuttle Radar Topography Mission (SRTM) GDEM, topographic maps at a scale of 1:25,000. Altitude estimates were obtained using RTK-dGPS (real time kinematic-Differential Global Positioning System) receivers. The results of the study showed that the second version of the ASTER DEM has much higher altitude errors, despite the increased resolution, than other DEM sources, as well as the first version of the ASTER DEM [6]. Similar work was carried out by American researchers in 2016: for comparison, versions 1-3 of the ASTER model, National Elevation Dataset (NED) model and SRTM model were taken. As a reference, a network of geodetic control points of the US National Geodetic Service was taken. The method of GPS reference points showed that the NED model has the lowest maximum, minimum and average errors (average error −0.29 m); the maximum errors are typical for ASTER models of all versions [7].
To assess the SRTM (X-SAR) model for applicability in hydrological modeling, a group of German scientists carried out work comparing the SRTM with the regional DGM25 DEM obtained by stereophotogrammetric interpretation of aerial photography [8]. Analysis of SRTM errors at several key sites showed average errors varying in the range receiver survey itself. When comparing the accuracy of the models obtained using TLS and photogrammetry on the reference stand, the error in determining the volumes of the latter was less than 2% [25]. For erosion bench work in the field performed in 2016 by American researchers, the average error was 0.03 cm with a standard deviation of 5 mm [26]. Similar errors were achieved in the same year by British researchers [27]. Thus, recent studies make it possible to conclude that if the methodology is followed, the use of UAVs can enable researchers to survey and create a DEM with similar accuracy to TLS on any territory, from the microrelief of agricultural land, to glacial landforms [28][29][30][31].
As noted above, previous studies have shown the accuracy of the DEM based on the UAV and, in fact, the Structure-from-Motion (SfM) algorithm [32]. Nevertheless, it is important to note the common factors, the presence of which ensures the achievement of high modeling accuracy. Thus, one of the factors is the use of ground control points [33][34][35][36][37]. As noted in Turner et al. [38], the use of ground control points led to a 6-10-fold increase in model accuracy relative to works in which ground control points were not used. The second factor, which is noted by researchers, is the use of the high-precision GNSS receiver, namely the receiver that makes it possible to receive network corrections in real-time kinematic mode (RTK) and L1 and L2 frequency signals to obtain the coordinates of ground control points. The newest GNSS-receivers of the largest manufacturers of geodesic equipment make it possible to reach the positioning error not exceeding 3-5 mm in height (Trimble R8, Trimble R10, LEICA GS18T, Javad Triumph-1, Sokkia GRX2), which is enough for tasks of microrelief reconstruction of lands used in agriculture [39].
The type of UAV used to build the DEM is also important. Globally, there are two types of UAVs-fixed wing and copters. Now there is the so-called hybrid type, which combines the vertical takeoff capability of multi-rotor copters with the stability and flight time of fixed-wing UAVs. The latter, in addition to the above described advantages, have one major disadvantage-a pronounced shutter effect associated with high average flight speed and affecting the planar displacements and overall "smearing" of the obtained models and orthophotomaps [40]. Multirotor UAVs are devoid of such disadvantages; however, they are very limited in the area covered. Currently, the most popular UAV for surveying is the DJI Phantom [41][42][43][44][45][46], which can directly perform surveying work for 15 min. The popularity of DJI Phantom is due to the low cost of the device, availability of components (propellers, servo-driven cameras, additional batteries), and good support from software manufacturers (Pix4d Capture, DroneDeploy, MapPilot, etc.). However, the operating time from a single battery is not enough to carry out work in an area of more than 1 sq. km. The existing works indicate that the methodology of small UAVs surveying large territories more than 10 sq. km is poorly developed [47].
The purpose of this study is to develop methodological recommendations for the survey of complex objects in a large area using low-cost UAVs. The above-mentioned factors determine the relevance, as well as the possible limitations, of the work performed.

Study Area
The study area was chosen as the key site located in the Sarycum area of the Dagestan State Nature Reserve, Republic of Dagestan, Russia ( Figure 1). The site "Sarycum barchans" is located in Kumtorkalinsky district, 18 km north-west of Makhachkala city, at the base of the northern slopes of the Narat-Tyube ridge, on the left bank of the Shura-Ozen river ( Figure 2). It takes the southwestern part of the Sarycum massif. The most striking feature is several large aeolian landforms, with a relative elevation difference of more than 100 m. mountain slopes, shrubby sandy steppe at the foot of barchans, shiblyak bushes on the slopes of the Narat-Tyube ridge, remains of pine forests and juniper woodlands, stony slopes and cliffs. The key area selected for the survey is represented by the Sarykum barkhan itself ( Figure 3) as well as the adjoining territory. The area selected for the survey was 15 sq.km. The amplitude of heights is 197 m from the minimum height in the delta of the Shura-Ozen river-46 m above sea level-to the maximum height of the barkhan-243 m.   According to earlier studies, the absolute altitude of the Sarycum aeolian complex is 243 m [48]. At the same time, the blowing pits are filling up. The open sands of Sarycum in some areas are intensively vegetated. The protected area and its adjacent territories represent an arid space on the border of lowlands and frontier ridges of Dagestan, a unique open sandy desert. Here, due to its physical and geographical specifics, a wide variety of biotopes have formed: a large scattered sandy massif, tree and shrub thickets at its foot, a river valley with meadow vegetation and thinned strips of floodplain tree and shrub vegetation, rocky dry mountain slopes with xerophytic shrubs, clay hilly submountain plain with sagebrush-grass vegetation, fescue-grass steppes on placors and gentle mountain slopes, shrubby sandy steppe at the foot of barchans, shiblyak bushes on the slopes of the Narat-Tyube ridge, remains of pine forests and juniper woodlands, stony slopes and cliffs. The key area selected for the survey is represented by the Sarykum barkhan itself ( Figure 3) as well as the adjoining territory. The area selected for the survey was 15 sq.km. The amplitude of heights is 197 m from the minimum height in the delta of the Shura-Ozen river-46 m above sea level-to the maximum height of the barkhan-243 m.

Field Research Methodology
The main amount of work related to the field collection of information was carried out between June 10 and 16, 2017. Until 2019, flights of UAVs up to 30 kg carried out over the territory of the Russian Federation did not require mandatory approval; however, to ensure airspace safety, altitude restrictions were adopted for the work, equal to 150 m above the point of takeoff. A DJI Phantom 4 UAV was used for surveying the study area (Table 1). A special methodology was developed for surveying such a complex object in terms of planar and altitude characteristics. Since the amplitude of heights as well as the area coverage does not allow capturing the whole object from one point in several flights, Drones 2021, 5, 7 6 of 17 this fact determined the necessity of moving the takeoff points. In order to provide the same illumination conditions, it is recommended to shoot such extended objects not by longitudinal or transverse tack throughout the object, but by small squares. Taking into account the maximum altitude above the take-off point, as well as the operating time of the UAV from one battery in the survey mode, the size of the corresponding "survey squares" is 750 by 750 m. The entire study area was divided into 26 such squares ( Figure 4), which were surveyed sequentially. The placement of the squares was based on the need to ensure 70% longitudinal and transverse overlap of the outermost images of the adjacent squares. In fact, this means that the borders of the squares were docked to each other, with a reserve making it possible to achieve the necessary overlap, even taking into account the changes in the height of the takeoff point.  The process of organizing flights in the field was performed according to the known scheme [11]. Capturing was carried out in an automated mode: the control of the flight, the set overlap of images, the height of the flight, the shutter response time of the camera was carried out in the mobile software Pix4D Capture. Takeoff was performed within the perimeter of the corresponding square. Control of the UAV operator's location within the The process of organizing flights in the field was performed according to the known scheme [11]. Capturing was carried out in an automated mode: the control of the flight, the set overlap of images, the height of the flight, the shutter response time of the camera was carried out in the mobile software Pix4D Capture. Takeoff was performed within the perimeter of the corresponding square. Control of the UAV operator's location within the square was performed using the RTK-GNSS receiver Trimble Geoexplorer 6000, using the same receiver the positioning of ground control points (GCP), was performed. The basic configuration of Trimble GeoExplorer 6000 includes a built-in GNSS-receiver and a sensitive antenna that can receive and process up to 200 channels of GPS and GLONASS at L1 C/A, L2C, L2E and L1 CT, L1 VT, L2 CT, L2 VT frequencies, respectively. The accuracy of coordinate measurement using an external antenna reaches 4 cm in the field (real-time kinematics (RTK)) and 1 cm in the postprocessing mode (postprocessing kinematics (PPK)). In our case, only the field RTK data from the internal antenna was used, obtained for a 70 × 70 cm black-and-white panel. GCPs of this size are clearly visible on images from an altitude of 150 m. A total 18 ground control points were used throughout the study area, 11 of which were used for model georeferencing and 7 for validation of positioning accuracy. As personal experience and numerous articles [49][50][51][52] have shown, the positioning accuracy of UAV survey results does not increase when 15 or more points are used for grounding, regardless of the size of the study area. The selection of positioning points was carried out on the basis of their location on representative elements of the relief-on the most accessible lowest area of the study area, the highest point of the ridge of a barchan, in points with a steep difference of heights on the upper and lower parts of these areas. Such differences are common for the steep slopes of the cliff rock, into which the Shura-Ozen River flows. In the northern part of the study area, since it is more or less flat, the positioning points were set on mottled areas with outcrops of buried sands, since such areas are particularly difficult to align in photogrammetric software, and we tried to compensate low contrast of the textured areas by the presence of a reference point. The locations of the geopositioning and verification points are presented in Figure 5.  The field results were processed using Agisoft Metashape software. This software has a user-friendly interface with rich functionality and the ability to process large data sets using NVIDIA CUDA cores. The field results were processed using Agisoft Metashape software. This software has a user-friendly interface with rich functionality and the ability to process large data sets using NVIDIA CUDA cores.

Results and Discussion
A total of 6,400 images, taken between June 10 and 16, 2017, were used to align the UAV flight results. Eleven GCPs were used for georeferencing. The RMS error of georeferencing of survey results was 0.15 m in X, 0.14 m in Y and 0.16 m in Z, which is an acceptable error according to the State Standard for conducting aerial survey works (GOST R 58854-2020) ( Table 2). To control the accuracy of geopositioning, 7 ground points located in different parts of the study area were used (Table 3). The analysis of the control points proved that the geopositioning of the flight results was performed within the acceptable range for this kind of work's level of accuracy-15-20 cm horizontally and 15-25 cm vertically, according to the national standard in the field of creating models based on the results of airborne survey, valid on the territory of the Russian Federation. The achieved accuracy is explained by the fact that the internal GNSS antenna of the Trimble Geoexplorer 6000 receiver was used for ground control point surveying, using which an accuracy of 0.1 m can be achieved. Taking into account the specifics of the object of research and the above-mentioned permissible accuracies for this kind of work, these limitations can be considered acceptable.
The processing of the results obtained when working on such a large area is not fundamentally different from working with local areas, but to achieve better results and reduce the probability of edge effects associated with the wide-angle lens of the DJI Phantom 4 camera, it is necessary to expand the common point search area. To do this, it is Drones 2021, 5, 7 9 of 17 necessary to set a knowingly coarser camera angle and position accuracy by increasing the camera angle by up to 50 degrees.
When dealing with such a large amount of photographic material, it is not so much the SfM algorithm embedded in a particular software, as the amount of RAM and, perhaps most importantly, the availability of a discrete graphics processing unit (GPU), that comes to the fore. Modern software is capable of using GPUs for computing both in OpenCL and CUDA computing modes. However, based on the results of tests run on systems with NVIDIA (CUDA) and AMD (OpenCL) acceleration units, it should be noted that modern software still works better with NVIDIA GPUs.
The main results of the UAV survey of the area are a digital elevation model and an orthophotoplane. The DEM creation, in our case, was performed by triangulation with linear interpolation over a dense point cloud. The point spacing of the cloud was 0.25 m with a total number of points equal to 124 million. The final resolution of the DEM under interpolation was 0.5 m (Figure 6). It should be noted that the topographic maps for this territory, allowing estimation of the relief of the territory with such detail, do not exist; the most detailed maps allow reconstruction of the relief with a resolution of 1 m at the time of 1963, when the Soviet topographic maps were clarified.  As noted earlier, the DEM was created by interpolating dense point cloud by triangulation with linear interpolation. The choice of interpolation method is based on the performed statistical analysis of residuals. For this purpose, based on the pre-sparse dense point cloud, the Kriging, Inverse Distance to a Power, Minimum Curvature, Nearest Neighbor, Radial Basis Function and, of course, Triangulation with Linear Interpolation methods were used.
Kriging is a geostatistical gridding method that has proven useful and popular in many fields. This method produces visually appealing maps from irregularly spaced data. Kriging attempts to express trends suggested in data, so that, for example, high points might be connected along a ridge rather than isolated by bull's-eye type contours.
With Inverse Distance to a Power (IDP), data are weighted during interpolation such that the influence of one point relative to another declines with distance from the grid node. Weighting is assigned to data through the use of a weighting power that controls As noted earlier, the DEM was created by interpolating dense point cloud by triangulation with linear interpolation. The choice of interpolation method is based on the performed statistical analysis of residuals. For this purpose, based on the pre-sparse dense point cloud, the Kriging, Inverse Distance to a Power, Minimum Curvature, Nearest Neighbor, Radial Basis Function and, of course, Triangulation with Linear Interpolation methods were used.
Kriging is a geostatistical gridding method that has proven useful and popular in many fields. This method produces visually appealing maps from irregularly spaced data.
Kriging attempts to express trends suggested in data, so that, for example, high points might be connected along a ridge rather than isolated by bull's-eye type contours.
With Inverse Distance to a Power (IDP), data are weighted during interpolation such that the influence of one point relative to another declines with distance from the grid node. Weighting is assigned to data through the use of a weighting power that controls how the weighting factors drop off as distance from a grid node increases. The greater the weighting power, the less effect points far from the grid node have during interpolation. As the power increases, the grid node value approaches the value of the nearest point. For a smaller power, the weights are more evenly distributed among the neighboring data points.
Minimum Curvature (MC) is widely used in the earth sciences. The interpolated surface generated by Minimum Curvature is analogous to a thin, linearly elastic plate passing through each of the data values with a minimum amount of bending. Minimum Curvature generates the smoothest possible surface while attempting to honor data as closely as possible.
The Nearest Neighbor (NN) gridding method assigns the value of the nearest point to each grid node. This method is useful when data are already evenly spaced, but need to be converted to a grid file. Alternatively, in cases where the data are nearly on a grid with only a few missing values, this method is effective for filling in the holes in the data. Since the point cloud was pre-generated with the same distance between neighboring points, this method could theoretically perform well in interpolation.
Radial Basis Function (RBF) interpolation is a diverse group of data interpolation methods. In terms of the ability to fit your data and to produce a smooth surface, the Multiquadric method is considered by many to be the best. All of the Radial Basis Function methods are exact interpolators, so they attempt to honor data. The operator can introduce a shaping factor to all the methods in an attempt to produce a smoother surface.
The Triangulation with Linear Interpolation (TIN) method uses the optimal Delaunay triangulation. The algorithm creates triangles by drawing lines between data points. The original points are connected in such a way that no triangle edges are intersected by other triangles. The result is a patchwork of triangular faces over the extent of the grid. This method is an exact interpolator.
In each case, 100 points were randomly selected from the initial point cloud, which did not participate in the interpolation and to which a comparison in specific nodes of the interpolated DEM was made. Two parameters, mean error and root mean square error, were used as basic statistics (Table 4). In all cases, dependence of the estimated value on the true value was plotted (Figure 7).
Despite the fact that the mean error for the Kriging method is minimal, the RMSE is smaller for the TIN method with, in general, a comparable mean error. Another important fact is that Kriging interpolation, regardless of the software in which it is performed, is resource-intensive and takes much more time. For example, the final Kriging method interpolation of dense point clouds took more than three hours on a computing unit powered by 8-core AMD Ryzen 7 CPU. For comparison, interpolation by the TIN method under similar conditions took 3 min. Thus, it is possible to conclude that this method is the most preferable at DEM creation on dense point clouds for large territories. maximum value of difference, average, standard deviation, value of 25 (Q1), 50 (median) and 75% percentiles (Q3). The normality of error distribution was also analyzed by plotting a histogram of frequencies ( Figure 8). The results of statistical analysis of the values of the difference of two DEMs in the selected points are presented in Table 5. Analyzing the obtained results, it can be concluded that the errors have a normal distribution, although the errors tend to move to negative values. In general, it can be stated that the model obtained using the UAV quite accurately describes the study area. It is necessary to note that there are no topographic maps that allow estimation of the relief of the territory with levels of detail that are comparable with the results of the survey from the UAV conducted on this area. The most detailed maps allow reconstruction of the relief with a resolution of 10 m at the time of 1963, when the Soviet topographic maps were being refined. This fact makes it difficult to compare the results of the survey with alternative sources of information about the topography of this territory. In addition, the dynamic relief in this study area creates difficulties-the Shura-Ozen River, cutting into the rock under the Coriolis force, erodes it, and the sand body has a tendency to shift in a south-east direction because of the winds.
These circumstances do not allow estimation of the accuracy of the model obtained from the UAV survey results by calculation of the models difference reliably. Therefore, to analyze the accuracy of elevation measurements only the areas with minimal influence of exogenous and anthropogenic processes were taken, respectively, which are the most stable-this is the floodplain part of the Shura-Ozen river and the flat part in front of the Narat-Tyube ridge. Within the stable areas, 490 points were randomly selected. The difference values for these points were calculated for two DEMs, one of which was a DEM based on the UAV survey, the second was based on the digitalized topographic map of scale 1:10,000. The following parameters were calculated as basic statistics: minimum and maximum value of difference, average, standard deviation, value of 25 (Q1), 50 (median) and 75% percentiles (Q3). The normality of error distribution was also analyzed by plotting a histogram of frequencies ( Figure 8). The results of statistical analysis of the values of the difference of two DEMs in the selected points are presented in Table 5.  Analyzing the obtained results, it can be concluded that the errors have a normal distribution, although the errors tend to move to negative values. In general, it can be stated that the model obtained using the UAV quite accurately describes the study area. Errors are mostly related to lower resolution of DEM, with which the modeling result is compared, relief variability, as well as small (on average 30-60 cm) difference in elevation systems between ellipsoid WGS84, in which UAV survey was performed and the Baltic system of heights, in which heights are given on topographic maps of the USSR, the sheet of which was used for construction of the DEM in 1963.
A textured model was created to visualize the results of terrain modeling. Since Agisoft Metashape is memory-intensive when creating mesh surfaces from a point cloud, the surface creation for texture mapping was done in the open source software Meshlab 2020.12. The textured model can be freely downloaded and viewed on the SketchFab page https://skfb.ly/6XMnQ (Figure 9).
The specific conditions of surveying of the selected object determined the difficulties encountered in the process of performing the work. These challenges can be thematically divided into several groups.

1.
The texture of the survey surface. Since June in Dagestan is characterized by hot and scorching sun, with usually zero cloudiness, aeolian forms, namely sands, create inconstant glare. In addition, the texture is extremely uniform in places. All these factors lead to the fact that the number of common points in the images tends to zero, and in places where it was still possible to detect a sufficient number of common points, there are specific noises and outliers in the dense point cloud. As a solution to this problem, it makes sense to increase the camera angle, which was described earlier.

2.
Difficult terrain. Since the altitude difference in the study area is 197 m, starting flight missions from one point becomes impossible, since the overlap of neighboring images decreases with increasing terrain height. In this case, there are no options other than to launch UAVs from the dune slopes and flat areas of the foothills.

3.
The hardware and software capabilities are undeveloped. In 2017, the DJI Phantom 4, in combination with Pix4D Capture, made it possible to survey small objects. The hardware and software capabilities differed depending on the operating system of the smartphone, but the basic principle remained the same. In 2020, much more advanced software features appeared, including taking into account terrain features when calculating flight altitude and overlapping images, as well as RTK solutions for the DJI Phantom 4. All of this would allow for a qualitatively different level of work; however, the cost of modifying the DJI Phantom 4 reduces its price appeal. Returning to the software, it is necessary to note that creation of the project for surveying a large territory "squares" is still not provided by software developers, so the only option for planning work of this extent is still either a manual survey or preparation of a series of kml or shp-files with borders of "squares" of interest in the GIS.
As noted earlier, some authors argue that surveying objects larger than 10 sq. km is unreasonable, which we categorically disagree with, since there is often no other way to obtain actual terrain information by other ways, while fixed-wing UAV surveys create inaccuracies and are more expensive. As the practice of work with low-cost multirotor UAVs has shown, such solutions make it possible to quickly obtain up-to-date and highly accurate information about surface elevations, as well as highly detailed orthophotomaps over a large area. For the Sarycum aeolian complex site, this is especially relevant, because, as also noted earlier, this natural monument is now actively vegetated, and the availability of ultra-detailed orthophotos allows making automated interpretation of dominant vegetation types and creating maps of vegetation overgrowth dynamics.
Drones 2021, 5, x FOR PEER REVIEW 15 of 18 and in places where it was still possible to detect a sufficient number of common points, there are specific noises and outliers in the dense point cloud. As a solution to this problem, it makes sense to increase the camera angle, which was described earlier. 2. Difficult terrain. Since the altitude difference in the study area is 197 m, starting flight missions from one point becomes impossible, since the overlap of neighboring images decreases with increasing terrain height. In this case, there are no options other than to launch UAVs from the dune slopes and flat areas of the foothills.

Conclusions
The obtained results of surveying the territory of the Sarycum area of the Dagestan Reserve showed that even low-cost unmanned aerial vehicles make it possible to create a digital elevation model of ultra-high resolution. For dynamic dune forms of relief, the availability of such information allows estimation of the intensity of aeolian processes at a qualitatively new level. At the same time, the works are largely automated and, as practice has shown, do not require large human resources. For the first time, a continuous survey was performed for this territory, as a result of which a high-resolution (0.5 m) DEM was produced for the area of 15 sq. km. Such work, taking into account the area of the survey, the characteristics of the object and the equipment used, has no analogues in the world.