Reconstruction of Earth Extreme Topography from UAV Structure from Motion Photogrammetry

UAV photogrammetry development during the last decade has allowed to catch information at a very high spatial and temporal resolution from terrains with very difficult or impossible human access. This paper deals with the application of these techniques to study and produce information of very extreme topography which is useful to plan works on this terrain or monitoring it over the time to study its evolution. The methodology stars with the execution of UAV flights on the cut slope studied, one with the cam horizontally oriented and other at 45o respect that orientation. Ground control points (GCPs) and check points (CPs) were measured for georeference and accuracy measurement purposes. Orthophoto was obtained projecting on a fitted plane to a studied surface. Moreover, since a digital surface model (DSM) is not able to represent faithfully that extreme morphology, information to project works or monitoring it has been derived from the point cloud generated during the photogrammetric process. An informatics program was developed to generate contour lines and cross sections derived from the point cloud, which was able to represent all terrain geometric characteristics, like several Z coordinates for a given planimetric (X, Y) point. Results yield a root mean square error (RMSE) in X, Y and Z directions of 0.053 m, 0.070 m and 0.061 m respectively. Furthermore, comparison between contour lines and cross sections generated from point cloud with the developed program on one hand and those generated from DSM on other hand, shown that the former are capable of representing terrain geometric characteristics that the latter cannot. The methodology proposed in this work has been shown as an adequate alternative to generate manageable information, as orthophoto, contour lines and cross sections, useful for the elaboration, for example, of projects for repairing or maintenance works of cut slopes with extreme topography.


Introduction
During the past decade there have been a quick technological development related to digital elevation modelling.For the most of geomorphic applications, topographic surveys have been mainly conducted using robotic total stations [1], or differential global navigation satellite system (GNSS) [2,3]).Nowadays, news technologies, such as terrestrial laser scanning (TLS) [4], aerial laser scanning (ALS) [5], softcopy photogrammetry [6,7], or aerial light detection and ranging (LiDAR) [8], have improved the accuracy of digital elevation models (DEMs) but they are often time-consuming and costly.
Landforms with complex topography can make these methodologies cannot be used and even prove dangerous for operators.Furthermore, in the most dynamic environments is necessary to have high temporal frequency at a very high spatial resolution images and DEMs for studying their evolution.This is especially necessary when the landform evolution may cause both human and economic disaster, as is the case of infrastructure built on these kinds of landforms.compared accuracies of photogrammetric projects differing in the image orientation: orthogonal images oriented to the terrain on the one hand, and classical vertical images orientation, on the other hand.They concluded that the first methodology is more adequate for this kind of morphology.
Therefore, in the light of the above, there has been a great development of UAV photogrammetry in recent years and is increasingly used in situations where classical photogrammetry is less efficient or simply not applicable.All this makes necessary to continue to deepen the development of specific methodologies to obtain accurate results using UAV photogrammetry in extreme situations, where classical photogrammetry has not even been able to arrive.
This paper provides a methodology to obtain photogrammetrically derived topographic information from UAV imagery for earth extreme and dynamic topography, where traditional techniques are dangerous or impossible to apply.The goals were to quantify the accuracy of the generated point cloud, and provide useful cartographic material for engineers, geologists and other technicians, of the terrain with extreme morphology.

Study area
The study area is a cut slope located at the N-340 road, in the province of Almería, southeast Spain, between Almería city and Aguadulce (Figure 1).This road was the longest in Spain (1248 km) and runs along the Mediterranean Sea coast from Cádiz to Barcelona.Nowadays there is a new parallel dual carriageway called A7 but N-340 has a heavy traffic, especially in sections that connect nearby cities.Like many coastal roads, N-340 has numerous cut slopes built in unstable rocky areas and there is a risk of landslides.The studied cut slope has 130 m long and the highest elevation difference between the road and the top of the cut slope is 70 m (Figure 2a).It has slopes almost vertical, even negative slopes, ie, part of the cut slope is suspended over the road.Furthermore, it suffered a landslide caused by torrential rains (Figure 2b) and it was necessary to have cartographic material to study its stability status and repair it.Traditional technics (Global Navigation Satellite Systems: GNSSs, terrestrial laser station, total station, …) were not possible to use here because the terrain morphology.So, the best option, and perhaps the only, was to carry out a UAV-photogrammetric project.

