Unmanned Aerial Vehicle Route Planning in the Presence of a Threat Environment Based on a Virtual Globe Platform

Abstract

Route planning is a key technology for an unmanned aerial vehicle (UAV) to fly reliably and safely in the presence of a threat environment. Existing route planning methods are mainly based on the simulation scene, whereas approaches based on the virtual globe platform have rarely been reported. In this paper, a new planning space for the virtual globe and the planner is proposed and a common threat model is constructed for threats including a no-fly zone, hazardous weather, radar coverage area, missile killing zone and dynamic threats. Additionally, an improved ant colony optimization (ACO) algorithm is developed to enhance route planning efficiency and terrain masking ability. Our route planning methods are optimized on the virtual globe platform for practicability. A route planning system and six types of planners were developed and implemented on the virtual globe platform. Finally, our evaluation results demonstrate that our optimum planner has better performance in terms of fuel consumption, terrain masking, and risk avoidance. Experiments also demonstrate that the method and system described in this paper can be used to perform global route planning and mission operations.

1. Introduction

A UAV is an aircraft without a pilot on board that can be remotely controlled or flown automatically based on a pre-planned route or automation system [1]. With the development of the aviation electronics industry, UAVs play increasingly important roles in military and civil fields [2,3].

Generally, route planning for a UAV is an optimization problem that aims to generate a feasible route based on the tasks. The problem is an NP-hard problem [3].

For different types of tasks, scholars have selected different route planning methods. The vector field method is used in static or dynamic target tracking [4,5]. The rapidly-exploring random-tree (RRT) method has been applied to the path planning problem of indoor robots and mini UAVs [6,7,8]. Genetic algorithms (GA) are used to solve the travelling salesman problems related to UAVs, such as maximum information collection [9]. The evolutionary algorithm (EA) is used for multi-constraint route planning in a simulation scenario [10,11,12]. The particle swarm optimizer (PSO) is used to solve the path planning problem of UAVs on the sea [13]. Improved ACO [14,15], A* and Theta* [16] algorithms are used for route planning in three-dimensional environments. To improve the investigation efficiency of unmanned bombs, scholars have proposed the Quantum Wind Driven [17] algorithm.

Scholars have proposed a variety of optimization and improvement methods for the route planning problem under different conditions and have solved these problems relatively well. For different scenarios, each algorithm also has its own limitations: the ability of RRT to avoid obstacles is unsatisfactory; the processing time of the A* algorithm will increase explosively as the planning scene enlarges; and the computational complexity of GA and EA algorithms is high [18]. Scholars tend to optimize the selected algorithm according to their own simulation scenarios but the results of the experiments seem to not be very objective. For example, the results obtained by PSO were much better than GA in [13], whereas the results of PSO were inferior to GA in [18].

The virtual globe platform has the great advantages of low cost and ease-of-use in data collection, browsing, visualization and other aspects [19]. For the path planning problem for a high-endurance UAV, the virtual globe platform is a good choice to realize the modelling and visualization of a very large environment. There are many common virtual globe platforms such as Skyline, Google Earth, Virtual Earth, World Wind, and ArcGlobe.

Moreover, UAVs work in dangerous enemy territory in military penetration tasks. Avoiding threats from an enemy is a key factor in the success of tasks. Route planning for penetration finds a feasible route between the start point and end point in the presence of a threat environment. The defence of medium and high altitude areas in air defence systems is improving because of the development of the Radar Netting Technique [20]. There is no opportunity to penetrate without stealth aircraft. However, there are many radar blind zones in low-altitude areas because of topography and the curvature of the earth, and low-level flight becomes an important way to penetrate. On the virtual globe platform, threat modelling and route planning are more real and effective and mission operations can be performed [21]. Using the interactive capabilities of the virtual globe platform, operations such as parameter adjustment, route editing and storage, and flight simulation can be realized and successfully applied to industrial fields.

On the virtual globe platform, the method of model construction, planning space partitions and the realization of the algorithm will be different from the previous studies in the following ways:

• When modelling the radar threat areas, the maximum coverage range of early-warning radar can be up to hundreds, even thousands of kilometres, with larger signal coverage in high-altitude areas than at low-altitudes. By using the virtual globe platform, we can construct a more reasonable radar threat model according to the radar equation and fully consider the influences of earth curvature and terrain masking.
• The scale of the terrain data is large on the virtual globe platform. In view of this problem, this paper proposes a multi-granularity planning space to achieve a balance between accuracy and efficiency.
• In low-altitude penetration, a good valley-following ability can effectively avoid a radar threat, including unknown radar threats. We propose a strengthened local valley-following algorithm for route planning.
• The planning space needs to be transformed between Cartesian systems and Geodetic systems. Because of the overhead problem generated by large-scale data and space transformation, we optimize some of the algorithm’s implementation details.
• Because there are certain errors in a UAV's navigation system and control system, this paper refers to industry standards, such as the performance based navigation (PBN) standard [22], to optimize and improve the robustness of the route to avoid collisions because of flight errors.

Section 2 describes the route planning problem and evaluation indexes in the risk environment. In Section 3 we propose a multi-granularity planning space on the virtual globe platform for route planning and consider multiple types of threats. In Section 4, we propose an improved ACO algorithm with valley-following and threat avoidance for route planning and some route optimization algorithms to make the route more effective and robust. Section 5 shows the practicability of the proposed methods through experiments. Section 6 presents conclusions and proposes future research work.

2. Route Planning Problem Descriptions

The location of a point $\mathrm{P}$ in geodetic space ${\mathrm{\Psi }}^3$ is described by longitude $\mathrm{x}$, latitude $\mathrm{y}$, and altitude $\mathrm{H}$. A UAV flies along the designated route point sequence ${\mathrm{R}}_{\mathrm{route}}=\left\{{\mathrm{r}}_1{,\mathrm{r}}_2,\cdots {,\mathrm{r}}_{\mathrm{n}}\right\},{\mathrm{r}}_{\mathrm{i}}\in {\mathrm{\Psi }}^3$, where ${\mathrm{r}}_1$ is the start point and ${\mathrm{r}}_{\mathrm{n}}$ is the end point. The goal of penetration routing is to obtain a feasible route with minimal fuel consumption, maximum terrain masking, and minimum risk.

Both the navigation systems and the control systems of UAVs contain certain deviations. The international civil aviation organization proposed the concept of the required navigation performance (RNP). RNP [22] describes the precision that can be attained during at least 95% of the flight time; the precision unit is the nautical mile (nmi). In Figure 1, the segment width of the route is defined as 4 × RNP. The planned route may not satisfy the performance constraints of UAVs.

We evaluate the planning route on four dimensions: minimal fuel consumption, maximum terrain masking, minimum risk and performance safety. Additionally, we evaluate the efficiency of the planner:

• ${\mathrm{L}}_{\mathrm{route}}$: Route length. Connect the points in ${\mathrm{R}}_{\mathrm{route}}$ successively and obtain a series of line segments, then calculate the total length of the segments to obtain the route length. To evaluate the result, we introduce the theoretically shortest length of the route, ${\mathrm{L}}_{\min }$. Then, we connect the ${\mathrm{r}}_1$ and ${\mathrm{r}}_{\mathrm{n}}$ of ${\mathrm{R}}_{\mathrm{route}}$ to obtain the straight line segment ${\mathrm{r}}_1{\mathrm{r}}_{\mathrm{n}}$. ${\mathrm{L}}_{\min }$ is the length of ${\mathrm{r}}_1{\mathrm{r}}_{\mathrm{n}}$.
• ${\mathrm{h}}_{\mathrm{avT}}$: Average terrain altitude passed by the route. ${\mathrm{h}}_{\mathrm{avT}}$ is calculated as follows: Obtain the interpolation points ${\{\mathrm{r}}_1{(\mathrm{x}}_1{,\mathrm{y}}_1{,\mathrm{H}}_1{),\mathrm{r}}_2{(\mathrm{x}}_2{,\mathrm{y}}_2{,\mathrm{H}}_2),\cdots {,\mathrm{r}}_{\mathrm{k}}{(\mathrm{x}}_{\mathrm{k}}{,\mathrm{y}}_{\mathrm{k}}{,\mathrm{H}}_{\mathrm{k}})\}$ with interval $\forall$ of ${\mathrm{R}}_{\mathrm{route}}$. We can then obtain points on the terrain ${\{\mathrm{r}}_1{(\mathrm{x}}_1{,\mathrm{y}}_1{,\mathrm{H}}_{1\mathrm{T}}{),\mathrm{r}}_2{(\mathrm{x}}_2{,\mathrm{y}}_2{,\mathrm{H}}_{2\mathrm{T}}),\cdots {,\mathrm{r}}_{\mathrm{k}}{(\mathrm{x}}_{\mathrm{k}}{,\mathrm{y}}_{\mathrm{k}}{,\mathrm{H}}_{\mathrm{kT}})\}$, and ${\mathrm{h}}_{\mathrm{avT}}={\sum }_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{H}}_{\mathrm{iT}}/\mathrm{k}$. To evaluate the result, we introduce a reference ${\mathrm{H}}_{\mathrm{avT}}$, the average terrain altitude passed by ${\mathrm{r}}_1{\mathrm{r}}_{\mathrm{n}}$. We can understand the terrain masking ability of the planner by contrasting the result with ${\mathrm{H}}_{\mathrm{avT}}$.
• ${\mathrm{L}}_{\mathrm{FA}}$: Route length that passes through the threat zones.
• ${\mathrm{L}}_{\mathrm{US}}$: Route length that does not meet the flight performance safety requirements.
• Processing time.

