Error Reduction in Vision-Based Multirotor Landing System

New applications are continuously appearing with drones as protagonists, but all of them share an essential critical maneuver—landing. New application requirements have led the study of novel landing strategies, in which vision systems have played and continue to play a key role. Generally, the new applications use the control and navigation systems embedded in the aircraft. However, the internal dynamics of these systems, initially focused on other tasks such as the smoothing trajectories between different waypoints, can trigger undesired behaviors. In this paper, we propose a landing system based on monocular vision and navigation information to estimate the helipad global position. In addition, the global estimation system includes a position error correction module by cylinder space transformation and a filtering system with a sliding window. To conclude, the landing system is evaluated with three quality metrics, showing how the proposed correction system together with stationary filtering improves the raw landing system.


Introduction
The growing demand for drone applications motivates the study of the support technologies for this type of small and powerful unmanned aerial vehicle (UAV). However, all new applications share an essential and critical maneuver-landing.
Generally, work in the literature about landing maneuvers, both for fixed and rotary wing UAVs, focuses on control strategies [1][2][3]. All of them require access to the internal vehicle states, the actuators or specific modes of the control or navigation system.
Under the precision landing concept is included all the solutions that approach this maneuver in an autonomous or supervised way, independent of the techniques and sensors used to estimate the vehicle states such as position or orientation, as well as its corresponding velocities and accelerations. The landing maneuver can be included within the navigation system where two main groups can be distinguished, outdoor and indoor navigation. Generally, outdoor navigation is based on the Global Navigation Satellite System (GNSS), but practically all current systems use fusion techniques that allow the integration of different strategies for estimating one or more vehicles' states, which is necessary for the control and/or navigation system. Some of the most common cases in navigation systems are the use of barometers and/or sonars to improve the accuracy of the altitude provided by GNSS, or the use of small zenith cameras to determine small horizontal displacements [4]. Other specific navigation techniques such as visual odometry or visual Simultaneous Location And Mapping (SLAM) [5] are beginning to be used in indoor navigation.
Thus, other landing works focus on improving the accuracy of instrument systems such as [6][7][8] or even context information such as safe landing zones, as in work by Shah Alam, Md et al. shown in work [9]. Some commercial solutions focus on the use of beacons that indicate the landing region, as can be seen in the work of J. Janousek and P. Marcon [7] using a commercial infrared light beacon. This type of system includes an parameters. This idea attempts to improve the image resolution quickly to improve the position estimation. The difference with respect to other works is that simultaneously descending and adjusting the positioning relies on the variable waypoint altitude adjustment without accessing the controller, so that without changing the internal descent controller of the PX4 [26] or adding new external control laws, the strategy allows smoothing the descent when approaching the target. This strategy allows taking further steps in the final phase of the approximation, improving the final estimate, and ensuring the stability of the system provided by the manufacturer. For this, we propose a function to modulate the default behavior of the PX4 controller which seeks a stationary descent at a constant speed.
In addition, this paper models the error of a precision landing system using a monocular vision system, context information from the helipad pattern, and the internal navigation system of the PX4 flight controller.
The vision system error modeling allows a fine calibration of helipad localization. This correction, together with a stationary filtering of the estimates by sliding time window and variable adjustment of the descent height, allow to reduce the spin effect produced by the estimation error of a vision system when integrated with the L 1 navigation algorithm, without access to the internal controller parameters or relative states of the aircraft that may require re-programming of the on-board computer.
To sum up, this paper presents two main contributions in precision landing by visionbased global position: the continuous adjustment of the approach to descent trajectory, and the improvement of the position by vision through systematic error adjustment and filtering.
The proposal strategy is evaluated after including an error model correction as well as different sliding time window filters. The work is developed on a hyper-realistic software in the loop (SITL) [33] simulation system with the PX4 flight controller and the AirSim [34] simulator.
Finally, the results of the study show the estimation error analysis and filtering of the estimates with sliding time window filters, minimizing the inter-waypoint noise spin effect generated by noisy waypoint transitions and improving three quality metrics of landing time, trajectory landing length, and landing accuracy without additional control law, enabling the use of the aircraft's guidance system as an alternative for the deployment of precision landing technology. This paper is organized as follows: Section 2 shows the problem formulation of helipad spatial estimation by monocular computer vision. Section 3 describes the landing strategy proposal and the global estimation module. The correction module design, the analysis of the complete landing system, and a description of the test environment can be found in Section 4. Finally, the conclusions are presented in Section 5.

