Mapping quality prediction for RTK/PPK-equipped micro-drones operating in complex natural environment

Drone mapping with GNSS-assisted photogrammetry is a highly efficient method for surveying smallor medium-sized areas. However, the mapping quality is not intuitively predictable, particularly in complex environments (with steep and cluttered terrain), in which the quality of the real-time kinematic (RTK) or postprocessed kinematic (PPK) positioning varies. We present a method to predict the mapping quality from the information that is available prior to the flight, such as the flight plan, expected flight time, approximate digital terrain model, prevailing surface texture, and embedded sensor characteristics. After detailing the important considerations, we also present the concept of global precision within the context of minimal and efficient ground control point placement in a complex terrain. Finally, we validate the proposed methodology by means of rigorous statistical testing against numerous experiments conducted under different mapping conditions.


Motivation
Unmanned aerial vehicles (UAVs) are becoming an important tool for surveyors, engineers, and scientists as the number of available costeffective and easy-to-use systems is rapidly increasing. These platforms offer an attractive alternative to mapping small areas with centimeterlevel resolution. Numerous successful applications have been reported: repetitive surveys of buildings, civil engineering structures or construction sites, as well as land monitoring and precision farming Colomina and Molina (2014).
The aerial images that are acquired in this manner need to be georeferenced accurately and precisely, either to localize specific features of interest in the images, such as cracks or corrosion spots in concrete structures, or for three-dimensional (3D) modeling. The typical airborne photogrammetric acquisition of a drone is depicted in Fig. 1, in which the flight plan is conceived to follow the terrain contours. Aerial triangulation, in which photographs are linked through image observations of so-called tie-points, is a common method. The 3D model data in object space are obtained using multiple ground control points (GCPs) Triggs et al. (1999). These points limit model distortions owing to error accumulation and are also used for quality control.
The process of initializing a dense and uniform ground control network is notoriously time and cost intensive (as physical signalization and separate surveying are required). Moreover, the process is sometimes complicated by the terrain, severely limiting the response time and the application of UAVs in difficult environments and/or under strict accuracy requirements . This is why initially research Rehak et al. (2014) and later industry senseFly (2015) and Mavinci (2016) focused on using airborne photogrammetry methods for UAVs, with the aim of removing or mitigating the need for GCPs. In assisted aerial triangulation, ground control is replaced with aerial control: a survey-grade GNSS receiver can provide centimeter-level accuracy positions for each camera station, provided that the reference and rover receivers share sufficient satellites (Fig. 1signal obstruction) in an appropriate geometric configuration (as in the right part of Fig. 2) and that the carrier-phase signal is of sufficient quality. This remains challenging onboard small UAVs Stubbs and Akos (2016). This technique is very popular at present in real-time kinematic (RTK) and post-processed kinematic (PPK)-enabled UAVs. Moreover, survey-grade GNSS receivers are now available on the market senseFly SA and Mavinci (2016). If an inertial measurement unit of sufficient quality is also available; for example, as in Mian et al. (2015), orientation control can be employed, which may be required in difficult mapping scenarios such as corridor mapping Rehak and Skaloud (2015). In recent years, variations of such techniques have been developed to deal with specific peculiarities of UAV-based aerial photogrammetry, including that of relative position/orientation control Blázquez and Colomina (2012) and Rehak and Skaloud (2016) and raw observation, and multi-sensor adjustment Cucci et al. (2017).

Challenges
Despite the extensive range of recent possibilities offered by both hardware and post-processing methodologies, the survey design remains crucial for the efficient deployment of UAVs: the operator needs to select factors such as suitable sensors, the configuration and number of GCPs, and the flight plan and its time schedule so that the accuracy requirements are fulfilled, while minimizing the survey time and cost as well as the chances of repeating all or part of the survey owing to unsatisfactory results.
In difficult mapping scenarios, such as corridor mapping or cluttered environments, in which high variations in the topography cause the image overlap and sampling distance to vary substantially, it becomes difficult, even for expert operators, to predict the quality of the final mapping products. For example, Fig. 1 shows how the 3D object point precision varies depending on the individual performance of directly observed trajectory information, individual image overlap and terrain texture. A direct consequence thereof is that either too many or too few GCPs are installed, leading in the first case to an unnecessary cost increase, or in the second case, requiring a repeat of all or part of the mission.
Another source of uncertainty in the survey outcomes arises from the quality of the aerial position control, which in the first instance affects drone guidance and in the second instance influences the image orientation quality. The latter aspect is partly driven by the signal-tonoise ratio (SNR) of the carrier-phase GNSS signal on two or more frequencies, which is determined by the drone hardware components and their assembly. Furthermore, the satellite constellation available at the time of the flight affects both the capability of obtaining a fixed GNSS position throughout the entire trajectory (operational safety), and the ability to resolve and maintain the double-difference carrier-phase ambiguities that are essential for observing the camera-coordinates with centimeter-level precision correctly.
The importance of GNSS constellation planning is well recognized among surveyors. However, it is typically performed for representative positions with arbitrary elevation masks, and the terrain is not necessarily taken into account ( Figs. 1 and 2). Furthermore, the fact that the satellite occlusions vary along the drone trajectory and as a function of time is neglected. Fig. 8 presents a situation in a cluttered environment in which the GNSS signal quality varies with time within the same flight.

