IMPROVING UAV TELEMETRY POSITIONING FOR DIRECT PHOTOGRAMMETRY

: A fundamental step of UAV photogrammetric processes is to collect Ground Control Points (GCPs) by means of geodetic-quality GNSS receivers or total stations, thus obtaining an absolutely oriented model with a centimetric accuracy. This procedure is usually time-consuming, expensive and potentially dangerous for operators who sometimes need to reach inaccessible areas. UAVs equipped with low-cost GNSS/IMU sensors can provide information about position and attitude of the images. This telemetry information is not enough for a photogrammetric restitution with a centimetric accuracy, but it can be usefully exploited when a lower accuracy is required. The algorithm proposed in this paper aims at improving the quality of this information, in order to introduce it into a direct-photogrammetric process, without collecting GCPs. In particular, the estimation of an optimal trajectory is obtained by combining the camera positions derived from UAV telemetry and from the relative orientation of the acquired images, by means of a least squares adjustment. Then, the resulting trajectory is used as a direct observation of the camera positions into a commercial software, thus replacing the information of GCPs. The algorithm has been tested on different datasets, comparing the classical photogrammetric solution (with GCPs) with the proposed one. These case-studies showed that using the improved trajectory as input to the commercial software (without GCPs) the reconstruction of the three-dimensional model can be improved with respect to the solution computed by using the UAV raw telemetry only.


INTRODUCTION
Unmanned Aerial Vehicles (UAVs) are nowadays a standard tool for photogrammetric surveys, with applications in e.g. mapping, land and building 3D modelling, and hazard monitoring (see e.g. Eisenbeiß, 2009;Rosnell et al., 2011;Lucieer et al., 2014;Nex and Remondino, 2014;Forlani et al., 2015;Ryan et al., 2015;Skarlatos et al., 2015;Turner et al., 2015;Pagliari et al., 2017;Avanzi et al., 2018). These applications generally rely on the collection of a set of Ground Control Points (GCPs), that is required to obtain an absolute orientation of the photogrammetric block and also to achieve centimetric accuracy in the Bundle Block Adjustment (BBA) (Grenzdörffer et al., 2008;Nocerino et al., 2013;Gini et al., 2012). In this procedure, tie points are generally automatically detected by the processing software and, therefore, no additional information is strictly required (Remondino et al., 2011). Despite this, the acquisition of GCPs, usually performed by a (geodetic quality) GNSS receiver or a total station, is a time consuming and expensive activity. Moreover, the operator collecting those points can be exposed to risks, especially in emergency or hazardous contexts .
In addition to the camera, the UAVs on the market that are currently used in photogrammetric applications commonly have on-board other low-cost sensors for the determination of the UAV position and attitude, such as GNSS receivers, barometers, radars, gyros, accelerometers, compasses. These sensors are mainly used for piloting purposes, allowing UAVs to travel along an a-priori designed path. The data collected by these sensors are usually recorded in the telemetry files (Daakir et al., 2017). These telemetry data can be exploited as an additional information into the BBA. However, due to the use of low-cost sensors, the quality of these data is generally quite * Corresponding author poor and, therefore, the telemetry cannot be efficiently used to reduce or eliminate the GCP information to get an absolute orientation and a 3D point cloud with centimetric or decimetric accuracy, unless real time kinematic (RTK) or post processed kinematic (PPK) capable GNSS is installed on-board (Rizaldy and Firdaus, 2012;Chiabrando et al., 2013;Eling et al., 2015;Albéri et al., 2017;Mian et al., 2015Mian et al., , 2016Wierzbicki, 2018). However, only few UAV models on the market currently include RTK or PPK receivers, making them more expensive.
To overcome this limitation, in the present work we implement a preprocessing strategy to improve the UAV positioning by combining the telemetry data with the camera positions computed by a relative orientation of the acquired block of images. This improved trajectory is then used as an input to the BBA, allowing to obtain a 3D model that differs at maximum by few tens of centimetres from the same model computed by also using GCPs. These statistics do not include bias effects due to the common use of stand-alone GNSS solutions in the UAV telemetry.