Problem Formulation
We consider the problem of a UAV landing on a certain landing pad using its internal autopilot waypoint guidance system and a monocular vision system with gimbal integrated in the UAV. Figure 1 shows the set of reference frames, where the superscripts {t, ph, c, z, g, b, n, e and g} correspond to the reference frames of the helipad (target), pinhole camera model, camera, gimbal socket for the camera, gimbal, body, North-East-Down (NED), Earth-Centered Earth-Fixed (ECEF), and global.
In this way, any given point p t t expressed in a flat pattern reference frame {t} and the homogeneous transformation n T t between the landing pad reference frame to the NED referent frame {n} can be expressed in the global reference frame {g} system p g t applying a set of transformations shown in Equations (1) and (2).
where superscript {j} over point p j i means the reference frame system, and the subscript i = {t, Re f } denotes the name of the point (target and reference). j T i means the homogeneous transformation between reference frame i and j. On the other hand, j F i refers to nonlinear transformations between reference frame i and j. p g Re f indicates the global position of the body (UAV) as a global reference point.
The set of reference frame systems involved in the transformations are shown in Figure 1  In this way, any given point expressed in a flat pattern reference frame { } and the homogeneous transformation . between the landing pad reference frame to the NED referent frame { } can be expressed in the global reference frame { } system applying a set of transformations shown in Equations (1) and (2).

⋅
(1) where superscript { } over point means the reference frame system, and the subscript = { , } denotes the name of the point (target and reference). . means the homogeneous transformation between reference frame and . On the other hand, . refers to nonlinear transformations between reference frame and .
indicates the global position of the body (UAV) as a global reference point.
The set of reference frame systems involved in the transformations are shown in Fig-Figure 1. General reference frame systems. Reference frames bottom-left to up: helipad (target), pinhole camera model (image plane), camera, gimbal socket, body, NED. Reference frames right up to down: NED, ECEF, Global.

Pattern (Helipad) Detection
We consider as the helipad a reference pattern defined by an Aruco pattern [35] with a certain number of bits, as part of a library B. As shown in the paper [20], the system identifies candidate square regions as Aruco markers, then encodes these regions and compares them to the pattern dictionary as desired.
The full process can be divided into the following steps: • Image conversion: Obtain an RGB image and transform it to grayscale. • Edge extraction: We understand as edge an intensity change boundary, some classical algorithms are Canny [36] and Sobel [37]. • Contour extraction: We understand a contour as a curve of points without gaps or jumps. Therefore, the objective is to identify if the edges found represent contours. An example of simple contour extraction can be given by a binarized image of an The transformation H is biunivocal and homogeneous, in other words, a point over a plane is a unique point over another plane and kH| k ∈ R and k = 0 is also the solution. This condition allows dividing the matrix H by the element h 33 , decreasing the dimension of terms to identify from 9 to 8. The correspondence between points (x i , y i ) ↔ x i , y i can be expressed in matrix form as b i = A i h, and their relationship is expressed as Equation (3) (more details in [38]). Knowing n pairs of points, the system of 2n equations and 8 unknowns is • Pattern library matching check: The binary code of the ROI is extracted, superimposing on the binarized and perspective-corrected image a grid of the same cell size as the searched one. Each grid cell receives a binary value according to if the corresponding color is black (zero) or white (one). The Hamming coding algorithm (ref) is applied to the extracted code to eliminate false negatives. This resulting code is compared with the selected pattern dictionary, filtering the regions identified as markers and belonging to the pattern dictionary from other regions. In addition, this step provides information about the marker id if the ROI belongs to the library.

Helipad Pose Estimation
For pose estimation, the Perspective-n-Point (PnP) problem [35] is formulated where the objective is to minimize the reprojection error Equation (7) of 3D points in the image plane {ph}. This problem is closely linked to a calibrated system, since it requires a camera model, pinhole, and a pattern that allows to relate identified features of an image with features of the pattern.
Given a point p c t ∈ R 3 belonging to the knowing pattern located in real-world 3D space and expressed in the camera reference frame {c}, it can be expressed in the image camera plane reference frame {ph} as p ph t ∈ R 2 . The relationship between the two reference frames is provided by the pinhole camera model in Equation (4).
where s is a scale factor and A intrinsic camera matrix [17]. The internal matrix A is composed of the focal distances f x , f y and the principal points c x , c y . The pinhole model can be improved with radial, tangential, or prism distortion corrections, adding n set of k i parameters to the model [39][40][41]. The set of internal parameters of the camera model can be expressed by the vector δ = f x , f y , c x , c y , k 1 , . . . , k n .
If the point is expressed in coordinates of the pattern reference frame p t t , there exists an extrinsic homogenous transformation c T t to relate the reference frame of the pattern to the camera reference frame Equation (5) is used.
where the transformation c T t = [ c R t |v c t ] is a rototranslation composed of the pattern's orientation c R t = R x (θ 1 )R y (θ 2 )R z (θ 3 ) to the camera and the pattern position vector to the camera v c t ∈ R 3 . Thus, the parameter vector to be identified to obtain the camera-pattern relationship is θ = (R, v) = (θ 1 , θ 2 , . . . , θ 6 ) ∈ R 6 .
Joining Equations (4) and (5) and adding distortion models, the camera model remains as a Function (6) that projects points p t t ∈ R 3 to p ph t ∈ R 2 points of the camera image plane.
Then, helipad pose estimation is the problem of minimizing the reprojection error Equation (7) of the observed helipad pattern features. One of the classic features to identify by computer vision are the corners. If the pattern is known, we know a priori the 3D position of these corners in the reference frame of the pattern.
where p t i ∈ and is a corner set of the pattern. O p t i ∈ R 2 is the corners obtained in the camera plane by a specific computer vision algorithm such as the Harris or Susan algorithm [22,42].
Furthermore, since all points p t i belong to a pattern plane, the z-component of all corners in the pattern frame will always be 0. This quality allows solving Equation (7) using specific methods such as the Infinitesimal Plane-Based Pose Estimation (IPPE) [43].
The estimation of internal camera parameters requires a learning phase modeled in Equation (7) as an optimization problem. In addition, identifying the six parameters to define the transformation c T t between the pattern calibration and the camera involves a similar process. Although both processes can be clustered as shown in Equation (7), the internal camera parameters δ will be constant for a particular vision system; however, the position of the pattern may change. For this reason, it is decoupled in two phases: on the one hand, a camera parameter learning (calibration) process, using a set of images of a known pattern to estimate internal camera parameters, and, on the other hand, the estimation of the helipad position for a certain image during flight.