Conventional approach
The commonly used method for UAV photogrammetry is to fly according to a flight plan and to post-process the recorded data to produce surface-related products such as 3D models and ortho-images. The images are first oriented either approximately by using navigation sensors or accurately using the computationally intensive process of tiepoint identification and matching. All inputs, namely the image coordinates, object coordinates, and possibly the camera position/attitude (GNSS + INS) are combined with appropriate weights within a bundle adjustment (BA), from which the mapping quality factors are derived. If the quality is unsatisfactory, the process needs to be restarted to acquire additional images, improve the aerial control or add more GCPs (a new target if new photographs are acquired or the measurement of existing identifiable targets).
The simulation of the geometry for the 3D reconstruction problem has already been explored, although it is not necessarily part of drone mission planning. For example, in Piatti and Lerma (2013), a method was described to simulate the geometry of a photogrammetric survey. In Buffa et al. (2016), noise was added to the simulated image measurements, and 3D reconstruction was executed and compared with the ground truth to quantify the expected survey precision. In this work, the reconstruction was executed only once, while in Luhmann (2011) and Tushev and Sukhovilov (2017), a Monte Carlo approach was used, exploiting several runs of the reconstruction step to enable estimation of the covariance matrix of the expected precision. In O. Mason (1994) and Olgague (1998), the expected precision covariance matrices were obtained via linear covariance propagation. However, the inclusion of aerial control adds another level of complexity to this process, especially if the spatio-temporal variations in the expected aerial control quality in obstructed environments are considered.

Proposition
In this study, we consider the fact that satellite visibility changes while the drone is moving during the mission planning phase (that is, both the satellites and the portion of visible sky are moving), and predict the capability of the drone guidance (sustaining the minimum number of visible satellites) as well as the success of the kinematic ambiguity resolution in real-time (RTK) or post-processing (PPK). Thereafter, we combine this information with either planned or realized (immediately after landing) image overlap, the prevailing texture of the scene, and the distribution of the GCPs to predict an accuracy map over the surveyed area that may be inspected by the operator via interactive visualization. The mission parameters can subsequently be altered to improve the estimated accuracy. If the predicted accuracy does not meet the operator requirements, stability analysis of the entire BA network is performed to determine the lowest-accuracy locations and to propose the optimal placement of additional GCP(s). Similar analyses are performed following the flight to display the predicted accuracy map with the actual camera perspective center position/attitude and its  quality prior to processing the actual images. The remainder of this paper is organized as follows: In Section 2, we present the methodology for simulating the required inputs. We detail the relationship between the geometry of visible satellites and the probability of successful ambiguity resolution, which is used to derive the aerial position uncertainties for each camera station and time. These are encapsulated into one parameter that serves as guidance for proposing the optimal mission time of the day. Moreover, we describe the additional inputs related to the probabilistic tie-point distribution based on the scene texture. Thereafter, we present the local and global quality factors that are useful for precision evaluation. Subsequently, in Section 3, we present an analysis of the network for its principal weaknesses and demonstrate the most efficient improvement in the mapping precision by placing additional GCP(s) within the indicated area(s). Finally, in Section 4, we present experiments conducted in mountainous environments to validate the proposed methodology.

Concept
As indicated in the left part of Fig. 3, the conventional loop of data acquisition, post-flight processing, and derivation of the quality control at the end of the survey process is not ideal, particularly in complex scenarios in which the aerial position control quality is likely to vary with time. Therefore, we aim to extend such tools with the spatiotemporal elements of the aerial control in an obstructed environment in combination with probabilistic measures for tie-point placement to simulate the acquisition process realistically and derive the expected quality prior to the flight. Overall, we aim to shorten the interaction loop before obtaining acceptable results regarding the conventional approach, in which the quality estimators are based on real data, possibly leading to repeated missions (left vs. right organogram of Fig. 3).
As schematically depicted in Fig. 4, the proposed method is based on the following inputs: i) the flight plan covering the area with the aimed overlap and resolution represented by the mean ground sampling distance; ii) the foreseen flight date and time; iii) a coarse digital terrain model (DTM) 1 ; and iv) the nominal camera interior orientation and whether this needs to be re-calibrated Fraser (1997). Based on these elements, simulated image and aerial control observations are generated, along with their expected precision (which varies temporally and spatially for aerial position control), and a full adjustment is run. This forms a basis for formulating realistic predictions of the mapping quality that can be expected from the actual survey at a specific date and time of the day, with details over the entire area. On this basis, the user can interactively modify the flight plan (for example, by selecting a different time, payload, flight line, and overlap) and ground control network (for example, by adding GCPs or changing their placements) until the quality requirements are fulfilled, thereby maximizing the probability of a successful data collection campaign and minimizing the costs thereof.
If the predicted quality is higher than the requirements, the operating procedure could be simplified, thereby enabling a reduction in the cost, prior to executing the flight and related survey. If the quality is unsatisfactory and improvement is sought through additional GCPs, further analysis of the network is performed to identify its weakest characteristics and areas in which the placement of GCPs would be the most effective.

Basic relation
The basic relation for 3D object restitution in photogrammetry is the collinearity in Eq. (1), which is in its most simple form based on the pinhole camera model. It relates the 3D coordinates of a tie-point P t m in a mapping frame (m) to its image coordinates x y , c c (measured from the principal point on a theoretical focal plane, the distance of which to the camera perspective center P c is 1) via a scale factor , rotation matrix R c m , and the 3D position of the camera perspective center in the mapping frame P c m : Additional parameters modeling departures from collinearity need to be assumed for real cameras, as in Section 2.4. The observations in relation to different elements of Eq. (1) are described in the following Sections 2.3 and 2.8.