3. Build the Planning Space

3.1. Virtual Globe Space

For the penetration tasks on the virtual globe, the planning space needs to be converted between Geodetic space ${\mathrm{\Psi }}^3$ and Cartesian space ${\mathrm{Ʊ}}^3$. ${\mathrm{\Psi }}^3$ treats the earth as a reference ellipsoid. $\mathrm{P}(\mathrm{x},\mathrm{y},\mathrm{H})\in {\mathrm{\Psi }}^3$ can be converted into $Q\left(\mathrm{X},\mathrm{Y},\mathrm{Z}\right)\in {\mathrm{Ʊ}}^3$ as follows:

$$\{\begin{array}{l}\mathrm{X}=\left(\mathrm{M}+\mathrm{H}\right)\cos \left(\mathrm{y}\right)\cos \left(\mathrm{x}\right)\\ \mathrm{Y}=\left(\mathrm{M}+\mathrm{H}\right)\cos \left(\mathrm{y}\right)\sin \left(\mathrm{x}\right)\\ \mathrm{Z}=\left[\mathrm{M}\left(1-\frac{{\mathrm{a}}^2-{\mathrm{b}}^2}{{\mathrm{a}}^2}\right)+\mathrm{H}\right]\sin \left(\mathrm{y}\right)\end{array}$$

(1)

where $\mathrm{M}$ represents the radius of curvature of the reference ellipsoid, computed as follows:

$$\mathrm{M}=\frac{\mathrm{a}}{\sqrt{1-\frac{{\mathrm{a}}^2-{\mathrm{b}}^2}{{\mathrm{a}}^2}{\sin }^2(\mathrm{y})}}$$

(2)

where $\mathrm{a}$ and $\mathrm{b}$ represent the length of the major axis and the minor axis of the reference ellipsoid, respectively.

$\mathrm{Q}(\mathrm{X},\mathrm{Y},\mathrm{Z})\in {\mathrm{Ʊ}}^3$ can be converted into $\mathrm{P}(\mathrm{x},\mathrm{y},\mathrm{H})\in {\mathrm{\Psi }}^3$ as follows:

$$\{\begin{array}{l}\mathrm{x}=\arctan (\frac{\mathrm{Y}}{\mathrm{X}})\\ \mathrm{y}=\arctan (\frac{\mathrm{Z}(\mathrm{M}+\mathrm{H})}{\sqrt{{(\mathrm{X}}^2+{\mathrm{Y}}^2)\left[\mathrm{M}\left(1-\frac{{\mathrm{a}}^2-{\mathrm{b}}^2}{{\mathrm{a}}^2}\right)+\mathrm{H}\right]}})\\ \mathrm{H}=\frac{\sqrt{{\mathrm{X}}^2+{\mathrm{Y}}^2}}{\cos \left(\mathrm{y}\right)}-\mathrm{M}\end{array}$$

(3)

The reference ellipsoid used by the Global Positioning System and virtual globe platforms is the World Geodetic System of 1984 ellipsoid (WGS-84 ellipsoid) [22].

The lengths of the same longitude interval in different latitudes and the same latitude interval in different longitudes are not fixed in ${\mathrm{\Psi }}^3$. For example, the distance between (90°, 29°, 4000 m) and (91°, 29°, 4000 m) is 97,422.04 m, whereas the distance between (90°, 26°, 4000 m) and (91°, 26°, 4000 m) is 100,114.77 m. In the later sections, the planning space is in ${\mathrm{\Psi }}^3$, but the distance measure, interpolation algorithm, equation calculation and other operations are carried out in ${\mathrm{Ʊ}}^3$.

The terrain is usually described by digital elevation model (DEM) data, whereas terrain texture is described by digital orthophoto map (DOM) data. The data of DEM and DOM can be represented as ${\{(\mathrm{x},\mathrm{y},\mathrm{z})}_{\mathrm{i}}\}$, where $\mathrm{x}$ is longitude and $\mathrm{y}$ is latitude. In DEM data, $\mathrm{z}$ is the altitude value of the corresponding position, whereas $\mathrm{z}$ is the pixel value for DOM data, as Figure 2a,b illustrate.

DEM and DOM describe discrete space, whereas the real world belongs to contiguous space. The information of any position can be obtained by using a spatial interpolation method through the adjacent points, and the bicubic interpolation algorithm is described as follows:

$${\mathrm{Z}}_{\mathrm{p}}=\mathrm{f}({\mathrm{x}}_{\mathrm{p}},{\mathrm{y}}_{\mathrm{p}})={{\displaystyle \sum }}_{\mathrm{i}=0}^3{\displaystyle \sum }_{\mathrm{j}=0}^{3-\mathrm{i}}{\mathrm{a}}_{\mathrm{ij}}{\mathrm{x}}_{\mathrm{p}}^{\mathrm{i}}{\mathrm{y}}_{\mathrm{p}}^{\mathrm{j}}$$

(4)

We can use the adjacent 4 × 4 points of point $\mathrm{A}$ in Figure 2a to calculate the coefficients of the interpolation function to determine the altitude value of point $\mathrm{A}$.

In Figure 2c, we can realize the visualization of 3D terrain by superimposing the DOM data on the DEM data. In this study, we used the improved NASA-WorldWind platform; the data source was the SRTM DEM data at 90-m resolution and the Landsat DOM data at 60-m resolution provided publicly by NASA, as shown in Figure 2d.

3.2. Threat Modelling

3.2.1. General Types of Threats

No-fly zones, hazardous weather, high-rise buildings and low-altitude control zones need to be considered when planning routes.

No-fly zones are regions that a UAV cannot fly into. A no-fly zone can be described as the enclosed area defined by the point set $\mathrm{A}=\{{\mathrm{A}}_1,{\mathrm{A}}_2,\cdots ,{\mathrm{A}}_{\mathrm{m}}\}$, m > 2, as shown in Figure 3a.

High buildings, slowly hazardous weather and low-altitude control zones are low-altitude regions that UAVs cannot fly into. Such threats can be described as $\{\{{\mathrm{A}}_1,{\mathrm{A}}_2,{\mathrm{A}}_3,\cdots ,{\mathrm{A}}_{\mathrm{m}}\},{\mathrm{H}}_{\max }\}$, m > 2, where ${\mathrm{H}}_{\max }$ is the upper bound. This type of threat is modelled as observed in Figure 3b,c. These threats can be dealt with as terrain data; we transform this type of threat area into terrain before route planning.

3.2.2. Killing Zone of an Air Defense Missile

The killing zone [23,24] is an important guideline for judging the campaign performance of Air Defense Missile Weapon Systems. A killing zone describes an area of space in which a missile can destroy a target at no less than a certain probability once the target enters the area. The mathematical model of a vertical killing zone is shown in Figure 4a. Horizontal lines define the high and low boundaries with heights of ${\mathrm{H}}_{\max }$ and ${\mathrm{H}}_{\min }$, respectively. ${\mathrm{D}}_{\mathrm{symin}}$ is the minimum slant range of the far boundary of the killing zone and ${\mathrm{D}}_{\mathrm{symax}}$ is the maximum slant range. The minimum slant range and the maximum height angle of the near boundary of the killing zone are ${\mathrm{D}}_{\mathrm{sjmin}}$ and ${\mathrm{A}}_{\max }$, respectively. The near boundary and far boundary are the arcs with $\mathrm{O}$ as the centre. Figure 4b shows the 3D visualization model of the missile killing zone.

3.2.3. Radar Threat Space

Radar plays a very important role in modern air defence systems; it is a key threat to the penetration of UAVs that needs to be avoided by route planning. Affected by the terrain masking and the earth curvature, it is difficult for radar to find a target in low-altitude flight.

The radar threat model should refer to the radar detection probability threshold that a UAV can maximally tolerate. Given the detection probability threshold and radar performance parameters, the maximum range of radar ${\mathrm{R}}_{\max }$ can be estimated by the radar equation [25]. Assuming that the radar is deployed at ${\mathrm{R}}_{\mathrm{a}}{(\mathrm{x}}_{\mathrm{a}}{,\mathrm{y}}_{\mathrm{a}}{,\mathrm{H}}_{\mathrm{a}})$ in ${\mathrm{\Psi }}^3$, we transform it into ${\mathrm{Ʊ}}^3$ as ${\mathrm{R}}_{\mathrm{a}}{(\mathrm{X}}_{\mathrm{a}}{,\mathrm{Y}}_{\mathrm{a}}{,\mathrm{Z}}_{\mathrm{a}})$, with the direction of pitch angle $\mathsf{\theta }$ and azimuth angle $\mathsf{\beta }$, and we can describe ${\mathrm{p}}_{\mathrm{k}}\left({\mathrm{X}}_{\mathrm{k}}{,\mathrm{Y}}_{\mathrm{k}}{,\mathrm{Z}}_{\mathrm{k}}\right)$ on the boundary of the radar coverage as:

$$\{\begin{array}{l}{\mathrm{X}}_{\mathrm{k}}={\mathrm{X}}_{\mathrm{a}}={\mathrm{R}}_{\max }{\mathrm{f}}^{\prime }\left(\mathsf{\theta }\right)\cos \left(\mathsf{\beta }\right)\\ {\mathrm{Y}}_{\mathrm{k}}={\mathrm{Y}}_{\mathrm{a}}={\mathrm{R}}_{\max }{\mathrm{f}}^{\prime }\left(\mathsf{\theta }\right)\sin \left(\mathsf{\beta }\right)\\ {\mathrm{Z}}_{\mathrm{k}}={\mathrm{Z}}_{\mathrm{a}}={\mathrm{R}}_{\max }{\mathrm{f}}^{\prime }\left(\mathsf{\theta }\right)\end{array}$$

(5)

where ${\mathrm{f}}^{\prime }(\mathsf{\theta })$ describes the beam pattern of the radar antenna in the direction of the pitch angle.

We use the Gaussian function to approximate the beam pattern as:

$${\mathrm{f}}^{\prime }(\mathsf{\theta })={\mathrm{e}}^{-4\ln \sqrt{2}{\mathsf{\theta }}^2{/\mathsf{\theta }}_{\mathrm{b}}}$$

(6)

where ${\mathsf{\theta }}_{\mathrm{b}}$ is the signal beam width.

The radar threat curved surface can thus be obtained. The probability of radar detection on this surface is constant. Outside the curved surface, the radar detection probability is less than the probability threshold and the UAV is regarded as safe. Within the curved surface, the radar detection probability is higher than the probability threshold and the UAV is not safe.

In low-altitude areas, radar signals may be masked by terrain. In Figure 5a, radar beams are partly masked by terrain and the shaded area represents the area that signals can reach. The visualization of the radar beam can be constructed as follows:

• Calculate the detection range of the radar with different sampled pitch angles by using Equations (5) and (6); this produces a series of sampling points such as $\mathrm{A}–\mathrm{G}$, as shown in Figure 5a.
• Connect the radar centre $\mathrm{O}$ to the sampling point and obtain a series of sampling points on this straight line by sampling with short intervals. Then, we transform these sampling points to ${\mathrm{\Psi }}^3$ and compare the altitudes of the sampling points with the terrain. Once the sampling point is below the terrain, we can determine that the subsequent region cannot be reached by the radar signal. In Figure 5a, according to the terrain visibility analysis of $\mathrm{OD}$, the boundary point $\mathrm{D}$ is moved to $\mathrm{M}$.

If we calculate the radar coverage directly, the earth curvature will introduce considerable error. As shown in Figure 5b, $\mathrm{O}$ is the earth centre, $\mathrm{BD}$ is the horizontal plane, the circular arc $\mathrm{CE}$ is the geoid, and point $\mathrm{B}$ is in the same horizontal plane as point $\mathrm{A}$. When $\mathrm{B}$ is very near to $\mathrm{A}$, the altitudes of $\mathrm{A}$ and $\mathrm{B}$ are approximately the same and may be ignored. However, when $\mathrm{B}$ is far from $\mathrm{A}$, the elevation difference between $\mathrm{A}$ and $\mathrm{B}$ due to the influence of the earth curvature, $\Delta \mathrm{H}$ in Figure 5b, cannot be ignored. $\mathrm{L}$ represents the length of $\mathrm{AB}$, $\mathrm{R}$ is the average curvature radius of the earth, and $\Delta \mathrm{H}$ can be calculated by:

$$\Delta \mathrm{H}=\sqrt{{\mathrm{R}}^2+{\mathrm{L}}^2}-\mathrm{R}$$

(7)

When the distance between two points is 100 km, the elevation difference is approximately 785 m because of the influence of the curvature. First, we need to convert from ${\mathrm{\Psi }}^3$ to ${\mathrm{Ʊ}}^3$ when calculating the boundary of the radar coverage. Then, we can calculate the radar coverage and perform visibility analysis. Finally, the results are converted back to ${\mathrm{\Psi }}^3$ and the errors caused by the curvature of the earth have been eliminated.

Figure 6 shows the 3D visualization of the radar network coverage under the influence of the earth’s curvature.

3.3. Dynamic Threat

The models presented in Section 3.2 can simulate most of the threats in normal circumstances. Sometimes the threat movement needs to be considered, for example, the forecast of hazardous weather such as a typhoon (hurricane) and other predictable dynamic threat areas. At this time, we need to introduce the dynamic threat model. Based on the general threat model, a time stamp $\mathrm{t}$ is introduced, the property of a dynamic threat at $\mathrm{t}$ is described by ${\mathrm{threat}}_{\mathrm{t}}=\{\{{\mathrm{A}}_1{,\mathrm{A}}_2{,\mathrm{A}}_3,\cdots {,\mathrm{A}}_{\mathrm{m}}{\},\mathrm{H}}_{\max },\mathrm{t}\},\mathrm{m}>2$. A full path dynamic threat can be described by $\left\{{\mathrm{threat}}_{{\mathrm{t}}_1}{,\mathrm{threat}}_{{\mathrm{t}}_2},\cdots {,\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{n}}}\right\}$, where the threat begins at ${\mathrm{t}}_1$ and disappears at ${\mathrm{t}}_{\mathrm{n}}$. As shown in Figure 7a,b, if we want to know the state of a threat at any time ${\mathrm{t}}_{\mathrm{k}}$, we can calculate it as follows:

• Search the interval [${\mathrm{t}}_i,{\mathrm{t}}_{i+1}]$ that satisfies ${\mathrm{t}}_{\mathrm{i}}\le {\mathrm{t}}_{\mathrm{k}}<{\mathrm{t}}_{\mathrm{i}+1}$ to find the ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}}}$ and ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}+1}}$. If not found, then there is no threat area at this time.
• Given the threat ${\mathrm{t}}_{\mathrm{k}}=2.4$ shown in Figure 7b, use linear interpolation for ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}}}$ and ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}+1}}$ to obtain ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{k}}}$. As shown in Figure 7b, we can use linear interpolation for the four vertices of ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}}}$ and ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{i}+1}}$ to obtain ${\mathrm{threat}}_{{\mathrm{t}}_{\mathrm{k}}}$.