THE PROPOSED METHOD
The telemetry preprocessing algorithm takes as input the UAV position observed by the on-board sensors (usually a GNSS receiver combined with a barometer) and the purely photogrammetric position of each acquired image i. This photogrammetric position is determined by BBA from the set of collected images and consists in the coordinates of the image projection centres in an arbitrary reference system (relative orientation). It can be computed by using commercial software, e.g. Agisoft Metashape Professional, without providing any GCP. Hereafter, we will refer to the two coordinate reference systems as absolute (a) and relative (r), respectively. Therefore, for a generic image i, the two "observed" quantities are given by the following equations: i is the observed position by the GNSS telemetry in the absolute-coordinate reference system, while ∆r P,(r) i is the vector joining two subsequent images i and i + 1 in the relative-coordinate reference system, computed from the purely photogrammetric trajectory estimated by the commercial software, i.e. ∆r P,(r) . The unknowns are the camera projection centre positions of two consecutive images i and i + 1 in the absolute-coordinate reference system, defined by the vectors r (a) i and r (a) i+1 ; the parameters of the rotation matrix R, e.g. the three Cardan angles ω, ϕ, and κ, and the scale factor k defining the transformation between the relative and the absolute reference systems. Therefore, considering a set of N images, the total number of unknowns is 3N + 4, where the first term represents the three components of the N position vectors r (a) i and the latter the three Cardan angles ω, ϕ, κ and the scale factor k. Note that, although the purely-relativephotogrammetric trajectory is available, to reduce the impact of possible drifts it is convenient writing Eq. 2 in terms of the position difference between subsequent images. Moreover, the difference between the camera projection centre and the point to which the telemetry is referred to (e.g. the phase centre of the GNSS antenna) is neglected because usually smaller than the accuracy of the telemetry itself.
To manage the non-linearity with respect to the Cardan angles ω, ϕ, and κ defining the rotation matrix R, we assume to know an approximated matrixR. The latter can be reckoned, together with the translation vectort and the scale factork, solving a least-squares system of equations modelling the Helmert transformation between the same set of points known in the absolute and relative reference systems by the single value decomposition (SVD) (Arun et al., 1987). This approach allows a very fast and non-iterative solution, i.e. no approximated values of the estimated parameters are required, but has the limit of disregarding the weights of the single observations, that in our case are fundamental. That is why we consider this solution as approximated apart from the scale factor k, that is fixed atk and not re-estimated. Therefore, given the approximated rotation matrixR, the matrix R can be written as: IfR is a good approximation of R, the three Cardan angles composing the δR matrix are small angles. Therefore, δR can be written as (Kraus, 1997): where δω, δϕ and δκ are the three small Cardan angles describing the correction to the approximated anglesω,φ andκ composing theR matrix.
Finally, according to Eqs. 3 and 4, Eq. 2 can be rewritten as: Recalling that the inverse of a rotation matrix is equal to its transpose, i.e.
finally, the optimal combined trajectory can be reckoned by solving the linearized least-squares adjustment including Eqs. 1 and 5 for all the N images of the block, once the proper variancecovariance matrix is assigned to the observation noises. Therefore, the observation vector y o becomes: As for the telemetry error, its variance-covariance matrix Cν i can be directly provided by the on-board sensors or can be empirically calibrated from a static acquisition before starting the flight. When empirically calibrated, we neglect the correlation between the three components of the observed position vector , thus the resulting covariance matrix is: As for the photogrammetric error, we first define: and then derive its variance-covariance matrix by applying the covariance propagation law as: As already mentioned, the photogrammetric observations are determined by differentiating the photogrammetric trajectory, computed by relative orientation. Therefore, to define the error covariance matrix Cη i of this kind of observations, it is required to propagate the covariance matrix of the position of the projection centre of each image C r P,(r) i , which is usually provided by commercial software. Applying again the covariance propagation law, we get: where I is a 3 × 3 identity matrix. Notice that Eq. 11 disregards the correlation between orientation parameters of different images, i.e. there are zeros outside the diagonal block, because this information is rarely provided by commercial software. Furthermore, considering only a couple of images at each time, the correlation introduced by using the photogrammetric position of one image in more than one observation is neglected, i.e. the fact that the observed relative position of the image i is used to compute both ∆r P,(r) i−1 and ∆r P,(r) i is neglected in terms of covariance propagation.
Introducing Eq. 11 into Eq. 10, we obtain: According to Eqs. 8 and 12 the cofactor matrix, used to weight the observations in the least-squares adjustment, can be finally written as: As we will see later, the values of the GNSS telemetry and photogrammetric variances and covariances composing the cofactor matrix could lead to statistical inconsistencies between the two solutions, especially in the vertical component. This is the case if drifts are present in the photogrammetric trajectory. A possible solution is to model the drift in the least-squares adjustment. However, no information on the drift is usually available. Therefore, assuming that the drift is not too large, a possible solution is to rescale by a factor λ the terms associated to the vertical component in the relative-photogrammetric-trajectory covariance matrix that is used to define Cε i in Eq. 12, i.e. to define a modified covariance matrix CrP,(a) i as: and the vectorr P,(a) i is the purely photogrammetric trajectory expressed in the absolute reference system by means of the approximated Helmert transform, i.e.: The factor λ can be determined by setting a statistical test to check whether the vertical coordinate of the GNSS telemetry trajectory z G,(a) i and the one of the photogrammetric trajectorỹ z P,(a) i have the same mean value for each image i: are the variances of the two vertical coordinates of the image i, taken from the Cν i and the CrP,(a) i matrices, respectively, and Z is a standardized normal distribution. Assuming a significance level α of the test, we can compute the threshold Zα used to determine whether the initial hypothesis is accepted or not. In other words, we are testing if the two vertical coordinates are statistically consistent with each other according to their standard deviations. In order to calibrate the value of the λ parameter, we assume that it cannot be smaller than 1 and that it has to be large enough to verify the hypothesis, i.e. to verify thatZi(λ) ≤ Zα, for at least (1 − α) N images, where N is the total number of images.
Once this λ covariance scale factor is empirically tuned, the modified version of Q, hereinafter Q , can be determined according to Eqs. 13 and 14 and introduced in the non-linear leastsquares adjustment solving Eqs. 1 and 5 for all the N images.
The approximated values of the absolute trajectory coordinates can be directly taken from the GNSS telemetry, while the three residual Cardan angles δω, δϕ and δκ can be initialized at 0. The solution of this system is the improved trajectory of the UAV, i.e. a set ofr (a) i vectors together with their covariance matrix Cr(a) i , that can be used as input in commercial software to perform a direct-photogrammetric solution.
A remark is worthwhile before concluding this section. In principle, the optimal exploitation of the telemetry trajectory can be obtained by introducing the corresponding observation equations into the BBA, as it is done e.g. in the CALGE software (Forlani, 1986), without the needing of passing through the trajectory derived from a relative orientation. However, this intermediate step offers the possibility of correcting possible drifts or other systematic errors that cannot be easily detected just looking at the input tie points. Moreover, the goal of this work is not to find the optimal combination of the data, but rather to define a feasible preprocessing procedure to improve the telemetry trajectory of a photogrammetric survey and the corresponding 3D restitution obtained from a (black box) commercial software, in this case Agisoft Metashape Professional.