Flight plan
The flight preparation is achieved by a mission planner such as Oborne (2019) and Gandor et al. (2015). The principle of such software is to assist the user in designing the trajectory that the drone is supposed to fly according to certain criteria. It requires an approximate DTM such as NASA and METI (2011) to achieve the desired parameters; for example, the scale (ground sampling distance) and lateral overlap, while ensuring that the drone remains below the maximum height fixed by the regulations.
Apart from providing the input for the drone guidance, the mission planner supplies the rough values of the camera poses for the simulation (Fig. 4). The positioning precision of the RTK/PPK is assumed to exhibit time-varying characteristics and is further detailed in Section 2.7.
(1) to the observed pixels ( , ) x y in a real image. It includes the sensor position in the camera frame and additional parameters related to departure from collinearity. A commonly used physical model is that of Brown Förstner and Wrobel (2015), which is defined in Eq. (2).
Where p p ( , ) x y is the position of the principal point, p d is the principal distance, K ( ) 1,2,3 and P ( ) 1,2 are the radial and tangential distortions parameters and = + r x y c c 2 2 2 . The camera parameter values can be obtained from a previous project; alternatively, the approximate magnitude of p d from the sensor datasheet can be used. The parameter uncertainty and need for re-calibration will influence the mapping precision. A reasonable selection considers the extreme cases (in which all parameters are either free or fixed: Fig. 4, steps 3/5), as well as an intermediate scenario (Fig. 4, step 4). The latter reflects the case in which the lens-related parameters may be temporarily stable (particularly on a rigidly mounted prime lens used at a similar temperature), while the values of p p , x y , and may vary slightly among flights on different days and need to be readjusted.

Ground control
In the absence of aerial control, at least three GCPs are required for absolute orientation in the mapping frame, whereas more may be necessary to improve the mapping quality. Their position and placement are typically part of the project design, for the purpose of which their approximate locations can be obtained by clicking on the digital map employed in the flight planning (Fig. 4, step 6). Alternatively, existing coordinates can be loaded in a file together with their accuracy values, which is dependent on surveying technology (such as GNSS and leveling) to be considered as the observations: The answer to the influence of the GCPs on the mapping accuracy in terms of the number, quality, and distribution is part of the analysis presented in Section 3.2.

Tie-points
We consider the tie-point as a general term for an image feature, the center of which is observed in at least two images by tie-point detection and a matching algorithm. Although the quality and quantity of the automated tie-point detection and matching vary slightly with the employed algorithms Miksik and Mikolajczyk (2012) and Rublee et al. (2011), the surface texture is the most important factor. Indeed, the tiepoint density on surfaces with a homogeneous texture (for example, water, fresh snow, and certain types of vegetation) may be low or even zero, whereas it may be extremely high in texture-rich environments (such as built-up areas). As the real texture is unknown at the prediction stage, the outcome of such a process can only be probabilistic. To determine such a probabilistic measure, we evaluated several flights at a 100 200 m height above ground level, in which we considered several texture categories and calculated their mean tie-point densities, as displayed in Table 1. Thereafter, we used the image residuals in these categories to estimate the mean variance of the image observations empirically.
It is reasonable to assume that the selection of the prevailing texture is guided by a priori knowledge regarding the terrain. Thereafter, the corresponding empirical values of the image tie-point density in Table 1 are converted into the tie-point density on the ground via flight parameters to generate tie-points on the approximate digital surface (Fig. 4, step 8) using a stochastic process (Fig. 4, step 9). These tie-points are subsequently merged with the other points (GCPs; check-points) that are manually entered by the user (Fig. 4, step 6) to generate image observations. The influence of the tie-point distribution and quality is discussed in Section 4.3.

Aerial position control
The relation between the GNSS-derived position and camera position is given by the following equation, where A is the 3D position given by the GNSS antenna, P c m is the camera position in the mapping frame, R c m is the camera orientation, and a c is the lever arm between the center of the GNSS antenna phase and the camera perspective center. In manned airborne missions, shift parameters (or drift parameters) could be added to the above equation to absorb effects due to incorrect determination of ambiguities within a block or per flight-line. In drone mapping the shorter distance to the base improves the reliability of ambiguity determination due to differential atmosphere, and thus improve greatly the accuracy. However, wrongly determined ambiguities could still happen if the geometric configuration is weak, and that due to proximity of obstacles could happen anytime within a flight-line. For these reasons, the modelling by shift (drift) parameters is no longer appropriate and it is better to plan the flight execution so that acceptable observation conditions are maintained throughout the whole trajectory. As mentioned previously, apart from the signal strength (SNR) in the code and carrier-phase observations, the precision of the GNSS positioning is mainly dependent on (i) the number of observed satellites and (ii) the satellite-to-receiver geometry. These geometrical factors can be analyzed in advance from the mission plan and coarse digital elevation model (DEM) over a time-span that is specified by the user (Fig. 4, step 10). For this purpose, the most recent information regarding the GNSS approximate satellite position is obtained from the almanac GSA (2020) or emphemerides. The time of each photo is predicted by the mission plan with respect to the specified mission start, considering its planned position and the nominal UAV speed. The positions of the GNSS satellites are computed from the almanac at each planned camera location and time (Fig. 4, step 11), while a ray-tracing algorithm allows for determining whether this satellite is in the line of sight; that is, not obstructed by the terrain Nowak (2017).
Once the visible satellite configuration is determined, the positioning covariance matrix is predicted by the calculation of the co-factor matrix scaled by the ranging precision (Fig. 4, step 12). The latter depends on the foreseen positioning mode used in Eq. (4): with carrier-phase differential corrections or without; that is, standalone positioning. The former can be achieved either during the flight (RTK) or thereafter (PPK). Although the noise level of the PPK is generally lower than that of the RTK, as discussed later (Section 4.5), the crucial aspect affecting the ranging is the capacity of resolving the carrier-phase ambiguities. This is related to the size of the ambiguity search space and its geometrical shape, which is proportional to the principal elements of the GNSS covariance matrix or its compound metric, known as the position dilution of precision (PDOP) Teunissen et al. (1999). We suggest employing the empirical probabilistic function proposed by Schaer et al. (2009), based on experience with many helicopter flights close to the terrain in a mountainous environment, to inflate the previously derived covariance matrix by means of such a function. This is illustrated in Fig. 5 as a function of the PDOP and number of observed satellites.
The color scheme of the probable positioning quality in each photograph is displayed in the flight plan, as illustrated in Fig. 8, to guide the feasibility of the foreseen flying time intuitively. If the result is not satisfactory (more red than green over a specific area), the flight time can be modified (by means of back-looping from Fig. 4, step 12 to Fig. 4, step 10) to analyze another interval.

