Visual Odometry in GPS Denied Zones for Fixed-Wing UAV with Reduced Accumulative Error Based on Satellite Imaginary

1. Introduction

Unmanned Aerial Vehicles (UAV) were initially developed for military purposes [1]. The goal was to create realistic aircraft capable of reconnaissance missions for the U.S. Air Force. Now, drones have evolved into safe and cost-effective systems equipped with multiple sensors, finding extensive use in civilian, commercial, and scientific domains. They are utilized for applications such as aerial photography, infrastructure inspection, and search and rescue.

There are two types of motorized UAVs:

• Multirotor drones. These UAV are characterized by their design with multiple rotor blades. These drones are similar in concept to traditional helicopters but are typically smaller. The main distinguishing ability of rotary-wing drones is their ability to hover and take off and land vertically (VTOL).
• Fixed-wing. Similar to traditional airplanes. Unlike rotary-wing drones, they achieve lift through the aerodynamic forces generated by their wings as they move through the air.

Multirotor drones have some advantages over fixed-wing drones, such as the ability to hover, take off and land vertically, and change direction quickly. However, they also have some disadvantages, such as lower speed, shorter range, and higher energy consumption. This work focuses on fixed wing UAVs for which there are multiple classification tables based on size, weight, range, endurance, etc. The Table 1 shows a classification from [2]; The UAV used for image capture is categorized as small. Because of its range and operating altitude, this kind of UAV is suitable for reconnaissance, and monitoring.

To ensure safe takeoff, flight, and landing, these vehicles are equipped with multiple onboard sensors, such as GPS (Global Positioning System), environmental sensors (barometers and humidity sensors), and Inertial Navigation Systems (INS), which combine GPS data with INS data for accurate and reliable position information. In addition, depending on the task, imaging sensors (RGB or hyperspectral cameras), Light Detection and Ranging (LIDAR), air quality sensors, etc., may also be included.

The UAVs are piloted by a human on a Ground Control Station (GCS). The UAV is equipped with various actuators to control the surfaces of the aircraft (such as ailerons, elevators, and rudder), and an autopilot system which is responsible for stabilizing the aircraft, following pre-programmed flight paths, and responding to commands from the GCS.

The UAV communicates with the GCS. This communication is crucial for transmitting commands, receiving flight data, and sending back telemetry. The UAV continuously sends telemetry data to the GCS and it includes:

• Altitude, speed, and heading information.
• Battery status and remaining flight time.
• Sensor data such as temperature, humidity, or camera feed.
• GPS coordinates for position tracking.

From the GCS, it is possible to choose between different flight modes, such as:

• Manual
• Altitude hold
• Loiter mode
• Autonomous (pre-programmed missions)

The overall structure of this paper is as follows: The problem of GPS-denied is stated in Section 2. Section 3 presents a review of the state of the art. In Section 4 we describe the algorithm and present the results. Finally, we discuss how your findings in Section 6.

2. Problem Statement

GPS signal is vital in order to achieve autonomy. Unfortunately, GPS signals can be intentionally blocked or disrupted in certain zones, including secure facilities, prisons, or sensitive government installations. There are multiple ways in which the GPS signal is intentionally restricted. For example, GPS jamming involves broadcasting radio signals on the same frequencies used by GPS satellites, overpowering the weaker GPS signals.

GPS spoofing is another way to restrict GPS. Here, a fake GPS signal that mimics legitimate signals is broadcast. GPS receivers are tricked into calculating incorrect position or velocity. Detecting GPS spoofing can be challenging, but several methods and technologies can help identify if a GPS signal is being spoofed [3,4,5,6,7,8,9,10].

Regardless of how the GPS signal is compromised, the mission may fail, since no other onboard sensor, like the Inertial Measurement Unit (IMU), is able to provide reliable position, velocity and altitude of the UAV since the state estimation drifts in time and becomes unusable after few seconds [11,12]. As a result, the UAV will take action depending on the flight mode or if the human pilot has visual contact with the vehicle, for instance, if there are no visual contact, most of the UAV have an automated routine to return to home, preventing it from completing its task. It is needed a solution where the UAV is able to keep flying or complete a mission even if the GPS signal is corrupted.

3. Related Work