3.4. Build the Grid Planning Space

Route planning on the virtual globe platform in this study is based on a graph search algorithm. We define a planning graph composed of nodes and edges. The graph needs to effectively express threats and terrain features.

In low-altitude penetration missions, threat models need to be marked. However, real-time sampling of a threat body at any point in the planning space will dramatically increase the delay time of the planner. This paper introduces a type of grid space with location information, elevation information, and other attribute information.

In Figure 8, the space is divided into a two-dimensional grid by equal intervals of latitude and longitude; each node in the grid has eight extension nodes. However, the actual distance represented by the equal intervals of latitude and longitude is different. In the virtual globe platform, the planning space is actually an irregular eight-connected grid space.

A UAV flies at the safe height $\mathrm{h}$ by default in this space. Each node stores the information including ${[\mathrm{lon},\mathrm{lat},\mathrm{alt},\mathrm{H}}_{\min }]$, where $\mathrm{lon}$ and $\mathrm{lat}$ represent the longitude and latitude of the node, respectively, $\mathrm{alt}$ represents the terrain altitude superposed by the high buildings and low-altitude control zones, and ${\mathrm{H}}_{\min }$ is the lowest bound of the radar coverage area or the missile killing zone in this area. If the node $\mathrm{h}$ that needs to be extended is higher than ${\mathrm{H}}_{\min }$, the node will be eliminated. As the figure shows, nodes $\mathrm{A}$, $\mathrm{B}$, $\mathrm{C}$ and $\mathrm{D}$ are in a radar coverage area. However, because of the terrain masking effect, points $\mathrm{B}$ and $\mathrm{C}$ can be used as planning candidate nodes with higher ${\mathrm{H}}_{\min }$, whereas nodes $\mathrm{A}$ and $\mathrm{D}$ with low ${\mathrm{H}}_{\min }$ will be directly eliminated. As the figure shows, node $\mathrm{E}$, $\mathrm{F}$, $\mathrm{G}$ and $\mathrm{H}$ are in a missile threat area and will be eliminated because of their low ${\mathrm{H}}_{\min }$. If the current point is located in the no-fly zone, ${\mathrm{H}}_{\min }$ is set to 0 and the point will be eliminated by the planner. Because of the planning space generated beforehand, the planner can quickly eliminate some infeasible nodes.