Precision map
The precision map is based on the predicted covariances of the individual tie-points in the object space. These are obtained by co-variance propagation of the previously described observation equations with the simulated measurements. This problem is part of the simultaneous optimization of the camera orientation and 3D location of the ground points from the image observations and aerial ground control; that is, the BA. In BA, the residuals associated with each observation are minimized as a function of the parameters by means of weighted nonlinear least-squares methods; for example, as described in Triggs et al. (1999) and Förstner and Wrobel (2015) or for specific aerial control of UAVs in Rehak and Skaloud (2016).
Firstly, a ray-tracing algorithm is applied to construct an adjacency list that links the camera orientations and tie-points in view in its image.

Table 1
Texture categories with mean tie-point density and precision. Fig. 5. Empirical GNSS precision positioning as function of number of visible satellites and co-factor elements. E. Cledat, et al. ISPRS Journal of Photogrammetry and Remote Sensing 167 (2020) 24-38 This set of observations, together with the terrestrial position (Section 2.5) and aerial position (Section 2.7) permits the photogrammetric network to be simulated (Fig. 4, step 13), with the weight matrix defined as the inverse of the covariance matrix of the observation vector . Assuming that the observation errors are zero mean, Gaussian distributed, and independent, a single standard deviation for each observation or class of observations can be specified. The parameters are the camera poses (position and attitude), tie-point positions, and intrinsic camera parameters (for example, principal distance, principal point, and/or lens distortion parameters).
These parameters are concatenated in the state vector x . Let f be the function describing the observation model; that is, the function such that = f x ( ). All of the collinearity equations of all tie-point measurements (Section 2.4), aerial position observations (Section 2.7), and terrestrial observations Section 2.5) are enclosed in this f function. Algorithms that are used to solve least-squares problems (such as Gauss-Newton, Levenberg-Marquardt or the dog-leg algorithm) require the Jacobian matrix of f with respect to x. This matrix, which is known as the design matrix, is denoted by A in the following. The Lie group theory Strasdat (2012) is used to handle the derivations of the rotation matrices describing the camera orientation.
In the scope of covariance propagation, the values of the parameters that are known a priori are used to construct A. Then, the covariance matrices of the parameters xx (Fig. 4,step 14), including those of the tie-points, are obtained by inverting the approximated Hessian matrix The inversion of the Hessian matrix is computationally demanding in terms of both memory and time, and it is therefore generally avoided in BA. As we are interested only in certain blocks close to the diagonal, the methods described in Rouet (2009) or Kümmerle et al. (2011) are preferable, as they allow for computation of only the required elements.
The extracted covariance matrices of the tie-points can be used to depict the 3D uncertainty, as illustrated in Fig. 1, but such a representation is less suitable for visualizing a point-cloud containing all tie-points. Therefore, we propose propagating the uncertainty from the tie-points (which can be represented by ellipsoids) to a mesh surface (Fig. 4, step 16), as illustrated in Fig. 6. The rigorous projection of 3D precision onto a surface is a non-trivial mathematical operation. Indeed, a simple approach using (for example, bilinear) interpolation of the covariance matrix elements could lead to values that differ from the truth by a factor of three, as indicated in Table A.5 of Appendix A. A preferable approach is to perform error propagation from the tie-points used as vertices in the digital surface model, as pursued in Skaloud and Schaer (2012). However, this strategy possibly underestimates the error in the face centers, and is therefore dependent on the DTM resolution. A numerical comparison employing this method indicated that the possible inaccuracy was still considerable (refer to Table A.5). Although a certain degree of imprecision owing to projection in the visualization is inevitable, using a method that considers the correct manifold of the covariance matrices, such as the log-Euclidean Riemannian metric (LERM) Arsigny et al. (2006), can maintain it at a tolerable level (for example,~5%). The details of such methodology are presented in Appendix A, while the visualization results are depicted in Fig. 6.

Precision analysis
The influences on the mapping precision are multiple. One aspect is related to errors in the observations (image, terrestrial, and aerial); another aspect is related to the prior knowledge (camera intrinsic incertitude) and global geometry (such as the number of parameters, observation redundancy, and corridor/block configuration). To detect the principal weaknesses and suggest an improvement, we propose first analyzing the situation locally at the scale of individual tie-points, and later at the level of the entire network.