Camera-Gimbal Frame
The camera is placed in a camera-gimbal socket, so it is necessary to include this referent frame {z}. As the camera-gimbal socket axis is equivalent with the general gimbal axes, but static, the z T c transformation is shown in Equation (8).
where R(x, θ) and R(z, ψ) represent θ and ψ rotations about the x and z axes of the camera reference system. As the axes of the gimbal and camera-gimbal socket are equivalent, the relationship between camera-gimbal socket and gimbal corresponds to the identity matrix g T z = I 4×4 .

Gimbal Body Frame
The gimbal's reference frame {g} to the UAV's body gravity center frame is defined as the composition of a roto-translation in Equation (9).
For this work, we consider the gimbal is static but located at the p b g position with R g (θ, φ, ψ) rotation to the body axes.

Body-NED Frame
The North-East-Down (NED) frame coordinates {n} to the UAV body gravity center {b}, n T b is equivalent to the body rotation at the angle defined by the yaw angle ψ to geographic north or azimuth, pitch attitude to horizon plane φ, and roll angle defined to gravity direction θ. These angles refer to the attitude and heading reference system (AHRS) frame of reference that groups magnetic, angular rate, and gravity (MARG) information. Generally, these systems usually include air data to provide altitude or wind speed information.

NED-ECEF-Global
The coordinate transformation between NED to the global reference frame {g} requires the use of the Earth-Centered Earth-Fixed (ECEF) reference system {e}, which allows us to apply the corresponding geodetic transformations to the terrestrial model and finally obtain the coordinates in global terms. In our case, we use a WGS84 (World Geodetic System 84) [44] datum.
The constant parameters of the WGS84 datum in Figure 2 refer to: r e semimajor axis (equatorial radius), r p semiminor axis (polar axis radius), ε first eccentricity and ε second eccentricity of the ellipsoid. It is important to differentiate the geocentric coordinates, referred to as the ECEF system, from the geodetic coordinates, referred to as the geodetic model (WGS84). This difference is provided by the geodetic model (datum) and is represented in the diagram on the right of Figure 2, where ϕ refers to geocentric latitude and ϕ refers to geodetic latitude.
Given a point p n t expressed in NED reference frame {n} of a local tangent plane (LTP) to a geodesic surface at a known point p g Re f = (λ, ϕ, h) T Re f , it can be expressed in ECEF coordinates {e} applying Equation (12). This equation corresponds to a translation in ECEF reference frame. However, to obtain p e Re f coordinates of our reference point in ECEF frame it is necessary to transform the global coordinates to ECEF applying Equation (14). The transformation between local coordinates and ECEF is given by the transformation Equation (13). In this work, we consider p g Re f = p g U AV . Given a point expressed in NED reference frame { } of a local tangent plane (LTP) to a geodesic surface at a known point = ( , , ℎ) , it can be expressed in ECEF coordinates { } applying Equation (12). This equation corresponds to a translation in ECEF reference frame. However, to obtain coordinates of our reference point in ECEF frame it is necessary to transform the global coordinates to ECEF applying Equation (14). The transformation between local coordinates and ECEF is given by the transformation Equation (13). In this work, we consider = . On the other hand, a given point expressed in ECEF can be expressed in global coordinates applying the transformation Equation (11).
where ( , , ℎ) means longitude, latitude, and altitude in WGS84 datum. ( , , ) are the coordinates in the ECEF reference frame. The transformation Equation (11) corresponds to Jijie Zhu's algorithm [45] analyzed and compared in [46]. On the other hand, a given point p e t expressed in ECEF can be expressed in global coordinates p g t applying the transformation Equation (11).
where (λ, ϕ, h) T means longitude, latitude, and altitude in WGS84 datum. (x e t , y e t , z e t ) T are the coordinates in the ECEF reference frame. The transformation Equation (11) corresponds to Jijie Zhu's algorithm [45] analyzed and compared in [46].

Proposal
In this section, first the landing strategy is described, then the method to determine the helipad's global position is detailed, and finally the error analysis of the helipad's position estimation is given.