Materials and methods
The methodology applied was based on digital UAV-photogrammetry techniques, experimented by the authors in previous works, with adaptations to the particular characteristics of the studied morphology.

Image Collection
A rotatory-wing UAV with eight rotors (OktoKopter) and MikroKopter (Moormerland, Germany) electronic boards was used to take the images.The OktoKopter has a payload of 2.5 kg and a motion-compensated gimbal for the sensor was mounted on it.The sensor was a Nikon 3100 reflex camera with a lens with 16 mm fixed focal length.This camera has a complementary metal-oxide semiconductor (CMOS) sensor of 14.8 effective megapixels with a size of 355.7 mm2 (23.1 × 15.4 mm). Figure 3 shows the entire system.The flight plan was programmed and loaded to the UAV using the MikroKopter-Tool [38] software.In order to photography the whole complicated cut slope surface, two flights were executed varying the photography axis: one with horizontal orientation and one with the axis tilted 45º downwards.During all these flights, the UAV remained in a vertical plane at a distance of approximately 50 m to the studied surface.As the image resolution was adjusted to 4240 × 2832 pixels, the ground sample distance (GSD) was 1.86 cm.The flight altitude varied from the road level to 100 m, the overlap between pictures was 90%.
Furthermore, 3D coordinates of several points scattered on the studied surface were measured with a total station without reflector.Altitude of these points could not be greater than around 35 m from the road level because the angle of the station telescope was very high and was not possible to look across the telescope.Of these points, those not located on the cut slope were measured using GNSSs.Two of these three points were used for orientating the measured points with the total station to a local coordinate system.The GNSSs measurements were made working with differential corrections in real-time kinematic (RTK) mode, with the base station on a geodesic pillar located closer than 1 km from the studied site.Both rover and base GNSS receivers were Trimble R6.For the RTK measurements, these dual-frequency geodetic instruments, which track the Global Position a b System (GPS) and Global Navigation Satellite System (GLONASS) signals simultaneously, have manufacturer's stated accuracy specifications of ±8 mm + 1 ppm RMS horizontal and ±15 mm + 1 ppm RMS vertical.Therefore, the maximum horizontal and vertical RMSs were ±9 and ±16 mm, respectively.The total station without reflector was a Stonex, STS22R model.This instrument has a measurement range, without reflector or white target, of approximately 100 m, and an accuracy of 5mm + 2 ppm in reflector-less fine mode.

Image Processing
The images were processed using the software Pix4Dmapper Pro, version 3. 1 [39].This package incorporates the SfM procedure described in [40] and has a three-step workflow.The first step of the images processing is to compute keypoints on them, which will be used to find matches between the images, and from these, to run an automatic aerial triangulation and bundle block adjustment.Images taken in this work were not geotagged.The internal calibration parameters, the relative camera position and orientation corresponding to each picture and the 3D relative coordinates of a sparse point cloud of the terrain are the results of this first step.To refer all these results to a local coordinate system the coordinates of measured GCPs and CPs are added to the project and the point are marked on the pictures where they appear.In this work, the coordinates of 18 points located on the cut slope were measured with the total station.It was very difficult to locate points on the cut slope which could be identified in the photos to georreference the point cloud.Vertex of characteristics shapes or metallic pieces of older works on the cut slope were used.Furthermore coordinates of 8 points scattered on the road and marked with targets of A4 format size (210x297 mm) were measured with total station and with GNSSs to refer the local coordinates system of the total station to the local system.Figure 4 shows the target used on the road (right) and one of the metallic piece used as target in the cut slope (left).The second step achieves a densification of the point cloud and a more detailed 3D model than that of step 1 is obtained.Finally, in the third step, a grid DSM can be generated at a specific grid size and the orthophoto is exported at a pre-selected resolution.The bundle adjustment can be carried out using only three GCPs, but it is advisable to use more than three to obtain optimal accuracy [35,41,42,43].
As the study surface is almost vertical, the orthophoto obtained when projecting on a horizontal plane does not offer valuable information, even confusing information, because there are zones where for a given X and Y coordinates, two or more Z coordinates can be present.To avoid this circumstance, a plane was adjusted to the cut slope surface and it was used to project for building the orthophoto.This plane was determined by fitting to the terrain point cloud obtained in the second step of the image processing, previously described.For this task, only the points centered on the interest area were taking into account.The resulting orthophoto will provide valuable information for the study of the cut slope surface: location and description of cracks, location of rocks with danger of falling into the road, take measurements or get information to plan works on the cut slope.