There are multiple solutions for the GPS-denied problem, divided mainly into two approaches. The first approach is based on information from other agents, such as another vehicle [13,14,15], Unattended Ground Sensors (UGS) [16], or other types of packet transmission devices that serve as reference points for the UAV [17,18,19,20]. However, these techniques require the prior installation of these agents along the flight path, which can be unfeasible, especially if the agents are antennas that are difficult to transport or if passage to certain areas of the trajectory is restricted.

Not relying on other agents adds another level of autonomy to the system. Other approaches to addressing the GPS-denied problem are based on Visual Information. Some approaches use Simultaneous Localization And Mapping (SLAM), which is a technique used to locate the UAV in unknown environments [21,22,23]. However, depending on the altitude and the range, achieving visual SLAM techniques may not be feasible. Most of these works, however, are only reliable when the problem of localization is at low altitude and low range [11].

The main problem of using the previous approach for high-altitude UAVs (> 1000 m) is the estimation of the scale. That is, the transformation between pixel displacements and meters. Moreover, the nonlinear transformation between meters and latitude/longitude coordinates necessitates the prior setting of the Universal Transverse Mercator (UTM) zone. This scale could be computed based on IMU information and any other altitude sensor, such as LIDAR or barometer. However, the limitation of the LIDAR lies in its range, while the limitation of the barometer lies in its resolution.

When the UAV operates at high altitudes, LIDAR ceases to function, and we must rely on the barometer to calculate sea-level altitude. Then, with the help of an elevation map, we can calculate the difference between altitude above sea level and terrain elevation, thus determining a scale. However, limitations in barometer resolution result in errors when converting pixel displacement to meters, leading subsequent readings to start from an incorrect latitude/longitude coordinate, thereby resulting in a misreading of the elevation map. Without assistance or prior knowledge (landmarks) of the area we are flying in, errors will accumulate over time.

This paper is based on visual information and landmarks. In [11], the authors start from a previously created map. The helicopter flies over the area of interest to collect image data and construct a feature map for localization. The area is divided into a set of discrete regions (nodes), which are related to each other. The problem with this approach is that the UAV needs to previously fly over the area, and it would not be possible if the GPS signal is corrupted since the beginning of the flight.

Our work is more similar to [12]. Here, the authors start with geo-referenced satellite or aerial images, and then they search for landmarks. However, one disadvantage is that the range is limited, ensuring that the scale between pixels and meters remains constant. Even if we know the scale, it is impossible to guarantee that the elevation of the terrain remains constant in real reconnaissance or surveillance applications.

The contribution of this work is an algorithm that reduces the accumulated error in the estimation of the GPS signal based on aerial and satellite images for a UAV flying at over 1000 meters altitude above the terrain and a range of over 15 kilometers.

4. Methodology

The algorithm consists of two phases:

• Visual Odometry: At the moment of GPS loss, the algorithm estimates a latitude/longitude coordinate based only on monocular vision. The accumulated error in this phase will depend on the precision with which the scale is calculated.
• Reduction of accumulated error: We search for Correspondences between the UAV image and the georeferenced map to correct the estimation errors from the previous phase.

Assumption A1. 

The aircraft is equipped with GPS signal from takeoff until some time after reaching the operating altitude.

Assumption A2. 

We know the terrain elevation map

4.1. Visual Odometry

The experiments where performed on a small fixed-wing UAV flying with an altitude greater than 1000 meters. 700 Photos were taken on a path of approximately 18 km. Figure 1 shows the real path obtained from the GPS.

It is assumed the UAV enters into a GPS denied zone at the start of the path in Figure 1. When GPS is lost, the plane takes the first photo at (latitude, longitude) = (20.64226,-103.8181417). The path ends at (20.641845,-103.774705) and Figure 2 shows both images.

Figure 3 describes this stage of the algorithm. It starts when the UAV loses GPS signal, it begins to calculate the location estimation based on two consecutive images. First, features are extracted from both images using Oriented Fast and Rotated Brief (ORB).

Although ORB is one of the worst descriptors in terms of scale invariance, we can assume that the scale does not change substantially between both images, as the UAV will not change its altitude drastically during its operation. On the contrary, it is one of the algorithms that finds the largest number of keypoints and also one of the fastest. An analysis among different descriptor algorithms (ORB, BRISK, KAZE, AKAZE, SURF, and SIFT) can be found in [24]. The issue with detecting so many keypoints is the time it will take to establish correspondences between the images. To avoid processing all of them, we will only use matches with the shortest distance. This distance is a measure of the feature point similarity and we will denote with $\rho$.