Local
The predicted precision of a tie-point is provided by the covariance matrix of this tie-point, which corresponds to the appropriate lines and columns extracted from xx . This covariance is influenced by numerous factors, including the precision of the image observations, number of images in which this tie-point is observed and the precision of their aerial control, the precision of the camera intrinsic parameters, and the location of the tie-points in the images.
The local analysis does not consider the relationship between tiepoints, and in its simplest form only inspects the block diagonals of xx , as employed in the precision map creation (Section 2.8). However, the values of the off-diagonal may be significant; for example, in a steep terrain mapped with a nadir-oriented camera. In such a situation, it is recommended either to i) perform principal component analysis of the covariance matrix and display its most significant component(s), or ii) rotate the covariance matrix to align it with the terrain principal slope and aspect when displaying the diagonal terms.

Global
Larger-scale analysis is useful for identifying the principal geometrical weaknesses of the entire photogrammetric network and when searching for their effective mitigation. We first review its derivation (this theorem is trivial when V is one vector of the canonical base). If xx is the covariance matrix of the parameters in our system, the maximum weakness of this system is in the direction V , where V x T has its maximum variance; thus, when V is the eigenvector associated with the maximum eigenvalue of xx .
Analogously to the analysis of structural stiffness in civil engineering Borre (1982), Meissl (1980) and Campos (2014), this V is associated with the first vibration mode of a beam structure; that is, the vibration mode for which the structure is the most likely to oscillate.
Note that a similar analysis is used in geodetic networks Baarda (1971) and Bjerhammar (1973) and forms part of modern adjustment software Lösler (2014) and Cattin (2007), including the registration of dense networks of terrestrial laser scanning Hullo (2016).
Let max be the largest eigenvalue of xx and V max be its associated eigenvector. Thus, V max is also the eigenvector associated with the minimum eigenvalue of the Hessian H , which is the inverse of max .
The eigenvalues of xx are sorted in descending order from largest to smallest: . Thereafter, a weakness mode m i can be defined as The application of such analysis to a simulated photogrammetric network over a flat and rectangular area with sub-optimally placed GCPs (close to one corner) is depicted in Fig. 7 for the first four modes. The middle part of each sub-plot presents the undistorted situation, whereas the other two parts depict the case with ± V i i . The modes 5 and beyond are insignificant.

Placement of ground control
The previously described analysis may be used as guidance for improving the mapping precision. This can be achieved in several manners; for example, by employing sensors of higher quality, pre-calibrating several intrinsic camera parameters, adding certain GCPs or modifying the flight plan (such as increasing the lateral and inline overlap, and adding flight lines with different orientations and heights).
Whereas modifying the flight plan is a relatively rapid operation, the addition of new GCPs may be considerably more demanding in terms of time and resources. For this reason, the subsequent discussion focuses on the optimization of the GCP locations.
Intuitively, new GCPs are placed at the weakest points of the network. Thus, the local analysis is not an ideal indicator, as the precision of the tie-points may vary significantly within a neighborhood. Moreover, the tie-points close to the mapping area border are usually visible only in certain images, which is one factor that makes them less precise. In such a situation, the errors committed in the (few) observations of the GCP image coordinates may cause larger systematic deformation.
Therefore, the selection of the GCP locations is based on global analysis. In the previously mentioned analogy with a beam structure, in which the nodes represent the camera poses and the beams represent the connections between tie-points in the object and image space, the placement of a new node (GCP) should have a maximum effect on increasing the structure resistance. In our case, it should minimize its deformation owing to the accumulation of random errors. Therefore, the suggested method consists of studying the eigenvector V 1 associated with the highest eigenvalue of xx . The tie-point with the highest displacement by this eigenvector is detected and is selected as the optimal location to add a new GCP. In Fig. 7, this corresponds to adding a GCP at the point for which the deformation from the first excitation mode is maximal, yet avoiding the areas close to the border of the area of interest, and selecting only points that are visible in at least three to four images. The details are explained as part of the