Point cloud management
Coordinates of dense point cloud obtained from photogrammetric process were referred to as UTM Zone 30N (European Terrestrial Reference System 1989, ETRS89) and the elevation as the mean sea level using the EGM08 geoid model.A DSM obtained from this point cloud would lose a lot of information due to the surface geometric characteristic mentioned above: regardless the reference plane taken into account, for a given X and Y coordinates, there will be areas where two or more Z coordinates are found, and it is not compatible with a DSM.Then, to get the most out of the information provided by the point cloud, an informatics program was developed to obtain, from the raw point cloud, contour lines on one hand and cross sections produced by vertical and perpendicular planes to the fitted plane on the other hand, at any point of the work area.This kind of terrain representation shows a very realistic vision of the surface shape and it is useful to plan and execute works on the terrain under study.The program had two main parts.In the first part, contour lines and vertical section are extracted from the point cloud.For this task, the program determines maximum and minimum X, Y and Z coordinates presented in the point cloud and the user can select the extreme coordinate to fix the box which will be studied.A dialog box (Figure 5a) shows the extreme coordinates of raw point cloud (blue color in the figure 5b) and it allows to fix the limits of the box which will be studied (green color in the figure 5b).Names of the file containing raw cloud data and processed data should be indicated in the bottom of the dialog box.The output file is structured by three columns, containing coordinates X, Y, and Z, respectively, of those points inside the box to be studied (green box in 5b).When this file is created, the second step starts showing a dialog box (Figure 5c) for generating horizontal and vertical cut slope sections (blue and red sections, respectively, in figure 5d).In this figure, ω is the fitted plane, ϕi and ϕi+1 the relative position of two consecutive horizontal planes which generate two consecutive contour lines, πi and πi+1 are two consecutive vertical planes which generate two consecutive cross sections.A and B are the points defining the intersection straight line between the horizontal plane with the minimum coordinate Z considered, and ω.So, Z range will be Zmax-Zmin and XY range will come from the distance between A and B (dAB).Number of horizontal sections will be (Zmax-Zmin)/Δv and number of vertical sections will be dAB/Δh.In this dialog box, "Width section" mains the maximum distance from one point of the cloud to the vertical cut plane (π1, π2, …) or horizontal cut plane (ϕ1, ϕ2, …), to be considered included in the section.This was taken into account to ensure a minimum number for a good definition of the section.In the top of the "Section generation" dialog box, the coordinates defining the cloud box to be studied, Z variation range and distance between A and B points (figure 5c) are showed.Then, distances in meters between consecutive vertical (Δv) and horizontal (Δh) section are fixed (in the Figure 5c, both of them are fixed to 1 m).Furthermore, the width vertical and horizontal sections are fixed to 0.01 m but this value can be changed.This is a critical adjustment because the accuracy of sections representation depends on this value.If it is too low, few points will be extracted and the section will be poorly defined.On the other hand, if this value is too high, a lot of points will be extracted and the section will be confusedly defined.So, to get an optimal value, an informatics program was developed to compare results taking into account several values of cutting plane width.
For a given Z, contour lines file are structured by two columns containing coordinates X, Y coordinates.Files containing cross sections are structured by two columns too, one for de distance measured on the intersection of πi with the horizontal plane, and the other column for the coordinate Z.One file is created for each section.The second part of the developed program allows display the generated sections.Its interface shows three graphical windows (figure 6): the main one for the orthophoto, projected on the fitted plane, other for the contour line, above the orthophoto windows, and the third for the cross section.When a click is done on the cut slope image, a cross appear on that, representing the intersection of the horizontal (red) and vertical (blue) planes with the terrain, and immediately the contour line and the cross section are drawn in its respective window.Extreme coordinates of each window, contour line elevation are showed.Furthermore, when the cursor is on the orthophoto, coordinates are showed in the bottom.