Landing Strategy
The landing strategy is responsible for telling the UAV navigation system the position to which it must go and the attitude it must have to align with the target (Algorithm 1). The position of the target is static, but the attitude and altitude to helipad vary overtime when the UAV attempts to land. To align the drone to the marker, it is necessary to determine the azimuth of the marker ψ n t . For this, we use the azimuth of the drone ψ n U AV and the orientation of the marker to the drone ψ b t . Figure 3 shows how the helipad azimuth can be obtained graphically by adding to the drone azimuth the orientation of the aircraft to the landing pad in Equation (16). When both systems are aligned, the marker azimuth will be equal to the drone azimuth.

Altitude Setpoint Strategy
In order to change the default controller descent behavior, it is possib behavior of the system in the transient state, i.e., before reaching the maxim

Altitude Setpoint Strategy
In order to change the default controller descent behavior, it is possible to use the behavior of the system in the transient state, i.e., before reaching the maximum velocity of stationary descent. Thus, if the new desired height is reached without having to reach the maximum descent speed, the behavior will be smooth, and if the destination point is far enough away, the controller saturates and descends with maximum constant speed behavior, without exceeding the internal controller parameters.
To define step points that allow a linear descent at constant speed α in an iterative loop, the new step point will correspond to the current height minus a certain parameter α.
Considering this process is iterative (discrete) with a sample time of ∆t, the previous equation can be expressed as follows: where t subscript means instant time. Solving the α term, it is verified that alpha corresponds to a speed term.
In our case, the aim is to design the h set (h U AV ) function such that the aircraft approaches with a smooth behavior to h 0 and at that point lands automatically with internal autopilot.
To perform this, we propose Function (23).
where β 1 is the weight of the exponential function and β 0 < 1, which allows slightly shifting the value of h 0 and to be able to switch to automatic landing mode. The β 1 = 0.1 value is set heuristically, while β 0 = 0.5 is set to shift 50% less than the switching height h 0 . The system must consider the relative flight altitude h U AV and the height of the UAV relative to the landing pad h b t . Figure 4 shows an approximate representation of the altitude-set function behavior Equation (20) vs. constant decreasing Equation (17).   Figure 4 shows the approximate descent behavior of our proposal versus a constant speed descent. In the final phase of the approximation, the descent becomes smoother than in the linear behavior. The Figure 4 behavior should be taken as an illustrative example of the desired behavior, not as a realistic simulation. The final behavior can be seen in  Figure 4 shows the approximate descent behavior of our proposal versus a constant speed descent. In the final phase of the approximation, the descent becomes smoother than in the linear behavior. The Figure 4 behavior should be taken as an illustrative example of the desired behavior, not as a realistic simulation. The final behavior can be seen in experimentation.

Filter
The states to be filtered are the global position of the helipad p g t = (λ t , ϕ t ) and its orientation to north or azimuth ψ n t . All these variables are static, since the landing pad is static; therefore, the filter model does not need to provide information for each new measurement, rather we need to know their stationary statistical values. For this we propose to generate a data buffer with memory. The size of the buffer defines the size of the filtering window L. The initialization saves L new measurements and then finds the mean or median of the buffered data. Finally, the buffer is reset.
To propagate the information over time, the sliding window does not overlap with previous values, but the value filtered at the previous instant is included as the first measurement in the clean buffer.
As for the filter memory, if the new values change substantially it will vanish in the long term, since the weight given to the past values is 1 L versus L−1 L for each new data, so the window size can be critical for cases where the target is moving. Finally, the size of the buffer/window L is linked to time thanks to the 1 Hz system sampling time to provide new measurements.

Helipad Global Position Estimation
The helipad global position estimation system is responsible for integrating the vision system, the heliport context information, the gimbal, and the UAV navigation states, to provide the landing strategy with the helipad global position. In addition, this includes a spatial correction system for the NED frame, which is the objective of study of this work. Figure 5 shows the diagram of the estimation system which is formulated in Equation (2). The system works as follows:

•
Aircraft global position p g U AV and attitude (θ, φ, ψ) b U AV is requested by the PX4 flight controller via MAVLink protocol [24] supported by the MAVSDK API [25].

•
The gimbal position p b g (x, y, z) and attitude (θ, φ, ψ) b g is requested by the AirSim simulation environment via UDP protocol described in Section 4.1. This information composes the b T z transform. • The vision system receives I image of W × H size and 3 RGB channels. The image is received via UDP protocol from the simulation system. In addition, the vision system has as input the context information from the helipad, the library (Lib) of the marker, the marker's identification number (Id ∈ Lib), and the real marker's size (MS) in meters. The library is characterized by the number of horizontal and vertical bits (squares) that form the geometry of the marker and the number of elements that make up the library. The vision system output provides a Boolean variable a ∈ B, that indicates if the landing pad has been detected or not. In addition, it provides the position of the landing pad to the camera p c t and the attitude (θ, φ, ψ) c t .

•
The camera pose estimation (p t c and t T c ) is gated by the PnP method integrated in the OpenCV Aruco library [47] from a previously pre-calibrated camera (Test environment).

•
The aircraft, gimbal, and landing pad position are combined in the set of n T b , b T z , and z T c transformations to obtain the positioning p n t and attitude (θ, φ, ψ) n t of the landing pad in NED frame.