Evaluation of check-point misclosures
Check-points are often used for the assessment of 3D models. However, when the check-points exhibit varying precision in different directions, which may be the case for coordinates determined by static carrier-phase GNSS on mountain slopes, their correct employment requires new considerations.
Let v be the check-point misclosure; that is, the difference between the GNSS-determined and photogrammetry-determined coordinates.
However, the previously presented method only predicts the covariance matrix cp of such a check-point for which the misclosure is v. Therefore, it is challenging to compare meaningfully the predicted covariance cp , representing the probability density, with punctual realization of the check-point misclosures v. Let dm be the covariance matrix of the direct measurement (likely obtained by static carrierphase differential GNSS 2 ).
The misclosure is a simple vector difference with compound covariance = + dif dm cp . The combination of the misclosures with their covariances in the square of the Mahalanobis norm v v T dif 1 theoretically follows the 2 distribution when the predicted precisions agree with the observed residuals.
However, if cp is underestimated (the prediction is too optimistic), the square Mahalanobis norm will be high, and if cp is overestimated (the prediction is too pessimistic), the square Mahalanobis norm will be low.
Let be the diagonal matrix containing the eigenvalues of dif 1 and V be the matrix of the associated eigenvectors, which is an orthonormal basis.
The standardized misclosures, which are defined as = v V v T , are constructed. These are unitless, with unitary standard deviation and uncorrelated coordinates. Note that the squared Mahalanobis norm of v is equal to the squared Euclidian norm of v . Thereafter, it is possible to aggregate these v , either for creating a histogram (for example, bottom plot of Fig. 15) or for proceeding with 2 tests.
For the latter, we formulated the hypothesis that the aggregated normalized misclosures follow a standard normal distribution. This hypothesis was not rejected by the one-sample Kolmogorov-Smirnov test Massey (1951) at the 5% significance level for 19 of the 26 studied cases. The majority of the problematic (rejected) cases were related to scenarios employing RTK-derived aerial control, indicating possible position biases owing to incorrect determination of (several) doubledifference carrier-phase ambiguities in several flight lines.
To proceed with the paradigm of global analysis, as presented in Section 3.2, we need to define a new check-point misclosure vector v as the concatenation of all individual check-point misclosures within one experiment. The notation v denotes the × n (3· ) 1 vector of the set of n check-points. Let cp be the theoretical covariance matrix of this v, which corresponds to the proper subset of rows and columns of the matrix xx . The principle of the subsequent analysis is to compare v to the spectral decomposition of its theoretical covariance matrix cp . Let i be the i th eigenvalue of cp when the eigenvalues are sorted in descending order and V i be the associated (normalized) eigenvector. The eigenvalue can be recovered from cp as a result of the following relation: Subsequently, we firstly have i , namely the predicted variance along the eigenvector V i , and secondly, the projection of the misclosures v on this eigenvector is V v· i . Note that the square of this projection could be interpreted as an a posteriori variance i , which is evaluated by substituting cp with the a posteriori variance of Fig. 9 depicts this principle in two dimensions. The true residual v (in red) is projected onto the eigenvectors V 1 and V 2 of its theoretical covariance matrix cp . The use of spectral analysis in the different experimental cases (Section 4) compares the actual realization of the error along an axis (v V · 1 or v V · 2 ) to the error prediction ( 1 or 2 ).

Experimental validation
The previously described methodology enables the prediction of the precision of the 3D model generated by the photogrammetric method with GNSS aerial position control under different spatial and temporal scenarios. The purpose of this section is to study the method with real data and validate it in various use cases so that the prediction corresponds to reality.

Test sites and equipment
For the experimental validation, we used small, lightweight (~1.2 kg) fixed-wing drones equipped with dual-frequency, dual-constellation (GPS & GLONASS) receivers, and small cameras with a~28 35 mm equivalent focal length. Two zones of mountainous terrain, located in the western part of Switzerland, were used for the evaluation: • Z1: A narrow valley that is surrounded by high mountain ridges in all directions but East. Numerous houses are located along the bottom of the valley. During winter 2015 to 2016 an exceptionally large avalanche fell down from the South ridge and destroyed part of the forest, the road, and several buildings. This zone first had to be mapped without GCPs (owing to the remaining risk), but repeated mapping later in summer allowed for the placement of 13 GCPs over a large part of the mapped area. The employed planes were an eBee RTK and an eBee+ from senseFly. The former carried a Cannon Power Shot S110 camera (CMOS sensor, × 4000 3000 pixels) configured for the raw image format.

Aerial position accuracy
The precision of the GNSS quality was predicted for the camera positions in Z1 using a flight plan evaluated at two different takeoff times, corresponding to "sub-optimal" and "close-to-optimal" satellite visibility conditions. This evaluation employed the satellite almanac and coarse DEM, the results of which are displayed in Fig. 8 using the color code of Fig. 5.
The actual 3D positioning quality obtained in PPK is illustrated in Fig. 10 in terms of an internal characteristic factor varying from 1 to 5 (best to worst). In agreement with the prediction, this indicates the decimeter level and centimeter level for the "sub-optimal" and 2 The positions of the GCPs were measured by terrestrial GNSS with particular care: each point was measured for at least min 10 , leading to a standard deviation of cm2 in planimetry and cm3 in altimetry. The covariance matrix of the ground measurement is Cledat, et al. ISPRS Journal of Photogrammetry and Remote Sensing 167 (2020) 24-38 "optimal" flight times, respectively. Note that, although the 2 nd flight exhibited better satellite configuration over most of the zone, it was shortened towards the South to maintain aircraft navigation safety, as critical satellite visibility was predicted at the final line. Further analysis involving other flights in the same zone compared the differences between the photograph positions obtained by GNSS (RTK/PPK) with those of the indirect orientation using all GCPs. These airborne position misclosures are presented in the plots of Fig. 11 for RTK and Fig. 12 for PPK. Note that the vertical component was generally worse than the horizontal one, as GNSS positioning is better in planimetry than in altimetry.
The correspondences between the residuals and posterior precision were quite strong for PPK, considering that the quality of the camera position derived via indirect positioning (the given precision of the camera pose center from the indirect orientation was considered in this statistical analysis) was not necessarily superior to that of GNSS and remained correlated with the other parameters (such as the attitude & interior orientation). This may explain certain residuals exceeding the normal probability levels.
Although such excess was less present for the RTK positioning, the precision of this positioning mode was lower than that of PPK by a factor of 10 in numerous photographs. This suggests that caution should be exercised when employing RTK technology for achieving centimeterlevel positioning in the airborne environment with frequent discontinuities in satellite tracking.