3.4.1. Threat Area Expansion

Considering the navigation error and the route optimization error, as shown in Figure 9a, the planned route in grid space may intersect the threatened area. We expand the boundary of the threat models when mapped to the planning space. The low boundary of the radar area and the missile killing zone are also reduced. Figure 9b shows that the route is safe after expanding the threat area.

3.4.2. Multi-Granularity Planning Space

The finer granularity of the grid space is because the more nodes in planning space, the more precise the planning results will be, although the planning will take more time. To obtain a balance between performance and precision for different application scenarios and navigation performance parameters of a UAV, different grid sizes need to be used.

We need resampling to construct different sizes of the grid using DEM data. Common resampling methods are Box Splines [26], interpolation [27] and the discrete wavelet transform [28]. This paper uses the bicubic interpolation method of Equation (4) to resample. We generate grid space with multi-granularity beforehand and choose different grid sizes according to the requirements of different tasks.

4. Route Planning Method

We propose an improved ACO-based planner that considers the local terrain environment to obtain the initial route. An ACO algorithm simulates the foraging process of ants in principle, and its developer used it to efficiently solve the Travelling Salesman Problem [29,30]. The pheromone model and the probability model are the core of this algorithm [31]. ACO algorithms have been widely used for pattern classification [32], cloud computing [33], network coding [34], robot path planning [35] and in other fields. We propose some optimizations to improve the ACO algorithm for the study of the virtual globe platform and penetration route planning.

4.1. Basic ACO Algorithm [29,30]

The ant colony begins searching for food without being told where the food is. In the process of moving, ants release pheromones into the environment. When there is no pheromone in the environment, the ants perform a random walk, and when the pheromone concentration is high, the ants will move along the pheromone path with a greater probability. If the ant finds food, it will move along its pheromone path to return to its nest. Along the shorter path, the ants make the roundtrip more quickly, which is more likely to attract more ants to walk along this path. The pheromone concentration will thus be further enhanced and eventually the ant colony will walk along the shortest path. The basic ACO algorithm simulates the pathfinding and feedback processes of the ant colony.

4.1.1. Roulette Wheel Selection Model

The ACO algorithm uses a roulette wheel selection model to simulate the walk choice of ants in nature. This model has a larger probability to choose adjacent nodes with higher pheromone concentration; it allows ants to make mistakes that have a small probability and retain the chance to find a better path. In the model, the ant colony has $\mathrm{m}$ ants. The transfer probability ${\mathrm{p}}_{\mathrm{ij}}^{\mathrm{k}}(\mathrm{t})$ of an ant $\mathrm{k}(0\le \mathrm{k}<\mathrm{m})$ moving from node $\mathrm{i}$ to an adjacent node $\mathrm{j}$ in the $\mathrm{t}$ round of iteration is expressed as:

$${\mathrm{p}}_{\mathrm{ij}}^{\mathrm{k}}(\mathrm{t})=\{\begin{array}{l}\frac{{\mathsf{\tau }}_{\mathrm{ij}}^{\mathsf{\alpha }}{(\mathrm{t})\mathsf{\eta }}_{\mathrm{ij}}^{\mathsf{\beta }}(\mathrm{t})}{{\sum }_{\mathrm{s}\in {\mathrm{J}}_{\mathrm{k}}(\mathrm{i})}{\mathsf{\tau }}_{\mathrm{is}}^{\mathsf{\alpha }}(\mathrm{t}){\mathsf{\eta }}_{\mathrm{is}}^{\mathsf{\beta }}(\mathrm{t})},\mathrm{if}\mathrm{j}\mathsf{ϵ}{\mathrm{J}}_{\mathrm{k}}(\mathrm{i})\\ 0,\mathrm{otherwise}\end{array}$$

(8)

where ${\mathsf{\tau }}_{\mathrm{ij}}\left(\mathrm{t}\right)$ is the residual pheromone of the edge $<\mathrm{i},\mathrm{j}>$. ${\mathsf{\eta }}_{\mathrm{ij}}\left(\mathrm{t}\right)$ is the heuristic function $\mathrm{from}$ node $\mathrm{i}$ to node $\mathrm{j}$; when solving the TSP, ${\mathsf{\eta }}_{\mathrm{ij}}\left(\mathrm{t}\right)={1/\mathrm{c}}_{\mathrm{ij}}$. ${\mathrm{c}}_{\mathrm{ij}}$ reflects the cost of the walk from node $\mathrm{i}$ to node $\mathrm{j}$. $\mathsf{\alpha }$ and $\mathsf{\beta }$ are the weights of the pheromone and the heuristic function, respectively. Adjacent nodes ${\mathrm{J}}_{\mathrm{k}}\left(\mathrm{i}\right)$ can be chosen by the ant $\mathrm{k}$ in node $\mathrm{i}$. According to the planning space defined above, each node has eight expanding nodes.

4.1.2. Pheromone Updating Model

The pheromone updating model provides a feedback mechanism in the ACO algorithm. After all the ants have completed an iteration round, pheromones will be left on the path traversed by the ants and the pheromone concentration is updated as follows:

$${\mathsf{\tau }}_{\mathrm{ij}}\left(\mathrm{t}+1\right)={\mathsf{\rho }\mathsf{\tau }}_{\mathrm{ij}}\left(\mathrm{t}\right)+{\displaystyle \sum }_{k=1}^m\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{k}}(\mathrm{t})$$

(9)

where $\mathsf{\rho }$ is the pheromone retention coefficient $(0\le \mathsf{\rho }<1)$. $\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{k}}(\mathrm{t})$ represents the pheromones left at the edge of $<\mathrm{i},\mathrm{j}>$ by the ant $\mathrm{k}$ at the iteration of round $\mathrm{t}$ and is calculated as follows:

$$\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{k}}\left(\mathrm{t}\right)=\{\begin{array}{l}\frac{1}{{\mathrm{L}}^{\mathrm{k}}},\mathrm{if}<\mathrm{i},\mathrm{j}>\in {\mathrm{Q}}_{\mathrm{k}}\\ 0,\mathrm{otherwise}\end{array}$$

(10)

where ${\mathrm{Q}}_{\mathrm{k}}$ is the set of edges traversed by ant $\mathrm{k}\mathrm{and}{\mathrm{L}}^{\mathrm{k}}$ is the total cost of the edges.

4.2. Parameter Optimization

In the basic ACO algorithm, the ants do not know the target location. Our application scenario knows the location of the target and we can thus introduce a heuristic function to guide the ants to accelerate the iterative process of the algorithm. In the pheromone feedback mechanism of the basic ACO algorithm, the pheromone left by the less successful ants may interfere with a better result, which easily causes the algorithm to fall into a local optimum. To solve the UAV problem described in this paper, we improved the heuristic cost function, pheromone update mechanisms, algorithm efficiency and other aspects of the algorithm. Additionally, to avoid the ants returning, we set up tabu lists to mark the nodes that the ants have traversed; the expanding nodes located in the threat area or in the tabu lists are excluded.

4.2.1. Cost Function

In view of the planning space built in Chapter 3, we introduce the cost function:

$$\{\begin{array}{l}\mathrm{W}={\displaystyle {\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}}{\mathrm{w}}_{\left(\mathrm{i}-1\right)\mathrm{i}}\\ {\mathrm{w}}_{\mathrm{ij}}=\frac{{\mathrm{l}}_{\mathrm{ijtraj}}}{{\mathrm{l}}_{\max }}+\epsilon \frac{{\mathrm{h}}_{\mathrm{ijtraj}}-{\mathrm{h}}_{\min }}{{\mathrm{h}}_{\max }-{\mathrm{h}}_{\min }}\end{array}$$

(11)

where ${\mathrm{w}}_{\mathrm{ij}}$ is the cost of the edge $<\mathrm{i},\mathrm{j}>$, ${\mathrm{l}}_{\mathrm{ijtraj}}$ is the Euclidean length of the edge $<\mathrm{i},\mathrm{j}>$, which approximately reflects the fuel consumption cost of the UAV, and ${\mathrm{h}}_{\mathrm{ijtraj}}$ is the average terrain elevation of node $\mathrm{i}$ and $\mathrm{j}$, which is ${(\mathrm{alt}}_{\mathrm{i}}+{\mathrm{alt}}_{\mathrm{j}})/2$ and reflects the altitude of the route. Because the units of the fuel consumption cost and height cost are not of the same order of magnitude, they are normalized. ${\mathrm{l}}_{\max }$ is the largest Euclidean distance of the two adjacent nodes in the grid. Because planning space is irregular, we set ${\mathrm{l}}_{\max }=\mathrm{grid}\mathrm{width}\times 160\mathrm{KM}$ when the grid width is determined. ${\mathrm{h}}_{\max }$ and ${\mathrm{h}}_{\min }$ are the highest and lowest elevation of the planning space, respectively. $\epsilon$ is the weight, which determines the valley-following performance of the planner and enables the UAV to fly as low as possible. Flying at low-altitude can improve the survival rate of the UAV if the deployment information of the enemy radar is unknown.

4.2.2. Heuristic Function and Valley-Following

In the basic ACO algorithm, the ants have difficulty reaching the target successfully because the ${\mathrm{c}}_{\mathrm{ij}}$ can only reflect the cost of node $\mathrm{i}$ to $\mathrm{j}$, and thus the direction to the end node is unknown. In this study, a heuristic cost function is introduced that can promote the ant to move towards the target:

$${\mathsf{\eta }}_{\mathrm{ij}}\left(\mathrm{t}\right)=2+\frac{{\mathrm{T}}_{\mathrm{i}}-{\mathrm{T}}_{\mathrm{j}}}{{\mathrm{l}}_{\max }}$$

(12)

where ${\mathrm{T}}_{\mathrm{j}}$ is the distance between node $\mathrm{j}$ and the target point and ${\mathrm{T}}_{\mathrm{i}}$ is the distance between node $\mathrm{i}$ and the end point. If ${\mathrm{T}}_{\mathrm{j}}{>\mathrm{T}}_{\mathrm{i}}$, the ants are moving farther and farther away from the target and will be assigned a smaller probability to select this node.

Equation (12) can promote the ants to walk to the destination node; however, the ants cannot use the local valley terrain information to achieve valley-following. When we analysed the terrain data, we found that the valley region was continuous, and the enhanced local valley-following algorithm is proposed accordingly. The eight adjacent nodes of the current point $\mathrm{i}$ are $\mathrm{A}–\mathrm{H}$, and the implementation method is as follows:

• Calculate the average terrain altitude of $\mathrm{A}–\mathrm{H}$ points.
• Compare the average altitude with the altitude of point $\mathrm{i}$; if greater than the threshold $\mathsf{\delta }$, the UAV at point $\mathrm{i}$ is located in the mountainous area.
• If the UAV is located in the mountainous region, the adjacent $\mathrm{A}–\mathrm{H}$ points are sorted according to their altitudes, and the $2\epsilon$ points with highest elevation are set to ${\mathsf{\eta }}_{\mathrm{ij}}\left(\mathrm{t}\right)=0$.
• If the UAV is located in the non-mountainous region, all the neighbouring nodes are reachable.

After this treatment, ${\mathsf{\eta }}_{\mathrm{ij}}$ can promote movement of the ants to the target point and make full use of local terrain information.

4.2.3. Pheromone Update Mechanisms

The updating mechanism of the basic ACO algorithm makes it easy for premature results to occur. The pheromone update strategy we use is as follows:

$${\mathsf{\tau }}_{\mathrm{ij}}\left(\mathrm{t}+1\right)={\mathsf{\rho }\mathsf{\tau }}_{\mathrm{ij}}\left(\mathrm{t}\right)+\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{best}*}(\mathrm{t})$$

(13)

where $\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{best}*}(\mathrm{t})$ is the pheromone increment of global optimal ant:

$$\Delta {\mathsf{\tau }}_{\mathrm{ij}}^{\mathrm{best}*}\left(\mathrm{t}\right)=\{\begin{array}{l}\frac{1}{{\mathrm{W}}^{\mathrm{best}*}},\mathrm{if}<\mathrm{i},\mathrm{j}>\in {\mathrm{Q}}_{\mathrm{best}}\\ 0,\mathrm{otherwise}\end{array}$$

(14)

where ${\mathrm{Q}}_{\mathrm{best}}$ is the set of edges traversed by the global optimal ant. ${\mathrm{W}}^{\mathrm{best}*}$ was calculated according to Equation (11) and is the cost of the global optimal route. To ensure that the nodes that the global optimal ant did not walk can be reached with a certain probability and to prevent premature results because of too much pheromone on the global optimal route, we limit the pheromone concentration to ${[\mathsf{\tau }}_{\min },{\mathsf{\tau }}_{\max }]$.

$$\{\begin{array}{l}{\mathsf{\tau }}_{\max }=\frac{1}{1-\mathsf{\rho }}\frac{1}{{\mathrm{W}}^{\mathrm{best}*}}\\ {\mathsf{\tau }}_{\min }=\frac{{\mathsf{\tau }}_{\max }}{\mathrm{n}}\end{array}$$

(15)

where $\mathrm{n}$ is the estimated node number of the global optimal route.

4.2.4. Algorithm Efficiency Improvement

Because of the introduction of the virtual globe platform, every distance calculation needs to transform between two coordinate systems. To reduce the time cost of the repeated distance computation of the ACO, we calculate the distance between the adjacent nodes in advance. This paper introduces a three-dimensional array of off-line storage; the third dimension of the array is used to store the distance from the current node $\mathrm{i}$ to the eight adjacent nodes. Additionally, a two-dimensional array is used to store the distance between $\mathrm{i}$ and the end node.

In the ACO algorithm, the worst performing ant may run a very biased path with a relatively small probability, thus affecting the overall efficiency of the algorithm. In this paper, the maximum number of search steps is $10\mathrm{n}$; when the search steps of ants exceed the maximum number, the result is discarded.

4.3. Dynamic Threat Avoidance

When considering dynamic threats, we need to introduce a dynamic threat avoidance algorithm. To avoid the dynamic threat, the time $\mathrm{t}$ of the ant at any node should be known. The ants carry out the following calculation after walking each step to save the current time.

Given a UAV flight speed of $\mathrm{v}$, the route point set that the ant walked is $\left\{{\mathrm{r}}_1,{\mathrm{r}}_2,\cdots ,{\mathrm{r}}_{\mathrm{k}}\right\},\mathrm{k}>1$; when k = 1, the UAV start time is $t_0$. If we know the time ${\mathrm{t}}_{\mathrm{k}-1}$ of the ant at node k − 1, when the ant walks to route point k we look up the distance of ${\mathrm{r}}_{\mathrm{k}-1}{\mathrm{r}}_{\mathrm{k}}$, named ${\mathrm{l}}_{{\mathrm{r}}_{\mathrm{k}-1}{\mathrm{r}}_{\mathrm{k}}}$, and the current time ${\mathrm{t}}_{\mathrm{k}}={\mathrm{t}}_{\mathrm{k}-1}+{\mathrm{l}}_{{\mathrm{r}}_{\mathrm{k}-1}{\mathrm{r}}_{\mathrm{k}}}/\mathrm{v}$. The time of the ant at any node can thus be determined.

The ants excluded the expanding nodes by the following steps:

• When the ant $\mathrm{k}$ at node $\mathrm{i}$ selects the eight expanding nodes, the nodes in the static threat and the ant’s tabu list are first excluded and we then can obtain the remaining expanding nodes ${\mathrm{J}}_{\mathrm{k}}\left(\mathrm{i}\right)$.
• By using the method mentioned above, we calculate the time the ant arrives at every expanding node $\mathrm{j}$, named ${\mathrm{t}}_{\mathrm{j}}$.
• We use the method mentioned in Section 3.3 to implement linear interpolation for every dynamic threat in the dynamic threat set and obtain the threat areas at ${\mathrm{t}}_{\mathrm{j}}$.
• If the node j is located in any threat area, we exclude this node.

The dynamic threat avoidance algorithm will consume many computing resources; therefore, we compile two route planning versions during use. When there is no dynamic threat in the planning scenario, the planning algorithm is used, named MACO, and when we need to consider the dynamic threat, we introduce the dynamic threat avoidance algorithm and name the planning algorithm MACOD.

4.4. Route Optimization

The route planned by the ACO algorithms above is composed of a series of broken line segments. To determine the feasible route, we need to optimize the route by a series of operations.

4.4.1. Route Compression

To avoid frequent turning of the UAV, it is necessary to compress the waypoints, as shown in Figure 10. When the broken line is compressed, it will be shorter than the original path planned by the grid search algorithm. By compressing the route, errors are introduced; however, the route is optimized to a certain extent.

When compressing the waypoints, the adjusted route should not cross the threatened area, and the average altitude of the route has not been significantly improved. After deleting waypoint ${\mathrm{r}}_{\mathrm{i}}$, the route should meet the criteria: (a) the adjusted route segments $\{{\mathrm{r}}_{\mathrm{i}-1},{\mathrm{r}}_{\mathrm{i}+1}\}$ should not be located in the threatened area and (b) the difference value of average altitude between ${\{\mathrm{r}}_{\mathrm{i}-1}{,\mathrm{r}}_{\mathrm{i}+1}\}$ and ${\{\mathrm{r}}_{\mathrm{i}-1}{,\mathrm{r}}_{\mathrm{i}}{,\mathrm{r}}_{\mathrm{i}+1}\}$should not exceed the specified threshold $\pounds$. The average altitude is calculated as follows: interpolate for the route with interval $\forall$ in ${\mathrm{Ʊ}}^3$, then transform the interval points to ${\mathrm{\Psi }}^3$ and obtain the terrain altitude for them.

This paper improves the Douglas-Peucker algorithm [36]; the improved algorithm steps are as follows:

• First, the algorithm compares the distance between each waypoint and the link line of two endpoints. If the distance is less than the tolerance $\mathsf{ϵ}$, judge the route segment to determine if it meets the conditions (a) and (b) after removing this waypoint. If satisfied, then delete this waypoint and if not, then keep this one.
• Select the point that has the farthest distance with the link line of two endpoints as the separation waypoint such that the route is divided into two segments, and then perform Step 1 recursively for these two segments until there is no waypoint to be deleted.

4.4.2. Smooth Route Turning

After compression, the route still consists of broken lines. The UAV cannot smoothly complete the turning process. There are some common segment fitting methods, such as the Dubins curve [15] and the Bézier curve [2], and some of the authors do not smooth the broken line [12]. In fact, a Fly-by Fix or Fly-over Fix [22] is usually used in PBN flight procedures for turning. In this study, we use the Fly-by Fix for UAV smooth turning. In Figure 11a, $\mathrm{DTA}$ is the length of the minimum straight route segment required by the Fly-by Fix; it is determined by the flight speed, turning angle and other parameters and is calculated as follows:

$$\{\begin{array}{l}\mathrm{R}=\frac{{(\min [{500,\mathrm{V}}_{\mathrm{KTAS}}+{\mathrm{V}}_{\mathrm{KTW}}])}^2}{\tan \left(\mathsf{\theta }\right)\times 68625.4}\\ \mathrm{DTA}=\mathrm{R}\times \tan \frac{\mathsf{\alpha }}{2}\end{array}$$

(16)

where $\mathrm{R}$ is the turning radius, ${\mathrm{V}}_{\mathrm{KTAS}}$ is the true airspeed, ${\mathrm{V}}_{\mathrm{KTW}}$ is the tailwind, and $\mathsf{\theta }$ is the slope angle.

UAV flight tends to avoid frequent rising or descending. If a UAV needs to fly from the waypoint $\mathrm{A}$ to waypoint $\mathrm{B}$ with different elevations, the UAV should rise/descend with maximum lifting angle to the horizontal line of $\mathrm{B}$ and then maintain level flight to $\mathrm{B}$.

4.4.3. Route Buffer Area

The route buffer area provides a vertical and horizontal buffer for the route for UAV safety. For the shadow regions shown in Figure 11b, the route buffer area regards the route as the centre line, with a width of 4 × RNP, and the height is a rectangular buffer needed to override the obstacle (ROC). We do interpolation for the entire buffer. If the altitude difference between the interpolation point and the terrain is less than ROC, the height of the waypoint needs to be increased so that the route buffer area does not intersect the terrain and the threat zone.

5. Results and Discussion

In this chapter, we first introduce the platform and parameters in Section 5.1. Pilot experiments with the route planner are reported in Section 5.2, Section 5.3 and Section 5.4. In Section 5.5, we compare six types of planners. To compare the different algorithms, we do not introduce dynamic threats and dynamic threat avoidance algorithms from Section 5.1 to Section 5.5. The planner in this paper is currently named MACO-pl. In Section 5.6, we introduce the dynamic threat and evaluate the performance of MACOD-pl.

5.1. Platform and Parameters

On the basis of the above methods, we developed a route planning system based on a virtual globe platform. This system provides interactive functions such as parameter input, route generation, route editing, route evaluation, threat modelling, data management, flight simulation and so on. Some interactive screenshots can be found in Appendix A. All of the route planning experiments that follow were completed on our route planning platform. The hardware environment was as follows: CPU: Core i7-4790; memory: 8GB DDR3 1600 MHz. The flight parameters of the UAV are shown in

Table 1. The UAV flight parameters.

${\mathbf{V}}_{\mathbf{KTAS}}$ (km/h)	${\mathbf{V}}_{\mathbf{KTW}}$ (km/h)	Maximum Lifting Angle	$\mathbf{h},\mathbf{ROC}$ (m)	RNP (nmi)	Turning Fix
350	0	20°	400	0.25	Fly-by Fix

.

The route buffer area width was 2 × RNP = 0.5 nmi, and ${\mathrm{h}}_{\max }$ and ${\mathrm{h}}_{\min }$ were 5000 and 0, respectively. The algorithm processing time was the average time of the three experiments. The implementation of the algorithm used different programming skills, programming languages and compilers, all of which have large impacts on performance; as a result, the processing time can only be used as a reference.

Much literature research and discussion was reviewed for the selection of parameters in the ACO algorithm [37]. The selection of parameters was based on a large number of experiments. For a specific application, the ACO algorithm requires many experiments to determine better parameters. Through experiments, we found that the processing time of the ACO algorithm was inversely proportional to the evaporation coefficient. With the increase of the evaporation coefficient, the positive feedback effect of pheromones was increased. However, the increase of the evaporation coefficient may cause the algorithm to sink into a local optimum. If the number of ants is too small, the algorithm may terminate prematurely; however, too many ants can increase the time of algorithm iterations. In this study, the parameters of the pheromone were: retention coefficient: $\mathsf{\rho }=0.7$; weight coefficient: $\mathsf{\alpha }=1$, $\mathsf{\beta }=4$; ant number: 30; iteration number: 200; node number n between the start point and end point; threshold $\mathsf{\delta }=10\mathrm{m}$; local uplift threshold in route compression: $\pounds =10\mathrm{m}$; and interpolation interval: $\forall =100\mathrm{m}$.

This experiment used five coordinate groups, as shown in

Table 2. Experimental coordinate groups.

Serial Number	Start Point	End Point	${\mathbf{L}}_{\mathbf{min}}$	${\mathbf{H}}_{\mathbf{avT}}$
	$(\mathbf{x},\mathbf{y},\mathbf{H})$	$(\mathbf{x},\mathbf{y},\mathbf{H})$	(km)	(m)
1	99.96463°	99.43347°	110.613	2339.93
	26.88202°	27.75616°
	1819.58 m	2530.68 m
2	90.86051°	94.46516°	350.051	4492.95
	29.30996°	29.42218°
	3573.99 m	2951.15 m
3	121.84255°	119.05692°	342.512	4.49
	30.98066°	32.95193°
	21.51 m	14.17 m
4	99.44990°	118.45611°	1909.435	1186.37
	27.81014°	25.09058°
	1723.22 m	397.35 m
5	109.78029°	120.41494°	2259.830	203.71
	18.85893°	37.01404°
	445.07 m	132.31 m

.

The planning space needs to be generated in advance, and the altitude property of the node can be stored off-line and can be used based on need. In this study, we generated a multi-granularity planning space; the pre-generated longitude range was [90, 125] and the pre-generated latitude range was [15, 45], corresponding to more than 10,000,000 square kilometres of planning area. At a grid width of 0.01°, the number of nodes was 10,500,000, and the grid space required approximately 5 min to be generated in the experimental machine. Before the implementation of the planning algorithm, we could load the corresponding offline block based on the size of the planning space.

Figure 12 shows the planning area of the first experimental coordinate group; the area is surrounded by the points (100.284, 27.622), (100.047, 26.723), (99.093, 26.920), (99.295, 27.812). Points $\mathrm{A}$ and $\mathrm{B}$ were the start and end points, respectively.

5.2. Get the Initial Route

The parameter $\mathsf{\epsilon }$ balances the valley-following ability and the length of the route. Different experimental results can be obtained by setting different coefficients. The initial altitude of the route is the sum of the node elevation and the ROC value. The experimental results are shown in Figure 13 and

Table 3. Comparison results for different coefficients.

Coefficients ($\mathit{\epsilon }$)	Planning Time (s)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	Waypoint Number	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
0	1.01	118.428	2488.07	86	+7.07%	+6.33%
0.5	1.09	119.074	2351.98	87	+7.65%	+0.51%
1	1.20	119.924	2152.15	88	+8.42%	−8.03%
1.5	1.29	119.995	2138.84	89	+8.48%	−8.59%
2	1.37	120.617	2058.38	90	+9.04%	−12.03%
2.5	1.55	123.085	2028.94	92	+11.28%	−13.29%
3	N/A

. It can be concluded that when $\epsilon$ was relatively high, the valley-following performance of the algorithm was better. However, a higher $\epsilon$ will increase the difficulty for an ant to reach the end point, which also increases the running time of the algorithm. When $\epsilon$ was 0, the algorithm calculated the shortest path; however, the path was longer than the shortest path in the continuous world because of the impact of irregular topography and the grid space. When $\epsilon \ge 3$, a feasible solution may not be obtained. Suitable values of $\mathsf{\epsilon }$ can be selected according to the application requirements. For the following experiments we chose $\epsilon =2$.

Figure 14 shows the convergence results of the MACO algorithm when $\epsilon =2$.

5.3. Route Optimization

5.3.1. Route Compression

The results of Figure 13c were first optimized by compression. As shown in Figure 15 and

Table 4. Comparison results for different tolerances.

$\mathsf{ϵ}$ (nmi)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	Waypoint Number	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
0.25	120.129	2063.33	34	+8.60%	−11.82%
0.50	117.172	2060.58	16	+5.93%	−11.94%
0.75	116.003	2141.67	9	+4.87%	−8.47%
1.00	115.061	2120.13	6	+4.02%	−9.39%

, different results were obtained by setting a different tolerance $\mathsf{ϵ}$.

Table 4 shows that the length of the compressed route may have been shorter than the shortest path obtained by the ACO algorithm in Table 3.

There are both advantages and disadvantages in the influence of tolerance, and the tolerance should be selected based on the actual needs.

5.3.2. Smooth Route Turning

As Figure 16 shows, we used Equation (16) and Fly-by Fix to address all UAV turning corners. With respect to the results of Figure 16, the route length after adjustment was 115.495 km with an average terrain height of 2142.80 m. The turning curve is a set of route control points with an interval of 50 m.

5.3.3. Route Buffer Areas

Figure 17a is the result of the introduction of route buffer zones. A route segment that does not meet the flight safety requirements is marked in red. Figure 17b is the result after height adjustment. After the introduction of the route buffer zone, the route length was adjusted to 115.609 km with an average height of 3228.90 m, and the terrain average height was 2142.50 m. Figure 18 is the vertical figure of route and terrain. The Terrain Altitude-curve is the terrain altitude that the planning route passes through, and the 2RNP Terrain Altitude-curve is the highest terrain altitude within the 2RNP range. The initial route was unsafe; we obtained a safe route that could meet the safety requirements after adjusting the elevation.

After route optimization, ${\mathrm{L}}_{\mathrm{US}}$ = 0; however, the average height of the route changed and the average terrain height passed through by the UAV also changed. Errors introduced by the optimization process may cause a route to pass through a threat area, as was mentioned in Section 3.4.1, and the threat area should thus be expanded when it is mapped into the planning space. The expansion width of the threat zone should be greater than the navigation performance error (0.5 nmi); the width of the corresponding expansion in the 0.01° grid width planning space was 1 node.

5.4. Multi-Granularity Results

In Section 3.4.2, we introduced a multi-granularity planning space suitable for different applications. Based on different grids, the experimental results are shown in

Table 5. Comparison of results with different granularity.

Grid Width	$\mathsf{ϵ}$ (nmi)	Processing Time (s)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
0.007°	0.75	2.15	115.117	1975.70	+4.07%	−15.57%
0.007°	0.50	2.15	116.899	1947.84	+5.68%	−16.76%
0.006°	0.75	2.73	114.404	1979.44	+3.43%	−15.41%
0.006°	0.50	2.73	117.060	1953.05	+5.83%	−16.53%
0.005°	0.75	3.21	114.086	1942.91	+3.14%	−16.97%
0.005°	0.50	3.21	115.603	1905.73	+4.51%	−18.56%
0.004°	0.75	4.24	114.403	1955.89	+3.43%	−16.41%
0.004°	0.50	4.24	115.523	1923.16	+4.44%	−17.81%
0.003°	0.75	6.18	114.498	1961.64	+3.51%	−16.17%
0.003°	0.50	6.18	115.300	1925.10	+4.24%	−17.73%
0.002°	0.75	7.37	115.033	1973.27	+4.00%	−15.67%
0.002°	0.50	7.37	115.711	1961.85	+4.61%	−16.16%

.

As the grid width was reduced and the planning space became closer to the real world, the route we obtained was more consistent with reality. The algorithm ran slower; however, when the grid width was reduced to a certain value limited by factors such as the input parameters and the experimental terrain data source, the results no longer improved. In the experiment scene shown in Figure 19, a better effect was achieved when grid width was 0.005° and $\mathsf{ϵ}$ was 0.50 nmi. The following experiments are based on these parameters.

5.5. Comparison Experiments

On our virtual globe platform we also implemented the route planner of the basic ACO pheromone model named ACO-pl and the Rank-based pheromone model named RAS-pl. The number of retaining ants of RAS-pl was 5. We used the local valley-following method described in Section 4.2.2 and the shortest distance from the target as the heuristic function to realize an A* based planner, named A*-pl.

Additionally, we implemented the planners based on PSO and EA, named PSO-pl and EA-pl, respectively. Because these two methods use the dimension-reduction encode method to achieve evolutionary computation, the local valley-following method in Section 4.2.2 was not applicable. This study used the minimal flight altitude method in [12] to update the average altitude of the route. Because of the interval difference between adjacent waypoints, interval interpolation for the route was needed when calculating the average altitude. In the following experiments, the inertial weight coefficient of PSO-pl was linear from 0.9 to 0.4. The iteration count of PSO-pl and EA-pl was 200 and the population size was 30.

When implementing the six types of planners, this study used data structures to store pre-calculated distances between the adjacent nodes and pre-calculated distances between each node and the end node. The adjacent waypoints in MACO-pl, ACO-pl, RAS-pl and A*-pl are adjacent nodes in the grid planning space, as shown in Figure 20a. The algorithm proposed and used in this study achieved computational acceleration. While the adjacent waypoints in PSO-pl and EA-pl are not necessarily adjacent nodes, as shown in Figure 20b, the waypoints P and Q can only move in two directions. At this point, P and Q are not adjacent nodes and thus distances must be calculated. Although the waypoints Q and H are adjacent nodes, we can obtain the distance between Q and H by a table lookup.

The comparison experiments used the same grid size and UAV parameters and were conducted for different coordinate groups. The results are shown in

Table 6. Planning results of 1–5 coordinate groups.

Serial Number	Method	Processing Time (s)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
1	MACO-pl	3.21	115.603	1905.73	+4.51%	−18.56%
1	ACO-pl	10.07	118.831	2262.32	+7.43%	−3.32%
1	RAS-pl	5.66	117.267	1914.28	+6.02%	−18.19%
1	A*-pl	2.16	116.875	1916.29	+5.66%	−18.10%
1	PSO-pl	9.78	118.679	2177.42	+7.29%	−6.94%
1	EA-pl	10.84	119.441	2190.32	+7.98%	−6.39%
2	MACO-pl	12.39	398.782	3471.04	+13.92%	−22.74%
2	ACO-pl	38.48	389.907	4074.93	+11.39%	−9.30%
2	RAS-pl	20.40	399.823	3623.61	+14.22%	−19.35%
2	A*-pl	148.14	396.934	3681.73	+13.39%	−18.05%
2	PSO-pl	40.02	386.679	4081.56	+10.46%	−9.15%
2	EA-pl	41.28	385.041	4094.39	+10.00%	−8.87%
3	MACO-pl	10.96	362.433	2.85	+5.82%	−36.53%
3	ACO-pl	36.42	368.345	2.98	+7.54%	−33.63%
3	RAS-pl	18.67	365.184	2.93	+6.62%	−34.74%
3	A*-pl	116.23	363.245	2.94	+6.05%	−34.52%
3	PSO-pl	39.57	366.176	4.01	+6.91%	−10.69%
3	EA-pl	38.95	364.504	3.95	+6.42%	−12.03%
4	MACO-pl	80.12	2115.333	1043.90	+10.78%	−12.01%
4	ACO-pl	201.15	2135.830	1097.04	+11.85%	−7.53%
4	RAS-pl	149.98	2136.106	1063.77	+11.87%	−10.33%
4	A*-pl	N/A
4	PSO-pl	278.26	2104.857	1137.25	+10.23%	−4.14%
4	EA-pl	299.38	2109.891	1131.32	+10.50%	−4.64%
5	MACO-pl	78.84	2371.012	84.59	+4.92%	−58.48%
5	ACO-pl	197.43	2401.023	116.56	+6.25%	−42.78%
5	RAS-pl	136.01	2390.174	85.99	+5.76%	−57.79%
5	A*-pl	N/A
5	PSO-pl	301.44	2398.322	127.76	+6.13%	−37.28%
5	EA-pl	279.64	2388.459	134.31	+5.69%	−34.07%

and Figure 21.

The pheromone evaporation model of ACO-pl easily resulted in the local route planning optimum. Compared with ACO-pl, RAS-pl used the results of higher ranking ants and had better exploration ability. The A* algorithm is highly efficient in solving a simple shortest path problem; however, the performance of this algorithm was reduced when the local terrain was introduced. A*-pl can have good solutions in small scenes; however, its open table became increasingly large when the scene became larger and more complex. Although we introduced a binary heap to speed up the query, it remained very difficult to find a solution. The valley-following abilities of EA-pl and PSO-pl were not as good; however, the performance for the height and length of these methods was relatively balanced. Moreover, EA-pl and PSO-pl were less efficient than MACO-pl because of frequent interpolation and distance calculation. Compared with the A* algorithm, intelligent planners have advantages in solving large-scale complex problems. We demonstrated by the comparisons that MACO-pl produced a better solution and that the algorithm was more efficient.

Figure 22 shows the results of MACO-pl. In the area of Figure 22a, the planner could track the trend of the terrain in the valley, which reflects better terrain tracking ability. Figure 22b is mostly a plains area, and the trend of the route was relatively straight. Figure 22c,d show the planning results in large-scale planning scenarios. Figure 23 shows some local details of these routes.

After marking the radar zone, no-fly zone, and missile threat zone, we did experiments on coordinate groups 1, 2, 4, and 5.

Table 7. Planning results in threat environments.

Serial Number	Method	Processing Time (s)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	${\mathbf{L}}_{\mathbf{FA}}$ (km)	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
1	MACO-pl	5.63	136.480	2725.64	0	+23.39%	+16.48%
1	ACO-pl	17.02	138.737	2870.48	0	+25.43%	+22.67%
1	RAS-pl	9.39	137.158	2757.80	0	+24.00%	+17.86%
1	A*-pl	5.27	136.497	2789.72	0	+23.40%	+19.22%
1	PSO-pl	29.84	138.528	2955.42	0	+25.24%	+26.30%
1	EA-pl	31.15	137.917	2918.81	0	+24.68%	+24.74%
2	MACO-pl	20.85	421.922	4379.01	0	+20.53%	−2.54%
2	ACO-pl	75.35	431.943	4687.66	0	+23.39%	+4.33%
2	RAS-pl	37.28	423.871	4470.49	0	+21.09%	−0.50%
2	A*-pl	326.92	422.775	4450.11	0	+20.78%	−0.95%
2	PSO-pl	91.49	430.563	4716.47	0	+23.00%	+4.97%
2	EA-pl	99.02	428.719	4671.99	0	+22.47%	+3.98%
4	MACO-pl	127.39	2250.836	1086.37	0	+17.88%	−8.43%
4	ACO-pl	391.50	2301.760	1143.74	0	+20.55%	−3.59%
4	RAS-pl	271.98	2284.541	1109.44	0	+19.64%	−6.49%
4	A*-pl	N/A
4	PSO-pl	401.32	2304.283	1217.07	0	+20.68%	+2.59%
4	EA-pl	384.56	2300.575	1203.64	0	+20.48%	+1.46%
5	MACO-pl	84.26	2438.912	86.22	0	+7.92%	−57.68%
5	ACO-pl	209.92	2494.026	119.14	0	+10.36%	−41.51%
5	RAS-pl	140.15	2477.678	86.85	0	+9.64%	−57.37%
5	A*-pl	N/A
5	PSO-pl	329.48	2489.205	129.49	0	+10.15%	−36.43%
5	EA-pl	286.77	2503.418	136.61	0	+10.78%	−32.94%

and Figure 24 show the results of different planners. We demonstrated that the MACO-pl had a better solution. Figure 25 shows local results by MACO-pl of four coordinate groups after plotting threatened areas. The planned routes could make use of the terrain condition to avoid the threatened area.

5.6. Experiments on Dynamic MACO

This section introduces the dynamic threat area discussed in Section 3.3. Based on MACOD-pl, we used coordinate group 1 to test the algorithm. Considering that the route optimization algorithms would bring approximately 1%–3% error to the initial route length, which is approximately 1–4 km under the current scenario, the maximum time error was approximately 40 s, and the moving speed of the dynamic threat area in the figure was approximately 12.6 m/s. The maximum error range was thus approximately 510 m. To ensure route safety, the boundary of the dynamic threat area was expanded by 0.01°. As shown in Figure 26 and

Table 8. The results of different start times under dynamic threats.

Start Time (s)	Planning Time (s)	${\mathbf{L}}_{\mathbf{route}}$ (km)	${\mathbf{h}}_{\mathbf{avT}}$ (m)	Compared with ${\mathbf{L}}_{\mathbf{min}}$	Compared with ${\mathbf{H}}_{\mathbf{avT}}$
0	10.01	133.850	2112.65	+21.01%	−9.71%
1000	12.14	135.633	2724.18	+22.62%	+16.42%
2000	12.89	136.204	2725.49	+23.14%	+16.48%
3000	9.49	115.377	1906.71	+4.31%	−18.51%

, when the start time of the UAV was not the same, the UAV could meet the threat at different times and places, and the planned routes were thus different. When the start time of the UAV was between 0 and 2000 s, the dynamic threat area would block the UAV from flying in the valley area. When the start time of the UAV was 3000 s, the dynamic threat area would pass through the valley area and the result of route planning was similar to the result without the threat. Because of the introduction of dynamic threats, the dynamic threat avoidance algorithm needed to consume many computing resources, and the operating efficiency of the MACOD-pl algorithm was much lower than that of the MACO-pl algorithm.

6. Conclusions and Future Work

This paper presented an improved ACO-based planner based on the virtual globe platform to solve the route planning problem in risk environments. We developed a route planning system and realized six different planners in the static threat. By experimental analysis, we demonstrated that our optimum planner had better performance in terms of fuel consumption, terrain masking, and risk avoidance. Finally, we demonstrated the effectiveness of MACOD-pl in the dynamic threat environment. Benefiting from the virtual globe platform, the method presented in this paper can achieve route planning globally and meet the current maximum combat radius of UAVs in battlefields. The planned route can effectively avoid various threats and, by changing the parameters, the planner is suitable for military penetration, search-and-rescue, and many other scenarios. Our route planning system can perform not only automatic planning but also route editing and storage, flight simulation experiments, and other tasks.

Many other types of planning algorithms exist; for example, some algorithms use a rotated coordinate frame to reduce the dimension of the problem [12,18]. Further research on how to apply other algorithms to the virtual globe platform and efficiently obtain better results would be worthwhile, as well as further research on the model of different dynamic threats and the efficient use of dynamic threat avoidance algorithms.