•
The correction module provides the p n t positioning and attitude (θ, φ, ψ) n t , tuned in NED coordinates. • Finally, the target position in NED frame p n t , together with the drone global position p g U AV and ellipsoid WGS84 approximation, are used to obtain the helipad global position in Equation (11). ing pad in NED frame.

•
The correction module provides the ′ positioning and attitude ( , , )′ , tuned in NED coordinates. • Finally, the target position in NED frame ′ , together with the drone global position and ellipsoid WGS84 approximation, are used to obtain the helipad global position in Equation (11).

Landing System Analysis
The aim of this section is to evaluate the proposed estimation system and to identify the necessary corrections to be incorporated in the "correction" module of Figure 5. To

Landing System Analysis
The aim of this section is to evaluate the proposed estimation system and to identify the necessary corrections to be incorporated in the "correction" module of Figure 5. To achieve this, first the test environment and the necessary parameters are detailed in the subsection Test environment. Next, the system estimation error is modeled to provide the landing system a correction module. The quality of the correction is evaluated using the root mean square error (RMSE) together with the variation in the data distribution in terms of data distribution structure, mean, and standard deviation.
Finally, a full landing system and classical linear decreasing descent are compared and evaluated with four quality metrics, which quantify the trajectory length, the time to land, and the accuracy of landing on the helipad.

Test Environment
In this work, we use a hyper-realistic test environment based on Software in The Loop (SITL). SITL systems are simulation architectures where virtual world environments interact to simulate object, vehicle, and sensor together with external systems such as a flight controller or ground station, among others. These environments are powerful testing tools for earlier phases of system integration, as they allow realistic results to be obtained without potentially dangerous and expensive risks.
In our case, AirSim [48] is used as world environment and PX4 flight controller configured as a quadcopter. The simulated physical model corresponds to the Iris quadcopter and the set of sensors, and their specifications are detailed in Table 1. The models of the simulated sensors can be found detailed in [34].  Figure 6 shows an SITL communication diagram between the main system modules in SITL. The GCS module refers to the ground control station, in our case QGround control [50]. GCS is used to help to download the .log files generated in the test missions. The vision-based estimation system requires knowledge of the internal camera parameters { , }. These parameters are obtained by standard calibration [39] using a chess pattern with nine rows, six columns and 20 cm sides of the squares. This pattern is integrated into the AirSim environment as a texture over a rectangular prism with 1.8 × 1.2 × 1.8 [ ] sides. To capture images, we implemented a system that automatically captures images while performing a spiral upward flight over the reference pattern. This allows obtaining a large set of images with the pattern from different positions. Figure 7 shows the image of the calibration pattern in the AirSim reference frame, random image of the image registration process used for calibration, and a sample of the reprojection error. The vision-based estimation system requires knowledge of the internal camera parameters {A, k i }. These parameters are obtained by standard calibration [39] using a chess pattern with nine rows, six columns and 20 cm sides of the squares. This pattern is integrated into the AirSim environment as a texture over a rectangular prism with 1.8 × 1.2 × 1. 8 [m] sides. To capture images, we implemented a system that automatically captures images while performing a spiral upward flight over the reference pattern. This allows obtaining a large set of images with the pattern from different positions. Figure 7 shows the image of the calibration pattern in the AirSim reference frame, random image of the image registration process used for calibration, and a sample of the reprojection error. grated into the AirSim environment as a texture over a rectangular prism with 1.8 × 1.2 × 1.8 [ ] sides. To capture images, we implemented a system that automatically captures images while performing a spiral upward flight over the reference pattern. This allows obtaining a large set of images with the pattern from different positions. Figure 7 shows the image of the calibration pattern in the AirSim reference frame, random image of the image registration process used for calibration, and a sample of the reprojection error. Finally, the internal camera parameters are shown in Table 2. The context information used for the experimentation is:  Finally, the internal camera parameters are shown in Table 2. The context information used for the experimentation is: Lib = 5 × 5 × 1000, Id = 68, and MS = 1 [m]. The landing system was developed in Python 3.6 with the AirSim [34] and MAVSDK [25] APIs. The experiments and the SITL environment were developed on a Windows Server 2019, 64 bits, hosted in AMD Ryzen 9 3900X 12-Core Processor CPU, 3.79 GHz with 64 GB RAM and 2×1TB SSD + 2×HDD 1.5TB of internal memory, graphic card Nvidia GeForce RTX 2060.