Tie-point density
It is interesting to note that, with the exception of (close to) uniform texture (for example, dense and high vegetation; water), the predicted mapping precision appeared not to be very sensitive with respect to the number of tie-points, as observed empirically together with its expected precision. For example, an increase or reduction by a factor of two of the tie-point density led to a difference of several percent in the predicted mapping precision. This follows from Gruber's rule, which states that only five points are sufficient for the relative orientation of two poses of a calibrated camera Nistér (2004). The number of points is higher for small, non-metric cameras, but several tens of well-distributed points are sufficient. Thus, the benefit of the added points exhibits asymptotic behavior. Although the selected inter-alpine environment of zones Z1 and Z2 included textural changes, from scree and grass to low and high vegetation, the observed variations in the point densities did not affect the orientation process. However, in certain areas with high vegetation, surface restitution was not possible (for example, the black areas within the internal part of Fig. 13).

Local analysis
This subsection studies the mapping prediction of different flight configurations in both zones, namely Z1 and Z2. In total, 26 scenarios, as described in Table 2, were derived from flights conducted above these two areas. In the corridor configuration, only two parallel flight lines were considered, where the so-called mono-block represents the "lawnmower sweeping" pattern with flying height adaptation per flight line, as illustrated in Fig. 13. A stair plan is a particular flight plan that staggers several smaller blocks. While the flying altitude remains constant within each block, it changes in a step among them. This plan is useful when mapping a complicated terrain shape, such as that of Z2 (that is, a narrow and steeply climbing valley).
In total, 13 control points were determined in Z1 and 20 were determined in Z2, which enabled different GCP/CP ratios and configurations to be considered. In the case of no GCP, all signalized points were considered as CPs. In the so-called bad GCP configuration, the barycenter of the GCPs differed substantially from the barycenter of the mapped area, while the remaining points were considered as CPs. However, the good GCP configuration considered the spread of GCPs over the entire terrain. To ensure that the evaluation was meaningful, the number of CPs was at least 5.
For the aerial control, standard GNSS positioning (~meter-level accuracy) was considered in the standalone case, whereas carrier-phase corrections were employed in the RTK and PPK, with the aim of In the specific case 13 (Fig. 14) as well as in all aggregated cases (Fig. 15) This demonstrates that the conducted experiments respected the basic properties of a normal distribution. However, global analysis needed to be performed to confirm such an assumption by means of additional statistical testing.

Global analysis
The comparison for case 1 performed in Z1 is detailed in Fig. 16. The prediction (indicated by the blue line) is the standard deviation evaluated along the eigenvector of the covariance matrix cp , computed as the square root of cp eigenvalues. The realization (indicated by the red line) is the projection of v onto the eigenvectors: v V · i . The predicted trend of decreasing the eigenvalues of cp was followed by the projection of the observed misclosures onto the corresponding eigenvectors. In terms of the deformation theory presented in Section 3.2, the misclosures were a linear combination of different deformation modes (schematically represented in Fig. 7), and in the studied case, the obtained coefficients were of the same order of magnitude as those calculated by means of the prediction.
As a similar agreement was observed in the other studied configurations, the observations were coherent with the predicted precision.  E. Cledat, et al. ISPRS Journal of Photogrammetry and Remote Sensing 167 (2020) 24-38 4.6. New GCP placement For the practical demonstration of the method relating to the optimal placement of the next GCPs (as presented in Section 3.3), case 5 in Z1 was used, where images of an inclined terrain in block configuration were oriented via three (sub-optimally placed) GCPs (Fig. 17).
In this case, we were searching for the best place to situate an additional GCP. For this purpose, we employed the following quality indicators (a b c , , , and d) and investigated different methods for the tiepoint precision aggregation. a assess the tie-point with maximum 3D standard-deviation. b is the root mean squared (RMS) of all tie-points standard-deviation, c is the mean of tie-points 3D standard-deviation and d is the largest eigen-value of the full covariance matrix of the tiepoints (and thus the highest mode).
We evaluated the quality indicators a b c , , , and d 13 with the three GCPs and placed them in the 1 st column of Table 3.
Firstly, we searched for a tie-point with the maximum incertitude (indicator a) and added the 4 th GCP exactly at this location. The coordinates (with respect to the tie-point centroid) of this new GCP are indicated in the first two lines of Table 3 and the updated quality indicators are displayed in the column "4 GCP f a [] ( )". As we will observe later, this strategy of selecting an adequate position for the 4 th GCP was not appropriate. Secondly, we found a tie-point with the maximum value in the eigenvector associated with the maximum eigenvalue of the global xx , and we placed the 4 th GCP at this point. The coordinates of this point and the quality indicators computed using this new GCP are presented in the final column of Table 3.
Thereafter, we performed a type of benchmarking to validate the suggested methods. We created a × 4 4 m grid over the mapped surface and placed the 4 th GCP at each node of this grid, evaluating the abovementioned quality indicators every time. We selected the grid place by minimizing a b c , , , and d separately and its coordinates are displayed in Table 4.
It can be observed from this analysis that, even if the quality indicator a is the most intuitive and straightforward, it should not be considered owing to its sensitivity to the random distribution of the tiepoints. Hence, we suggest focusing on the quality indicators b c , , and d. To compare the quality indicators computed with three GCPs and the same quality factor computed with four GCPs, we could define the following improvement factor: where the empty white square symbol: could be replaced by any of the quality indicators a b c , , or d.
GCPs 3 is the indicator calculated for three GCPs, whereas min is the indicator calculated for the optimal position of the new GCP (in bold in Table 4). This ratio is defined as null if the new GCP is not usable; for example, if it does not appear in at least two images, and 100% if the GCP placement is at the best selection, as evaluated by the (time-consuming) grid search.
For any quality indicators b c , or d, this improvement factor is approximately 45% for the second column of Table 3; that is, if a new GCP is placed on the tie-point with the maximum standard deviation. However, this improvement factor is approximately 95% for the third column of Table 3; that is, if 4 th GCP is placed on the highest value of the eigenvector associated with the highest eigenvalue of xx . The quality evaluation that was performed using this exhaustive search demonstrated practically that the method for placing a new GCP, as described in Section 3.3, is nearly optimal (at a 95% level) with respect to the considered global mapping quality indicator b c , or d.