EXPERIMENTAL CASE STUDIES
The proposed method was tested on the field to assess its performance. To this aim two surveys were performed in the following areas: (1) the bank between the Muzza Channel and the Adda River in Cassano D'Adda (Milan, Italy) and (2) the Rossia flight field in Gossolengo (Piacenza, Italy). The results of these two experiments will be presented and commented in the next subsections.

Cassano d'Adda experiment
The photogrammetric survey of the bank between the Muzza Channel and the Adda River in Cassano D'Adda (Milan, Italy) was performed by using a DJI Matrice 210 quadcopter, mounting a DJI Zenmuse X5S camera with a micro 4/3 sensor and a DJI MFT 15 mm f 1.7 ASPH lens on a 3-axis gimbal. The mean flight altitude was 50 m and the mean overlapping was 70%, acquiring about 380 photos over an area of 350 m × 150 m. A set of 17 GCPs, homogeneously spread over the surveying area according to natural obstacles, were acquired by a geodetic GNSS receiver working in network real time kinematic (NRTK) mode, leading to an accuracy of few centimetres in the positioning. First of all, a reference solution including also GCPs was computed by BBA through the Agisoft Metashape Professional software. At this stage, the images are processed at full resolution (High Accuracy option in the software) to find the key points used by the automatic matching algorithm to detect tie points (TPs), and the image coordinates of GCPs were manually collimated. After the BBA, also the 3D dense point cloud was computed by exploiting the full resolution images, see Figure 1. As for the camera calibration, a self-calibration based on the available GCPs was directly applied inside the Metashape software, while performing the BBA. Then, this reference calibration was fixed throughout the further steps of the experiment, to reduce the impact of the camera calibration in the comparisons. Once the reference solution was obtained, it was firstly used to assess the telemetry data quality. The histograms and the statistics of the coordinate residuals are reported in Figure 2 and Table 1, respectively. Thanks to the on-board barometer, the accuracy was higher in the vertical direction (with an empirical standard deviation of about half a meter) than in the horizontal one (with an empirical standard deviation of few meters). As for the angular observations measured by the on-board gyros, the empirical standard deviation with respect to the reference   solution is in the order of 1 • with not negligible biases (up to 5 • ). For this reason, they are not taken into account in this investigation. These empirical accuracies are then used in the subsequent steps, e.g. to define the Cν i matrix (see Eq. 8), since the considered UAV does not provide information about the accuracy of the data recorded in the telemetry file.
Before proceeding to the telemetry preprocess, the consistency between the telemetry and the photogrammetric-relative trajectory was checked by applying the testing procedure shown in Eq. 17. The telemetry and photogrammetric-relative (converted in the absolute reference system by means of the approximated Helmert transform, see Eq. 16) trajectories are shown in Figure 3. Here, a drift in the vertical component of the photogrammetric-relative solution is clearly visible. Therefore, making the two trajectories statistically compatible required a scale factorλ = 20, for a test significance level α = 5%. Figure 3. Telemetry trajectory (blue) and photogrammetric-relative trajectory (red). The latter was converted in the absolute reference system by means of the approximated Helmert transformation. The vertical axis is 4 times scaled with respect to the horizontal ones.
The UAV telemetry was then preprocessed before the BBA, by implementing the method described in Section 2. Here, the error variances of the telemetry coordinates were considered according to the results in Figure 2, while the error covariances of the photogrammetric-relative solution were taken from the Metashape software, after applying the estimated scale factorλ. The cofactor matrix Q was then obtained according to Eqs. 8 and 14. Table 2 shows the statistics of the differences between the image positions estimated by the preprocessed telemetry and those computed in the reference solution. There is a clear improvement in terms of standard deviation of the preprocessed trajectory when comparing it with the one observed by the onboard sensors (see also statistics in Tables 1). On the other hand, the mean error cannot be reduced due to the quality of the onboard GNSS receiver and the stand-alone method in processing its data. This preprocessed trajectory, namely the set ofri vectors, together with its estimation error, namely the set of Cr i matrices, was introduced again into the Metashape software to compute the 3D point cloud without using GCPs in the BBA. Hereafter, this solution is referred to as the "adjusted" one.
In order to evaluate the improvement carried by the proposed preprocessing algorithm into the "adjusted" solution, other two alternative solutions were computed. Firstly, a 3D point cloud was computed by the Metashape software using as input only the original telemetry (without GCPs) in the BBA. The telemetry was weighted according to the standard deviations of Table 1. Hereafter, this solution is referred to as the "telemetry" one. Secondly, a 3D point cloud was computed by the Metashape software without any telemetry and GCPs in the BBA, i.e. by simply applying the approach used to compute the photogrammetric-relative trajectory. To correctly scale and locate the resulting 3D model, the transformation described by Eq. 16 was applied. Hereafter, this solution is referred to as the "relative" one.
The "adjusted", "telemetry", and "relative" solutions, i.e. the three direct-photogrammetric solutions, were compared with the reference one. For the first comparison the GCPs were used as control points (CPs) in the three direct-photogrammetric solutions. In particular, we computed the residuals between their estimated coordinates and the ones observed by the geodetic GNSS receiver.  Table 4. Statistics of the differences between the image centres of projection estimated by the "adjusted" (Adj), "telemetry" (Tel), and "relative" (Ref) solutions and the ones of the reference solution.
the reference one. However, as expected, it was not possible to remove the bias related to the low quality of the on-board GNSS receiver.
The second comparison is about the quality of the estimated external orientation parameters by the BBA. The focus is on the estimated image centres of projection that were compared with those in the reference solution. Table 4 shows that the "adjusted" solution again introduces an improvement with respect to the "telemetry" and "relative" solutions, reducing the standard deviations of the residuals. A bias is present, due to the fact that the "absolute" georeferencing procedure relies on the low-cost on-board GNSS receiver.
Finally, the computed 3D point clouds were compared in terms of distances with respect to the reference one. For each point of the direct-photogrammetric point clouds, the Euclidian distance with respect to the closest point in the reference cloud was computed. To correctly perform this comparison, the bias in the georeferencing was removed, by translating the point clouds by the mean value of the CP residuals reported in Table 3. Figure 4 shows that the "adjusted" solution leads to a big improvement "adjusted" "trajectory" "relative" Figure 4. Cloud-to-cloud distances of the "adjusted", "telemetry", and "relative" point clouds with respect to the reference one. On the right side of each colourbar, the histogram showing the distance distribution is depicted. The colourbar is saturated to 1.2 m (disregarding less than 1%, 3%, and 15% of points, respectively) The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences, Volume XLIII- B2-2021XXIV ISPRS Congress (2021 with respect to the other two. In fact, here more than 68% of the points has a distance below 10 cm and more than 96% of the points below 30 cm, with a mean distance of about 9 cm, thus confirming a relative geometrical accuracy at the decimetre level. This is not true for the other two cases, where only more than 12% ("telemetry" solution) and more than 32% ("relative" solution) of the points have a distance below 30 cm, and the mean distance is always greater than 60 cm.

Rossia experiment
To survey the Rossia flight field in Gossolengo (Piacenza, Italy) a flight with a DJI Matrice 210 v2 quadcopter was performed, mounting a DJI Zenmuse X5S camera with a micro 4/3 sensor and a DJI MFT 15 mm f 1.7 ASPH lens on a 3-axis gimbal. The flight was designed in a double grid configuration, with the mean flight altitude of 35 m and the mean overlapping of 60% in the longitudinal direction and 70% in the transversal one. About 160 photos were acquired over an area of 150 m×120 m and 10 GCPs, homogeneously spread over the surveying area, were acquired by a Leica MS60 multistation , leading to an accuracy better than a centimetre (see Figure 5). A reference solution was computed by the Metashape software, exploiting full resolution images in both the BBA and 3D dense point cloud determination, applying the same procedure of the Cassano d'Adda experiment (see Section 3.1), except for the self-calibration. In this case, the camera was calibrated by an independent block of images on a smaller area (about 20 m × 20 m) close to the take-off point. Here, another set of 12 GCPs placed on a regular grid was measured by the Leica MS60 multistation, obtaining millimetric accuracy. Then, a set of photos with different heights and camera attitudes were taken and processed by the Metashape software, exploiting the self-calibration feature based on GCPs. The camera calibration parameters computed at this stage were fixed for all the other solutions computed throughout this experiment (including the reference one). Notice that, the proposed configuration reproduces directly on the field the classical camera calibration by a checkerboard.
Also in this experiment, the UAV did not provide information about the telemetry accuracy. Therefore, as a first step, the recorded telemetry was compared with the trajectory obtained by the reference solution. The statistics of this comparison are reported in Table 5, showing a better quality of the on-board instruments (in terms of standard deviation) than the one of the UAV used for the previous experiment.
East [m]  It has to be stressed that this approach to determine the telemetry accuracy is not applicable from the practical point of view. In fact, if a solution computed by GCPs is available, there is no need for a direct-photogrammetric solution. Therefore, another approach was adopted. Before performing the flight, the UAV was kept still at the same ground position for few minutes, once the telemetry logging was activated. Evaluating the variability of the recorded position around its empirical mean value, the average accuracy of the telemetry for the whole flight can be defined. These estimates of the accuracy are reported in Table 6. If compared with the ones of  Table 6. Standard deviations of the UAV position from telemetry, holding the UAV on ground for three minutes before the take-off.
According to the proposed algorithm, the telemetry trajectory was preprocessed, showing a quality improvement with respect to the raw UAV telemetry, as it can be seen by comparing the standard deviations of Tables 5 and 7. In particular, the preprocessing algorithm is able to halve the standard deviations of the estimated trajectory with respect to the reference one, while it is not able to reduce the biases due to the absolute positioning performed by the on-board GNSS only.  Table 7. Statistics of the differences between the image positions from the preprocessed telemetry and those estimated in the reference solution.
Then, the preprocessed trajectory was used as input to compute the 3D point cloud exploiting a BBA without any GCP in the Metashape software. This solution will be referred to as "adjusted". Also in this experiment, a "telemetry" solution was computed, i.e. a 3D point cloud exploiting a BBA without any GCP, providing the raw UAV telemetry instead of the preprocessed one as input to the Metashape software. These two solutions were compared to the reference one in terms of differences in the estimated coordinates of the GCPs (that were just used as CPs in the BBA) and in terms of point cloud differences, applying the comparison methodology explained in Section 3.1.
The comparisons in terms of residuals on CPs (see Table 8) show that the "adjusted" and "telemetry" solutions have similar standard deviations at the decimetre level and that, as expected, it was not possible to remove the bias related to the low quality of the on-board GNSS receiver.
The comparisons in terms of point cloud distances (see Figure 6) confirm that the two direct-photogrammetric solutions have the same quality, reaching the decimetre level of accuracy with more than 80% of the points with a distance smaller than 30 cm, without considering the already mentioned bias due to the low quality of the on-board GNSS receiver.
These results show that the two direct-photogrammetric solutions are almost equivalent in terms of estimated 3D model, reaching an accuracy at the decimetre level. This equivalency is mainly caused by the UAV raw telemetry quality that in the current experiment is better than half a meter, allowing a decimetric accuracy in the 3D model reconstruction, if correctly weighted inside commercial software packages (here Agisoft Metashape Professional).  Table 8. Statistics of the differences between the GCPs/CPs coordinates estimated by the reference (Ref), "adjusted" (Adj), and "telemetry" (Tel) solutions and the ones observed by the Leica MS60 multistation.
"adjusted" "trajectory" Figure 6. Cloud-to-cloud distances of the "adjusted" and "telemetry" point clouds with respect to the reference one. On the right side of each colourbar, the histogram showing the distance distribution is depicted. The colourbar is saturated to 1.2 m (disregarding less than 1% of the points in both cases).

CONCLUSIONS
In this work, a preprocessing algorithm for the UAV raw telemetry was proposed. The aim is to improve the telemetry quality by combining it with a purely photogrammetric trajectory obtained by processing the acquired images without the using GCPs. This preprocessed trajectory was used as input for a direct-photogrammetric BBA into commercial software packages, instead of the raw one.
The two performed experiments showed that this approach leads to a general improvement of the trajectory accuracy, up to the decimetre level. As for the quality of 3D geometrical models, e.g. point clouds, determined by commercial software packages exploiting the improved telemetry, a relative decimetric accuracy becomes achievable, also correcting the effect of photogrammetric drifts, usually present when GCPs are not considered. Nevertheless, if the quality of the raw UAV telemetry is already at the level of few decimetres, the improvement could not be so significant.
The developed algorithm can be very useful especially when the GCP collection is dangerous, e.g. in emergency and risky scenarios, too expensive, or when centimetric accuracy is not really required, e.g. in hydraulic modelling. It is a cheaper alternative to the use of RTK/PPK capable GNSS receivers on-board the UAV, that usually increase the cost of the instrumentation.
Finally, it has to be underlined that the proposed method does not work in case of linear acquisitions because of rotation ambiguities. Since linear acquisitions are sometimes performed in UAV photogrammetry, this is a drawback to be overcome. In this respect, a possible future improvement of the developed algorithm consists in also estimating angular corrections to the telemetry data by using the information from the relative photogrammetric solution.