NED Error Modeling
To evaluate the estimation error, we propose to analyze the estimation data provided by the vision system over twenty static flights located at seventeen different positions at the same relative altitude above the ground, 10 meters.
The selected positions correspond to five different headings centered, 45 • {northeast, southeast, southwest, and northwest} and four different distances {2, 3, 4, 5} √ 2 to the takeoff origin where the helipad is located. In each position is recorded a total of 1000 p t U AV samples. Looking at Figure  We consider the aircraft control system is asymptotically stable so that in steady state its position converges to the reference one. Thus, we consider as ground truth the reference positions for the steady flight.
The error position for each of the components is given by Equation (21).
where e(x)  We consider the aircraft control system is asymptotically stable so that in steady state its position converges to the reference one. Thus, we consider as ground truth the reference positions for the steady flight.
The error position for each of the components is given by Equation (21). Looking at the errors (Figure 9), the error distribution increases with increasing distance from the north-east plane origin (0, 0). This means the error position depends on the position in the NE plane.  Regarding altitude error, Figure 10a,b show for each twenty register positions a distribution with four "modes". In Figure 10a these modes show as four scatter clusters and in Figure 10b as four peaks in each twenty distributions. In addition, the mean and median of the total error distribution, Figure 10, are displaced from the origin, showing a bias in altitude.
The different colors in Figure 10b show each of the twenty records, all of them showing four modes and centered on the same error terms. In this work, we focus on bias correction of mean and median; however, modeling the error in altitude is outside the scope of this paper. Regarding altitude error, Figure 10a,b show for each twenty register positions a distribution with four "modes". In Figure 10a these modes show as four scatter clusters and in Figure 10b as four peaks in each twenty distributions. In addition, the mean and median of the total error distribution, Figure 10, are displaced from the origin, showing a bias in altitude.
The different colors in Figure 10b show each of the twenty records, all of them showing four modes and centered on the same error terms. In this work, we focus on bias correction of mean and median; however, modeling the error in altitude is outside the scope of this paper.
The altitude error distribution may be a consequence of the internal discretizing of the simulator in the image render, so that the vision system, when segmenting the ROI of the helipad, extracts its contour with a size variation. This would be explainable according to the pinhole model of Equation (4) and PnP Formulation (7), since its scale factor is constant, but the size of the ROI and corner positions would change.
Regarding altitude error, Figure 10a,b show for each twenty register positions a distribution with four "modes". In Figure 10a these modes show as four scatter clusters and in Figure 10b as four peaks in each twenty distributions. In addition, the mean and median of the total error distribution, Figure 10, are displaced from the origin, showing a bias in altitude.
The different colors in Figure 10b show each of the twenty records, all of them showing four modes and centered on the same error terms. In this work, we focus on bias correction of mean and median; however, modeling the error in altitude is outside the scope of this paper. The altitude error distribution may be a consequence of the internal discretizing of the simulator in the image render, so that the vision system, when segmenting the ROI of the helipad, extracts its contour with a size variation. This would be explainable according

Polar Space Error Analysis
The visual results in Figures 8 and 9 show an apparent angular and radial bias of the helipad global position estimation system. We change the cartesian space to the cylindrical space defined by Equation (22), where the terms E, N, D are the coordinates in the NED referent frame.
When plotting data in the new space in Figure 11, it can be seen how the data set is apparently clustered around a constant bias in the angle and radial error. to the pinhole model of Equation (4) and PnP Formulation (7), since its scale factor is constant, but the size of the ROI and corner positions would change.

Polar Space Error Analysis
The visual results in Figures 8 and 9 show an apparent angular and radial bias of the helipad global position estimation system. We change the cartesian space to the cylindrical space defined by Equation (22), where the terms , , are the coordinates in the NED referent frame.
When plotting data in the new space in Figure 11, it can be seen how the data set is apparently clustered around a constant bias in the angle and radial error. However, when showing the distance error (radial) behavior versus distance, Figure  12 shows a high linear correlation between distance and radial error. The angular error is also tested for linear dependence on distance, but the correlation does not exceed 35% of variance score in Equation (23), so it has been discarded. However, when showing the distance error (radial) behavior versus distance, Figure 12 shows a high linear correlation between distance and radial error. The angular error is also tested for linear dependence on distance, but the correlation does not exceed 35% of variance score r 2 in Equation (23), so it has been discarded.
where Var(x −x) indicates the variance of the error betweenx model estimation and x data.
Where ( − ) indicates the variance of the error between model estimation and data. Figure 12 shows with "+" the radial error centers of each of the measurement positions and the linear model fitted (blue line) by least squares to these points. This model has a slope = 2.4923 × 10 and independent term = 2.778497 × 10 [ ]. The linear model obtains around 94% of the variance score in Equation (23). Given the results in Table 3 and Figure 12, the error bias in polar space can be tuned by the model of Equation (24), where the apostrophe over coordinates means corrected coordinate, subscript zero start value, and the bias of coordinate .

Error Correction in NED Space
Given , , and coordinates of a point in the NED frame and knowing the correction in a cylindrical space in Equation (24), the objective is to return to the NED space. For this purpose, we apply the transformation Equation (25).

=.
( ) where its terms , are taken as shown in Equation (26):  Figure 12 shows with "+" the radial error centers of each of the measurement positions and the linear model fitted (blue line) by least squares to these points. This model has a slope α r = 2.4923 × 10 −2 and independent term β r = 2.778497 × 10 −2 [m]. The linear model obtains around 94% of the variance score in Equation (23).
Given the results in Table 3 and Figure 12, the error bias in polar space can be tuned by the model of Equation (24), where the apostrophe over coordinates means corrected coordinate, subscript zero start value, and β i the bias of coordinate i.