The first step is to consider the first two consecutive images and extract ORB features. We filter these matches to keep only the best ones. The best matches will be those whose similarity distance $\rho$ between features is less than a threshold. This threshold will depend on the images and the descriptor algorithm used. The heuristic used for this process is as follows: A test is conducted by plotting the best 200 matches, and the result is shown in the Figure 4.

Those matches were sorted based on similarity $\rho$. This means that after 200 matches we have no matching error. The distance of similarity between features in the last match (the worst) in this experiment is 691. based on this information, we set a threshold of 300, trying to ensure that we only have good matches.

Once we know the correspondence between both point clouds, we compute the centroid of them

$$C_i=\left(\frac{1}{n}\sum _{j=1}^n{x}_{ij},\frac{1}{n}\sum _{j=1}^n{y}_{ij}\right)$$

(1)

where $C_i$ is the centroid of the 2D points in the image i, x and y are the pixel coordinates of the feature j in image i and n is the number of good matches.

The idea is to estimate the UAV displacement from the displacement of one centroid to the next one. It is important to note that we have access to the angle of rotation since this is information obtained from the IMU and it is available even the GPS has been corrupted. Then, it reduces to calculating the UAV displacement in meters from the displacement of the centroids in pixel units.

There is an issue of ambiguity when trying to convert from pixels to meters. It may occur that a large displacement between the centroids is due to either moving very fast or having the keypoints very close to the UAV. Since altitude (estimated from the barometer) is another piece of information we have access to even without GPS, we use the barometer to heuristically scale s the altitude and use this value to convert from pixel units to meters. The scale s is calculated with

$$s=h\delta$$

(2)

where s is the scale to transform from pixel to meters. $\delta$ is a constant heuristically selected and h represents the distance from the UAV to the ground, this is

$$h=a-e$$

(3)

where a is the UAV altitude and e the terrain elevation, taken from the digital elevation model (DEM) which is displayed in Figure 6.

Now it is needed to convert the last GPS signal from latitude/longitude to UTM coordinates. This is done using the open-source Geospatial Data Abstraction Library (GDAL). We start from the last coordinates obtained from GPS (before it has been denied) and convert them to UTM coordinates (Figure 7). In UTM system, the Earth is divided into 60 longitudinal zones, each spanning 6 degrees of longitude. These zones are numbered consecutively from west to east, starting at 180°W. Each zone is designated by a number from 1 to 60. Within each zone, a Transverse Mercator projection is used. This means that the surface of the earth is projected onto a cylinder that is tangent along a meridian within the zone.

UTM coordinates are typically referenced to a specific geodetic datum, such as WGS84 (World Geodetic System 1984), which provides a consistent reference surface for mapping. In the UTM system, locations on the Earth’s surface are specified using easting and northing coordinates, measured in meters. UTM provides relatively high accuracy for mapping and navigation within a specific zone. However, accuracy may decrease near the edges of zones due to distortion introduced by the projection. For a more detailed explanation about UTM coordinate system refer to [25,26,27].

The displacement of the UAV in meters is given by

$$\left[\begin{array}{c}{d}_x^{\left(UTM\right)}\\ d_y^{\left(UTM\right)}\end{array}\right]=\left[\begin{array}{c}s{d}^{\left(px\right)}cos\left(\psi \right)\\ s{d}^{\left(px\right)}sin\left(\psi \right)\end{array}\right]$$

(4)

where $d_x^{\left(UTM\right)}$ and $d_y^{\left(UTM\right)}$ represent the displacement in meters, $d^{\left(px\right)}$ is the displacement between the centroids in pixel units. $\psi$ is the yaw angle from the IMU and s the scale.

Once we have the last real latitude/longitude coordinates in UTM (meters), we can add the displacement of the centroids (in x and y direction), which has been already transformed into meters. Afterwards, the new UTM coordinate is converted back into latitude/longitude coordinates.

Up to this point, there are several issues involved. The first one is that the scale is calculated based on the altitude of the UAV relative to the terrain. Depending on the resolution of the elevation map and the resolution of the barometer, this scale may be computed incorrectly. If the scale is incorrect, then the conversion between pixels and meters will have an error, and the new coordinate will be calculated based on this incorrect scale.