3.4.Accuracy assessment
For every CP, the accuracy assessment in Easting (X), Northing (Y) and height (Z), was performed by comparing the CPs measured coordinates with the interpolated coordinates from the four nearest points of the dense cloud generated in the photogrammetric process, resulting in RMSEX, RMSEY and RMSEZ accuracy measures, respectively: where: n: number of CPs.Xsi, Ysi and Zsi: X, Y and Z coordinates measured with the total station for the ith CP.Xci, Yci and Zci: X, Y and Z coordinates of the closest cloud's point to (Xsi, Ysi, Zsi).
Interpolation method was inverse distance to square power.An informatic program was developed in Visual Basic V6.0 language for carrying out this task.

Results and discussion
Extreme coordinates of the prism containing the raw point cloud were (540117, 4074712, 0) and (540453, 4074967, 143), which means dimensions of 336 × 255 × 143 m and 11048542 points were generated in the photogrammetric process.As the interest area was smaller than all covered area, this was reduced using the informatics program (Figure 4a and 4b).The extreme coordinates of the study area were (Xmin = 540220, Ymin = 4074750, Zmin = 20) and (Xmax = 540350, Ymax = 4074800, Zmax = 90), which means dimensions of 130 × 50 × 90 m and the number of point was 2933590.In this way, all tasks to be carried out on the points cloud will be much faster by having to handle fewest points.which was used as projection plane to make the orthophoto.Planes for cross sections were perpendicular to the AB straight line (Figure 7), which is the intersection of ω with the horizontal plane Z=20 m: = 4153968.907+ 87.714× (5) Table 1 shows the coordinates of the 18 points measured on the cut slope and the 8 measured on the road (Xm, Ym, Zm).Points 1, 9, 14, 19 and 22 were used as GCPs and the rest were used as CPs.Furthermore, points 23 and 26 were used for orienting the coordinates measured with total station.In the table 1 and Figure 7 can be observed than the maximum altitude of the measured points was 51.677 m.It was not possible to measure points located at higher altitudes doubt to the limitations exposed in the material and methods section.Futhermore, table 1 shows estimated coordinates from points cloud as described in 3.4 section (Xe, Ye, Ze), error for each coordinate (EX, EY, ZZ) and RMSE in X, Y and Z directions.Error ranges were from -0.085 to 0.085 m (0.170 m) for X component, from -0.131 to 0.099 m (0.230 m) for Y component and from -0.100 to 0.118 m (0.218 m) for Z component, indicating similar ranges for each component.RMSE for X, Y and Z directions were 0.053 m, 0.070 m and 0.061 m respectively which are similar or better to those reported in other works carried out in similar conditions.So, Lucieer et al. [32] used Agisoft PhotoScan Profesional 0.85 software for the 3D terrain reconstruction and reached a geometric accuracy of 0.060 m in planimetry and 0.044 m in the Z component.They worked on an Antartic surface, the flight altitude was 50 m above ground level and the camera was similar to that used in this work.Turner et al. [44] worked in similar conditions than Lucieer et al. [32] and reported accuracies of 0.100 for planimetric component and 0.150 for Z component.Agüera-Vega et al. [43] used similar UAV-system than that used in this work to study the accuracy of digital surface models and orthophotos derived from unmanned aerial vehicle photogrammetry on terrains with different morphologies and comparing results from different flight altitudes and number of GCPs.At 50 m flight altitude and 5 GCPs, which are similar conditions to those of our work, they report RMSE values around 0.050 m for X, Y and Z components.Hugenholtz et al. [45] reported vertical RMSE values of 0.106 m and 0.097 m derived from a photogrammetric project on a stockpile, with a rotatory-wing UAV and 100 m flight altitude, before and after a portion of this pile was removed.In complex-morphology terrains, Carvajal-Ramírez et al. [37] yield planimetric and vertical accuracies of 0.058 m and 0.100 m respectively working in a road cut slope.They used a rotatory wing UAV and carried out two flights for the photogrammetry project, one with the cam in vertical position and the second with the camera oriented 45º with the terrain.Table 1.Coordinates of measured points (Xm, Ym, Zm), estimated coordinates (Xe, Ye, Ze), error (EX, EY, ZZ) and RMSE in X, Y and Z directions.All coordinates are referred to as UTM Zone 30N (ETRS89) and the elevation as the mean sea level using the EGM08 geoid model. 1  As regard the orthophoto, figure 8 shows the orthophoto obtained using a horizontal plane (a) and that obtained using the fitted plane (b).As the surface under study was the cut slope, orthophoto projected on a traditional horizontal plane shows less information than that projected on the plane used in this work because this offers the optimal vision to detect, describe and to measure cracks, location of rocks with danger of falling into the road, take measurements or get information to plan works on the cut slope.Furthermore, in the orthophoto of figure 8a terrain located at the top of the cut slope hidden to another located at the bottom.This can be observed around the center of the picture.Nevertheless, orthophoto of figure 8b shows with much more detail than figure 8a the area under study.In addition, contour lines and cross sections yielded from the dense point cloud give valuable information to the engineers, geologists and other technicians for proposing the necessary works on the cut slope to repair it and avoid a new landslide.For this terrain, if a horizontal plane of projection is used, the DSM would not have represented the reality since for a given coordinates X, Y, a single value of Z coordinate would be given by the DSM, but this is not the case for certain zones of the studied terrain.The developed informatics program has made possible to obtain contour lines directly from point cloud, so these represent the terrain truthfully.
As stated in section 3.3, a critical adjustment of this program was to fix the cutting planes width for generating contour lines and cross sections.Extraction of intersection points between point cloud and cutting plane was carried out taking into account 0.5 cm, 1 cm, 2 cm and 5 cm plane widths and studying by comparison results using the developed informatics program.Figure 9 shows the results for three cross sections (CS1, CS2 and CS3).The number in the bottom of each section is the number of point extracted from the point cloud.When 0.5 cm width is fixed, observed accuracy of a b cross section representation is not enough because exist quite gaps between extracted points.When 5 cm width is fixed, appear representation sections with an elevate number of extracted points which make confuse the section.Nevertheless, differences between sections representation taking into account 1 cm and 2 cm width were not remarkable, both values yield well defined representation sections.Finally, for our study, 1 cm width was fixed to generate the contour lines and cross sections.Figure 10 shows the contour levels yielded from DSM (a), and from point cloud (b).In the first case, no contour line crosses another one, but in the second case it can be observed than several contour lines cross other lines.Similary, the developed program generates the cross sections from point cloud yielded in the photogrammetry project.Figure 11 shows several cross sections in which it can be observed that for a given planimetric point X, Y, there may be more than one value of Z because the terrain morphology which can be represented if these sections are generated from point cloud (blue lines in figure 11), but no if these are generated from a DSM (red lines in figure 11).