Error Correction in NED Space
Given N, E, and D coordinates of a point p n i in the NED frame and knowing the correction in a cylindrical space in Equation (24), the objective is to return to the NED space. For this purpose, we apply the transformation Equation (25).
where its termsr, θ are taken as shown in Equation (26): where the hat over N, E, and D in Equation (26) means the NED coordinate with cylindrical corrections. β {θ,r,D} are the biases of the radial, angular, and altitude terms, respectively. v β θ components and θ mean new position and new angular position after β θ rotation correction. Figure 13 shows the error distributions of the raw position estimation and error distribution after applying the correction Equation (26) where the hat over , , and in Equation (26) means the NED coordinate with cylindrical corrections. { , , } are the biases of the radial, angular, and altitude terms, respectively.
components and mean new position and new angular position after rotation correction. Figure 13 shows the error distributions of the raw position estimation and error distribution after applying the correction Equation (26) with the parameters of Table 1.

North
East Down For each drawing in Figure 13, the distributions of the real data are shown in blue and with a blue line their error distribution function [51]. The orange line shows a Gaussian distribution equivalent to the real data in Equation (27). For the NED coordinate, three different figures are represented: raw data Figure 13a-c, data after cylindrical correction using the mean of the raw data Figure 13d-f, and data after cylindrical correction For each drawing in Figure 13, the distributions of the real data are shown in blue and with a blue line their error distribution function [51]. The orange line shows a Gaussian distribution equivalent to the real data in Equation (27). For the NED coordinate, three different figures are represented: raw data Figure 13a-c, data after cylindrical correction using the mean of the raw data Figure 13d-f, and data after cylindrical correction according to the median of the raw data Figure 13g-i. In addition, the mean value and the standard deviation of each case represented by µ and σ, respectively, are indicated on each graph in their legend. Figure 13 shows how correction Equation (25) modifies the structure of the error distribution for the cases of N and E coordinates, when comparing the blue with the orange lines. It can be seen in Figure 13a,b how the initial distribution is close to a Gaussian distribution Equation (27) and Figure 13d-h.
The effect on the Gaussian approximations in mean and standard deviation represented in Figure 13 is quantified in Table 4.  Table 4 shows the effect of mean and standard deviation on the data when applying the correction Equation (25) using the mean and median bias value indicated in Table 3.
The standard deviation rows of Table 4 decrease one order of magnitude in all coordinates when applying the correction Equation (25). The same effect can be seen in Table 5 when the RMSE Equation (28) of each NED coordinate is calculated. This metric can be considered as an indicator of accuracy.