The second problem is that this sequential algorithm relies on the previous location. If there was an error in the previous calculation, then we will be reading an incorrect elevation from the elevation map. If this reading is incorrect, then the scale for the next iteration will also have an error. In addition, if there are few keypoints, it is possible that we are computing the displacement based on some incorrect matches. In case the number of keypoints does not exceed a threshold, we will read the new yaw $\psi$ from the IMU, the distance between centroids $d^{\left(px\right)}$ that we had in the previous iteration (Figure 3).

Figure 8 shows the behavior of the algorithm using different values of $\delta$ to compute (2).

As expected, there is no single optimal value for the constant $\delta$ that is suitable for all flight segments.

4.2. Reduction of Accumulated Error

So far, the presented algorithm is able to estimate the UAV location, but the prediction error accumulates over time. To address this issue, we introduce a novel algorithm based on satellite images (geo-referenced) and the flight plan.

Using satellite imaging for localization is not a new idea [29,30,31]. In [29], the authors presented a method to index ortho-map databases based on image features. This method allow high speed search of geo-referenced images. On the other hand, in [30] an extension of the Cross-View Matching Network (CVM-Net) for aerial imaging was presented. This paper offers a Markov localization framework that enforces temporal consistency. Finally, in [31] a localization method for a multi-rotor is discussed, where an Autoencoder is used to compress images to an embedding space, enabling the inner product kernel to compare images.

The above-mentioned projects can achieve high accuracy for their intended tasks, but they are not suitable for a tactical fixed-wing UAV which can navigate around 200 Km. In the presented case we have the following limitations:

• Power consumption of the Hardware: We are allow to supply small embedded hardwares like a Raspberry Pi 4 [32] or an Orange Pi 5 [33], to not compromise the entire system. This prevents us from using many of the techniques based on deep learning architectures.
• Onboard processing: Working in GPS-Denied zones commonly imply noisy or denied communication, requiring the UAV to compute GPS estimation onboard.
• Altitude estimation: low altitude Multi-rotors can achieved high precision due to laser or ultrasonic sensors, but for a tactical fixed-wing UAV we have to estimate altitude with a barometer sensor, that can have up to hundred meter of error.
• High range missions: Tactical UAV can reach up to 200 Km range of operation. Then the algorithm has to work with long paths.

To reduce the cumulative error in an efficient way we present the following pre-processing scheme (Figure 9), we detail this scheme in the following subsections:

4.2.1. Calculate Flight Plan

We provide the user two options: they can either upload a previous flight plan (a list of latitude-longitude coordinates), or we can calculate a Flight plan based on the elevation map, a starting point, an end point, and a working altitude. Firstly, we load the elevation map as an image and operate this image to calculate the safe zone based on the working altitude of the plane.

After, we use the Rapidly-exploring random trees (RRT*) [34] to calculate the safe coordinates in the map. RRT* is a random algorithm that generates random extensions from the start point. RRT* can be used in the non structure environments It can be applied in non-structured environments if we can calculate intersections between a known point and a newly generated point.

In Figure 10, we show an example of the RRT* calculation over a real elevation map of a mountainous area in Hidalgo, Mexico. For better understanding, we include in Figure 11 a second example drawn in 3D for the same map but with different starting and ending points.

4.2.2. Filter Key-Points with the Flight Plan