Conclusions
The use of UAVs as platforms that can be coupled with different types of sensors to capture information at very high geometric and temporal resolution is being implemented from several years ago.This paper proposes a methodology for reconstruct extreme topography that has been developed in order to obtain useful information for engineers, geologists and other technicians which are interested in the study of this type of these surfaces.The equipment used in this study, composed of a light UAV and a non-metric camera has been shown as the only system capable of achieving that goal since for extreme topography, terrestrial instruments as global navigation satellite system (GNSS), terrestrial laser scanner and aerial or space surveys are no suitable.In this case, terrestrial instruments are useful for the measurement of GCPs and CPs.
The use of SfM tecniques, incorporated in moderns UAV-photogrammetry commercial software, allows us to obtain a dense point cloud representing the studied surface from which, Digital Surfaces Models and orthophotos can be derived.Point density has been enough to represent accurately the terrain.So, RMSEs found in X, Y and Z direction have been 0.053 m, 0.070 m and 0.061 m, respectively, similar to the best accuracies reported in other works for similar conditions.Traditionally in topography, the reference plane to project the orthophoto has been horizontal but there are times when it is necessary to project on a non-horizontal plane, such as surfaces with a very high slope.In this sense, the methodology proposed in this work, consisting of projecting on the plane that best fits the surface under study, has been shown as an adequate alternative since the derived orthophoto allows to carry out different types of studies.
In addition, DSM may not reliably represent surface geometry and is necessary to work with cloud point to generate useful information.For this purpose, the developed software in this work has proved suitable to generate manageable information, as contour levels and cross sections, useful for the elaboration, for example, of projects for repairing or maintenance works of cut slopes with extreme topography.Figure 12 shows the cut slope aspect after the initial works projected using the information derived from the methodology described in this study.Point density derived from UAV-photogrammetry has allowed us to represent the terrain from contour lines and cross sections, which is the used format to elaborate the plans to carry out works on terrains.For each cutting plane that generates a contour line or cross section, it has been sufficient to consider one centimeter width to obtain accuracy representation of these lines.