Landing Evaluation
To evaluate the landing system, we propose to test twenty landing missions from the same position at (10, 10, −20) NED meters to the helipad. The twenty flights are divided into four groups corresponding to using the raw landing system without correction (without) Equation (2), applying the bias correction Equation (25) (Bias), with bias correction and a mean filter (Section 3.1.3) with a sliding time window (Mean&Bias) and a median filter together with the bias correction (Median&Bias). For each mission, we use as quality metrics the landing trajectory distance Equation (29), time to land Equation (30), and landing accuracy. In addition, the results are compared with classical linear descent setpoints with α = 0.7 [m/s].
where the k term means temporal step starting in t start instant and ending in t land moment. p n k represents the global position at k instant in the NED reference frame. For the flights' analysis, we use the information obtained from the logs recorded by the flight controller in each flight and unloaded with the ground control station (GCS). In particular, we focus on the global positioning of the UAV. This positioning is given by the EKF2 fusion system [52] integrated in PX4 and with the specific sensor parameters indicated in the SITL [33] (Section 4.1). To ensure that we evaluate exclusively the landing phase of the UAV, we study the trajectories from the instant t start where the height gradient is detected to be negative, and the altitude is less than 99.8% of the altitude desired. The final instant t land is obtained when the helipad is reached by the same method as t start .
Finally, third quality metric, landing accuracy, is obtained with the RMSE of the last ten position samples of the NED coordinates (without altitude). In this way, we ensure that we are on the ground with the same value plus a precision error. For this, we take as ground truth the control position of the marker. Figure 14 shows the behavior of the aircraft when activating the landing system with the four modes corresponding to the landing system without correction (Figure 14a), landing with bias correction (Figure 14b), landing with bias correction and mean filter (Figure 14c,d), and landing with bias correction and median filter. The five trajectories in each of the figures show a different flight, using the corresponding landing mode in each case. The total number of flights is five for each landing mode, i.e., twenty flights.
phase of the UAV, we study the trajectories from the instant where the height gradient is detected to be negative, and the altitude is less than 99.8% of the altitude desired. The final instant is obtained when the helipad is reached by the same method as .
Finally, third quality metric, landing accuracy, is obtained with the RMSE of the last ten position samples of the NED coordinates (without altitude). In this way, we ensure that we are on the ground with the same value plus a precision error. For this, we take as ground truth the control position of the marker. Figure 14 shows the behavior of the aircraft when activating the landing system with the four modes corresponding to the landing system without correction (Figure 14a), landing with bias correction (Figure 14b), landing with bias correction and mean filter ( Figure  14c,d), and landing with bias correction and median filter. The five trajectories in each of the figures show a different flight, using the corresponding landing mode in each case. The total number of flights is five for each landing mode, i.e., twenty flights.  Figure 15 illustrates the temporal behavior of each analysis group, showing the three global position components: latitude, longitude, and relative altitude. In addition, our proposal is comparing with linear descent, paying special attention to altitude evolution Figure 15e,f. In this case, estimated altitude and setpoints are shown for the exponential proposal and linear descending.
It can be seen in all cases in Figure 15 how the position of the landing pad, defined by a blue dashed horizontal line, is obtained in the stationary state. The effect of the error in precision is shown with an oscillating behavior (blue dotted line). However, when bias Equation (26) and filtering (Section 3.1.3) corrections are applied, the effect is damped.   It can be seen in all cases in Figure 15 how the position of the landing pad, defined by a blue dashed horizontal line, is obtained in the stationary state. The effect of the error in precision is shown with an oscillating behavior (blue dotted line). However, when bias Equation (26) and filtering (Section 3.1.3) corrections are applied, the effect is damped. The linear decrease approach converges in latitude and longitude (Figure 15b,d), but the system finds it difficult to dampen the spin effect. In the case of an exponential decrease, the inter-waypoint spin effect is considerably less than linear (Figure 15a-d). In both cases, the linear and exponential altitude decreasing approaches ( Figure 15) show when error correction is applied and filtering estimation inter-waypoint noise spin effect is smoothed. Figure 15e,f show a small break at the end of the trajectory that exceeds the reference and picks it up again. To identify the instant to obtain the landing surface, we keep the first term that satisfies the gradient and proximity conditions explained above.
This effect may be a consequence of a decoupling of the fusion system in the estimation of the height, for example for giving more weight to the estimation than to the measurements in the EKF2 filter. Therefore, when identifying the instant of reaching the pad, we are left with the first term that meets the conditions of gradient and proximity explained above.
The total of twenty flights with exponential decreasing approach and the other twenty flights with linear decreasing approaches are summarized in Tables 6 and 7. These It can be seen in all cases in Figure 15 how the position of the landing pad, defined by a blue dashed horizontal line, is obtained in the stationary state. The effect of the error in precision is shown with an oscillating behavior (blue dotted line). However, when bias Equation (26) and filtering (Section 3.1.3) corrections are applied, the effect is damped.
The linear decrease approach converges in latitude and longitude (Figure 15b,d), but the system finds it difficult to dampen the spin effect. In the case of an exponential decrease, the inter-waypoint spin effect is considerably less than linear (Figure 15a-d). In both cases, the linear and exponential altitude decreasing approaches ( Figure 15) show when error correction is applied and filtering estimation inter-waypoint noise spin effect is smoothed. Figure 15e,f show a small break at the end of the trajectory that exceeds the reference and picks it up again. To identify the instant to obtain the landing surface, we keep the first term that satisfies the gradient and proximity conditions explained above.
This effect may be a consequence of a decoupling of the fusion system in the estimation of the height, for example for giving more weight to the estimation than to the measurements in the EKF2 filter. Therefore, when identifying the instant of reaching the pad, we are left with the first term that meets the conditions of gradient and proximity explained above.
The total of twenty flights with exponential decreasing approach and the other twenty flights with linear decreasing approaches are summarized in Tables 6 and 7. These tables are built with the mean and median values of all flights, so the quality metrics means the mean and median values of all flights.  Tables 6 and 7 show the results of the exponential and linear decreasing approaches grouped for easy comparison. The previous tables show for all cases that the exponential descent proposal improves the results of the linear descent, emphasizing that in the bestcase scenario of the linear approach (Mena&Bias), the trajectory distance is reduced by 32%, the time to land by 61%, and the RMSE of the precision landing by 12%.
In both tables, the mode that provides the minimum mean landing trajectory distance is the bias correction together with the median filter. The minimum mean time to land is provided by the bias correction together with the median filter and the minimum mean RMSE is again provided by the bias correction together with the mean filter. Finally, the similarity between the mean and median values in Tables 6 and 7 are a good indicator of normal distribution.

Conclusions
Through the study of the global position estimation system error of Section 4, an angular and radial bias was identified. In addition, it was shown how the error distribution increases its degree of dispersion with the distance to the origin. This position error was initially modeled in a cylindrical space in Equation (24) and transferred to the NED reference space under the transformation Equations (25) and (26). The corrections over cylindrical space produced a structural transformation of the position error distribution, approximating its distributions to the Gaussian error.
On the other hand, it was verified how using a system aimed at smoothing trajectories between waypoints can produce a spin effect if the new waypoints are updated with a frequency such that the UAV cannot obtain the previous target and these new waypoints correspond to the same position, but with high uncertainty. Therefore, the path planning system with path smoothing between waypoints can work as an error amplifier performing a circular trajectory.
Finally, we conclude that the combination of an exponential altitude decrease, together with the correction of systematic estimation error and a sliding time window filtering, improves all three quality metrics proposed and reduces the effect of the inter-waypoint noise spin effect. These results facilitate the development of new applications that require a lightweight but robust precision landing strategy.