Using a satellite image for reducing the accumulative error can be time-expensive, to avoid this, we apply a filter of key-points based on the flight plan. First we calculate the keypoints of a high-resolution image (this is very time-expensive but all this process is applied offline before the plane departs.

In contrast with the previous stage, where we used ORB, we are now implementing AKAZE [35] for feature extraction, which is an accelerated version of KAZE [36]. The reason for this choice is that AKAZE exhibits one of the best rotation invariance in contrast with other descriptors [24]. This is important because the UAV is not always heading north. The result of implementing AKAZE is shown in Figure 12.

The cardinality of key-points in Figure 12 (b) is 68061, this will imply a huge consumption of resources. To mitigate this, we filter these key-points by using the flight plan that we show in Figure 13 (a). In this case, the user provides a path instead of calculating it. Consider in Figure 14 (a), the histogram of the angular distance (in the latitude-longitude space) from each key-point to a point in the path (we use a linear interpolation of the path).

If we use a threshold of 0.005 over the distance to the path, we can obtain the filter key-points (Figure 13 (b)). In this case we count 25256 key-points, which is still a huge dataset to operate in real time. To have a better filter over the key-points we fit an exponential probability model, shown in (5), to the histogram of angular distance by using the maximum likelihood estimation (MLE) algorithm. In Figure 14 (b), we show the fitted exponential function.

$$f\left(x\right)=\lambda e^{\lambda x}$$

(5)

We use the fitted exponential function to randomly sample key-points (without replacement), where key-points closer to the flight plan are more probable to be sample. This process has the virtue that can choose the cardinality of the filtered key-points and then bound the time needed for correction.

In Figure 15 (a), we illustrate a simulation of the probability distribution over the map. Notice that the probability of using a particular site depends on this probability function and the quantity of key-points. In Figure 15 (b), we present the final filtered key-points, limited to 5000. We find that this number of key-points allow us to achieve a real time application.

4.2.3. Spacial Indexing

Once we get the filtered points, we want to calculate efficiently the closest point to the current position of the plane. Comparing the current position of the UAV with each possible satellite key-point is an $O\left(n\right)$ algorithm. This can be inefficient for real time applications.

the problem of rapidly find 2D or 3D points near a target point is known as spatial indexing [37]. This is a highly studied problem and there are some efficient solutions like K-dimensional tree (KD-tree) [38], Geohash [39], and Quad-tree [40]. In this paper, we implemented Quad-tree to get efficiently the K-nearest neighbors to search for a good match to reset position.

In Figure 16, we offer a general scheme for the Quad-tree algorithm. As its name suggest, Quad-tree is a tree data structure which each node has exactly four children. This allows to search points and query closest point in $O\left(log\right(n\left)\right)$. The tree separates a two dimensional space in four parts, each par can subsequently be separated in other four parts.

We have validated that Quad-tree is enough to archive real time for our particular case, but we reserve a comparison with other spatial indexing techniques for a future research.

Finally, if we find a satellite keypoint in the current image, we use (6) to apply the coordinates correction with the offset between the keypoint in pixel coordinates and the center of the image (Figure 17). Then we transform this offset to a Latitude/Longitude offset with GDAL.

$$\left[\begin{array}{c}\mathrm{Latitud}{\mathrm{e}}_{corr}\\ \mathrm{Longitud}{\mathrm{e}}_{corr}\end{array}\right]=\left[\begin{array}{c}\mathrm{Latitud}{\mathrm{e}}_{\left(kpt\right)}\\ \mathrm{Longitud}{\mathrm{e}}_{\left(kpt\right)}\end{array}\right]-\phi \left(s{R}_{\psi }\left[\begin{array}{c}x-c_x\\ y-c_y\end{array}\right]\right)$$

(6)

where $\mathrm{Latitud}{\mathrm{e}}_{\left(kpt\right)}$ and $\mathrm{Longitud}{\mathrm{e}}_{\left(kpt\right)}$ are the keypoint coordinates in the satellite image. x and y are the keypoint coordinates in the image. $c_x$ and $c_y$ are the image center. $\psi$ is the yaw orientation and $R_{\psi }$ is the rotation matrix, s the scale and $\phi \left(·\right)$ is the transformation from UTM to Latitude/Longitude.

Figure 18 is utilized to provide a concise overview of this algorithm. There is a preprocessing step where we obtain keypoints from the satellite image. If, during the flight, we find one of these keypoints, we correct our geolocalization, otherwise, we estimate the coordinates with the visual odometry phase.

5. Results

Figure 19 shows the results of the correction algorithm. As can be seen, the error is not incremental in contrast with previous approach where only visual odometry was used.

Since the output is very noisy, it is decided to filter it with a Kalman Filter and the output is shown in Figure 20.

Finally, we show the error plots of the first stage (visual odometry only) and the correction stage in Figure 21 and in Table 2.

6. Conclusions

In this work, an algorithm that estimates the geolocation of a UAV in the absence of GPS signal was presented. Several factors impact this problem, such as the altitude of the UAV relative to the terrain, the geolocation where the experiment takes place, and the distance of the route. It was demonstrated that there is no way to optimize the scale for pixel-to-meter transformation, so the error will accumulate as the aircraft moves forward. The results show that it is feasible to use satellite image features to extract landmarks that help the UAV estimate its location, even though the images are not updated, have different intensity levels, and varying resolutions.

The algorithm filters the keypoints detected in the satellite image to ensure computational lightness and onboard execution capability to ensure autonomy.

It is demonstrated that this scheme prevents the error to increment over time even in long range missions.