Figure 1 .
Figure 1.Geographical location of the study area.The red arrow indicates the location of the study area.Red line represents the N-340 road and the blue line represents the A7 dual carriageway.

Figure 3 .
Figure 3. UAV equipment used in this work.

Figure 4 .
Figure 4. Detail of the target used for the measurements of GCPs and CPs located on the cut slope (left) and that used for those located on the road (right).

PreprintsFigure 5 .
Figure 5. Program developed to obtain horizontal and vertical sections files: (a) dialog box for delimiting the work area; (b) graphical representation: blue prism indicates the whole space covered by the pictures and green prism indicates the selected area; (c) dialog box for generating level and cross; (d) graphical representation of the sections generation.

Figure 6 .
Figure 6.Interface of the informatics program developed to draw contour levels and vertical sections from the point cloud.Red straight line represents the horizontal plane and the blue red line the vertical plane cutting the point cloud.Contour line is represented in the top window by red dots and cross section is represented in the right window by blue dots.

Figure 7 .
Figure 7. Schematic (left) and realistic (right) representation of the spatial distribution of the measured points as GCPs (blue pots), CPs (red dots), and points for orientating total station measurements (green dots).

Figure 8 .
Figure 8. Orthophoto obtained using a horizontal plane (a) and that obtained using the fitted plane (b).

Figure 9 .
Figure 9. Example of three cross sections (CS1, CS2 and CS3) obtained with the informatics program developed, taking into account a cutting plane width of 0.5 cm, 1 cm, 2 cm and 5cm.

Figure 10 .
Figure 10.Contour levels yielded from DSM (a), and from point cloud (b).

PreprintsFigure 11 .
Figure 11.Example of four cross sections generated from DSM (red lines) and from point cloud (blue lines).

Figure 12 .
Figure12.Cut slope aspect after the initial works projected using the information derived from the methodology described in this study.
points used as GCPs.Rest of points were used as CP.2points used for orienting total station measurements.