Conclusions and perspectives
We have proposed a method for assessing the quality of drone mapping in a complex natural environment prior to the mapping being performed. The method is based on a flight plan and consists of simulating all of the observables required in the real mapping adjustment, considering a particular time and date, the observations, and the geometry of the photogrammetric network, using a priori knowledge regarding the terrain, drone payload, and GNSS satellites. Compared to the traditional approach, its application avoids the need for mission repetition frequency, which increases with the terrain complexity ( Fig. 3 left vs. right).
The comparison of the quality prediction with the actual mapping accuracy in various geometrical configurations as well as the quality of the airborne GNSS positioning in the mountain environment exhibited satisfactory agreement in the sense that the statistical testing did not reject the hypothesis (at 5% significance) of the observed misclosures fitting the prediction.
Practical experience also demonstrated that the PPK approach is preferable over the RTK technology in an environment in which frequent occlusion of satellite signals occurs owing to either the drone motion or its surroundings. For the sake of navigation safety, the method of predicting the satellite positioning quality is worth considering not only at camera stations, but also (possibly with different criteria) over the entire drone trajectory, including the takeoff and landing zones.
The signalization and surveying of GCPs requires substantial effort in drone mapping, particularly in a complex terrain. Therefore, a method for the most effective placement of new GCPs was designed and evaluated. When actual realization of GCPs in the optimal location is not practical (for example, when the environment is too steep or is not accessible), the user may alter other geometrical or temporal parameters of the flight, and eventually change the equipment to maximize the likelihood that the executed mission will satisfy the desired precision criteria. Therefore, it is important that the prediction is not only correct, but also is executed rapidly and can offer intuitive interpretation. In our experience, the adopted metric of covariance interpolation and visualization satisfies this requirement.
In the future, we plan to develop an algorithm that can provide improved determination of the tie-point density and tie-point quality based on the (a priori) ground texture. This algorithm could make use of modern computer vision and machine learning. Moreover, the precision prediction could be extended to optimize the entire flight plan and overall GCP placement.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgement
This work was supported by Innosuisse Grant No. 18442.1 PFIW-IW. This support is greatly appreciated.   Appendix A. Interpolation of estimated precision

A.1. Bilinear
The first (and naive) method for computing the precision of any point on a surface is to proceed to bilinear interpolation of the covariance matrix. Let P P , a b , and P c be three vertices of a given triangle of the Delaunay triangulation and let P d be a point inside this triangle. Then, there exist and in [0, 1] with + < 1, such that = + + P P P P (1 ) d a b c . The bilinear interpolation leads to the following: However, it is not mathematically correct to interpolate covariance matrices element wise. Skaloud and Schaer (2012) suggested performing error propagation from the three observed vertices of the triangle to the point P d . Below is a modified version of the presented error propagation, taking into account the eventual correlations between the different points in a rigorous, fully 3D calculation. This algorithm was designed for estimating the precision of a grid cell center from a relatively regular and dense laser-scanning pattern. It is less suitable for randomly (and more sparsely) distributed tie-points, because it may overestimate the precision towards the middle of the triangles (when P d does not necessarily lie in the perfect plane defined by P P , a b , and P c ).

A.3. Log interpolation
As the covariance matrices are not directly associated with the conventional metric of Euclidean space, it is improper to proceed with linear combinations of their elements as previously suggested, although this appears to be practical. However, the logarithms of such matrices are part of Euclidean space Arsigny et al. (2006). Let the log of a covariance matrix be the matrix in which the eigenvalues are substituted by their logarithms. If = V V T , where is the diagonal matrix of the eigenvalues, and V is the matrix of the eigenvectors, we can define log ( ) as the diagonal matrix containing the logarithms of the diagonal elements of , and = log V log V ( ) ( ) T . Moreover, it is possible to demonstrate that = exp log ( ( )) using the power series definition of the exponential function. Therefore, the suggested interpolation of the covariance matrix logarithm is expressed as follows: .

A.4. Comparison
To assess the suitability of the methods for obtaining dd , we execute the algorithm that evaluates the covariance matrices of the tie-points. Thereafter, we consider three precision maps for display purposes, which are created by the three methods described above.
We run the code a second time to obtain different tie-point configurations, each with its predicted covariance matrix. This second prediction enables validation of the interpolation method selected for the first prediction. For each tie-point of the second prediction, the covariance matrix will be compared to the value provided by the interpolated map based on the first prediction. There are several manners in which the × 3 3 covariance matrices can be compared. One of these is the relative LERM Arsigny et al. (2006). The relative LERM of the matrices 1 and 2 is expressed by the formula below, where • F is the Frobenius norm (the square root of the sum of the squared elements of the matrix). It is possible to compute the relative horizontal error and relative vertical error, which are defined by the following two equations.

11%
31% 6% To aggregate the statistics computed for all tie-points, we compute the maximum difference, the mean of these differences, and the square root of the mean of the squared differences in Table A.5. In both criteria presented in this table, it can be observed that the logarithmic interpolation is the most appropriate method for our problem, as it allows for interpolating the covariance matrix at any position, a fortiori, with the smallest error.