International Journal of

Geo-Information Article

Unmanned Aerial Vehicle Route Planning in the Presence of a Threat Environment Based on a Virtual Globe Platform Ming Zhang 1 , Chen Su 1 , Yuan Liu 2 , Mingyuan Hu 2 and Yuesheng Zhu 1, * 1 2

*

Shenzhen Graduate School, Peking University, Shenzhen 518055, China; [email protected] (M.Z.); [email protected] (C.S.) Institute of Space and Earth Information Science, The Chinese University of Hong Kong, Hong Kong; [email protected] (Y.L.); [email protected] (M.H.) Correspondence: [email protected]; Tel.: +86-755-2603-5352

Academic Editor: Wolfgang Kainz Received: 5 July 2016; Accepted: 28 September 2016; Published: 10 October 2016

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. Keywords: unmanned aerial vehicle; route planning; virtual globe; ant colony optimization; 3D environments

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–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–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.

ISPRS Int. J. Geo-Inf. 2016, 5, 184; doi:10.3390/ijgi5100184

www.mdpi.com/journal/ijgi

ISPRS Int. J. Geo-Inf. 2016, 5, 184

2 of 30

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.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

3 of 30

ISPRS Int. J.Planning Geo-Inf. 2016, 5, 184 2. Route Problem Descriptions

3 of 30

3 The P in space Ψ3 isΨdescribed by longitude x, latitude y, and altitude The location locationofofa apoint point P geodetic in geodetic space is described by longitude x, latitude y, and 3 H. A UAV flies along the designated route point sequence R = r , r , · · · , r , r ∈ Ψ , { } route 1 R 2 route = n{r1 ,r2i,⋯,rn }, riwhere altitude H . A UAV flies along the designated route point sequence ∈ Ψ3 , rwhere point and rn is the point. The goal of The penetration is to obtain feasible routea 1 is therstart andend rn is the end point. goal of routing penetration routinga is to obtain 1 is the start point with minimal fuel consumption, maximum terrain masking, and minimum risk. feasible route with minimal fuel consumption, maximum terrain masking, and minimum risk. Both the navigation navigationsystems systemsand andthe the control systems of UAVs contain certain deviations. Both the control systems of UAVs contain certain deviations. The The international aviation organization proposed conceptofof the the required required navigation international civilcivil aviation organization proposed thetheconcept navigation performance (RNP). RNP RNP [22] [22]describes describesthe theprecision precisionthat thatcan can attained during at least of performance (RNP). bebe attained during at least 95%95% of the the flight time; the precision unit is the nautical mile (nmi). In Figure 1, the segment width of the route flight time; the precision unit is the nautical mile (nmi). In Figure 1, the segment width of the route is is defined × RNP. The planned route may not satisfy performance constraints UAVs. defined asas 4 ×4 RNP. The planned route may not satisfy thethe performance constraints ofof UAVs.

Figure 1. Plan view for RNP segment width. Figure 1. Plan view for RNP segment width.

We evaluate the planning route on four dimensions: minimal fuel consumption, maximum Wemasking, evaluateminimum the planning on four dimensions: minimal fuel maximum terrain risk route and performance safety. Additionally, we consumption, evaluate the efficiency of terrain masking, minimum risk and performance safety. Additionally, we evaluate the efficiency of the planner: the planner:  Lroute : Route length. Connect the points in Rroute successively and obtain a series of line length. Connect the points successively and obtain a series of line • Lsegments, route : Route then calculate the total lengthinofRroute the segments to obtain the route length. Tosegments, evaluate then calculate the total length of the segments to obtain the route length. To evaluate the result, the result, we introduce the theoretically shortest length of the route, Lmin . Then, we connect the we introduce the theoretically shortest length of the route, L . Then, we connect the r and r of n min 1 r1 and rn of Rroute to obtain the straight line segment r1 rn . Lmin is the length of r1 rn . R to obtain the straight line segment r1 rn . Lmin is the length of r1 rn .  hroute avT : Average terrain altitude passed by the route. havT is calculated as follows: Obtain the • hinterpolation : Average terrain{raltitude passed by the route. havT is calculated as follows: Obtain the avT points 1 (x1 ,y1 ,H1 ),r2 (x2 ,y2 ,H2 ),⋯,rk (xk ,yk ,Hk )} with interval ∀ of Rroute . We can interpolation points {r1 (x 1 , y1 , H1 ), r2 (x 2 , y2 , H2 ) , · · · , rk (x k , yk , Hk )} with interval ∀ of Rroute . then obtain points on the terrain {r1 (x1 ,y1 ,H1T ),r2 (x2 ,y2 ,H2T ),⋯,rk (xk ,yk ,HkT )} , and We can then obtain points on the terrain {r1 (x 1 , y1 , H1T ), r2 (x 2 , y2 , H2T ) , · · · , rk (x k , yk , HkT )}, havT = ∑ki=1 HiT ⁄k . To evaluate the result, we introduce a reference HavT , the average terrain and havT = ∑ki=1 HiT /k. To evaluate the result, we introduce a reference H , the average altitude passed by r1 rn . We can understand the terrain masking ability ofavT the planner by terrain altitude passed by r1 rn . We can understand the terrain masking ability of the planner by contrasting the result with HavT . contrasting the result with HavT .  LFA : Route length that passes through the threat zones. • LFA : Route length that passes through the threat zones.  LUS : Route length that does not meet the flight performance safety requirements. • LProcessing US : Route length time. that does not meet the flight performance safety requirements. • Processing time. 3. Build the Planning Space 3. Build the Planning Space 3.1. Virtual Globe Space 3.1. Virtual Globe Space For the penetration tasks on the virtual globe, the planning space needs to be converted between For the penetration tasks on the virtual globe, the3planning space needs to be converted between 3 Geodetic space spaceΨ3Ψand and Cartesian space Ʊ3 . Ψ treatsasthe earth as a reference 3 . Ψ3 treats Geodetic Cartesian space U the earth a reference ellipsoid. P (x, y,ellipsoid. H) ∈ Ψ3 3 3 P(x,y,H)∈Ψ can be converted into 𝑄(X, Y, Z) ∈ Ʊ as follows: 3 can be converted into Q (X, Y, Z) ∈ U as follows:  X = (M + H) cos(y) cos(x) (M + H) cos (y) cos (x)   XY= = (M + H) cos(y) sin(x) Y = h(M + H) cos (y2) sin (xi) aa22 − −b2b   ZZ= H sin(y) sin (y) M(11 −− a22 )++ H] = [M a {

(1) (1)

where asas follows: where M Mrepresents representsthe theradius radiusofofcurvature curvatureofofthe thereference referenceellipsoid, ellipsoid,computed computed follows:

= q MM =

aa 2

2

2a − bb2sin2 (y) √1 1−−a − a2 sin2 (y) a2

(2) (2)

where a and b represent the length of the major axis and the minor axis of the reference ellipsoid, respectively.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

4 of 30

ISPRS J. Geo-Inf. 2016, 5, 184 the length of the major axis and the minor axis of the 4reference of 30 where a Int. and b represent ellipsoid, respectively. Q(X, Y, Z)∈Ʊ3 can be converted into P(x,y,H)∈Ψ3 as follows: Q (X, Y, Z) ∈ U3 can be converted into P (x, y, H) ∈ Ψ3 as follows:

          

Y x = arctan( ) Y X( ) x = arctan X

Z(M H)H) Z(M+ +

h   i)) 2 M 1 −2 a2 − b22 + H X + Y ( ) a −a2b 2 2

y == arctan( arctan( r

2

√ √(X + Y ) [M (1 − X2 + Y2 H = cos − M (y)

a2

) + H]

(3)

(3)

√X2 + Y2 H= − M { The reference ellipsoid used bycos(y) the Global Positioning System and virtual globe platforms is the

World Geodetic System of 1984used ellipsoid ellipsoid) [22]. and virtual globe platforms is the The reference ellipsoid by the (WGS-84 Global Positioning System The lengths of System the same longitude interval in ellipsoid) different[22]. latitudes and the same latitude interval World Geodetic of 1984 ellipsoid (WGS-84 3 ◦ , 4000 m) Thelongitudes lengths of the interval different latitudes and the same latitude interval in different aresame not longitude fixed in Ψ . Forinexample, the distance between (90◦ , 29 3 ◦ ◦ ◦ ◦ in different longitudes are not fixed in Ψ . For example, the distance between (90°, 29°, 4000 m) and◦ , 26◦ , and (91 , 29 , 4000 m) is 97,422.04 m, whereas the distance between (90 , 26 , 4000 m) and (91 3 m) and (91°, 26°, 4000 (91°, 29°, 4000 m) is 97,422.04 m, whereas the distance between (90°, 26°, 4000 4000 m) is 100,114.77 m. In the later sections, the planning space is in Ψ , but the distance measure, m) is 100,114.77 m. Inequation the latercalculation sections, theand planning space is in are Ψ3 , carried but the out distance interpolation algorithm, other operations in U3 .measure, interpolation algorithm, equation calculation and other operations are carried out in Ʊ3 . The terrain is usually described by digital elevation model (DEM) data, whereas terrain texture is The terrain is usually described by digital elevation model (DEM) data, whereas terrain texture described by digital orthophoto map (DOM) data. The data of DEM and DOM can be represented is described by digital orthophoto map (DOM) data. The data of DEM and DOM can be represented as {(as x, y, z)i }, },where y is latitude.InInDEM DEM data, is the altitude {(x,y,z) wherex xisislongitude longitude and and y is latitude. data, z isz the altitude valuevalue of theof the i corresponding position, whereas z is the pixel value for DOM data, as Figure 2a,b illustrate. corresponding position, whereas z is the pixel value for DOM data, as Figure 2a,b illustrate.

(a)

(b)

(c)

(d)

Figure 2. (a) Description of regular grid DEM; (b) DOM data; (c) DEM and DOM data visualization;

Figure 2. (a) Description of regular grid DEM; (b) DOM data; (c) DEM and DOM data visualization; (d) Virtual globe platform. (d) Virtual globe platform.

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

5 of 30 5 of 30

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

3−i i ji j ∑ Zp Z= (xpp,,y yp)) ==∑ 3i= xpp yp p =ff(x 0 aij xapijy i=0 j=0 j=0 3

∑

∑

(4) (4)

We We can can use use the the adjacent adjacent 44 ××44points pointsofofpoint pointA A in in Figure Figure 2a 2a to to calculate calculate the the coefficients coefficients of of the the interpolation A. interpolation function function to to determine determine the the altitude altitude value value of of point point A. In 2c, we we can canrealize realizethe thevisualization visualizationofof3D 3D terrain superimposing DOM data on In Figure Figure 2c, terrain byby superimposing thethe DOM data on the the DEM data. In this study, we used the improved NASA-WorldWind platform; the data source was DEM data. In this study, we used the improved NASA-WorldWind platform; the data source was the the SRTM DEM at 90-m resolution and the Landsat DOM 60-m resolution SRTM DEM datadata at 90-m resolution and the Landsat DOM data at data 60-mat resolution providedprovided publicly publicly by NASA, as shown in Figure 2d. by NASA, as shown in Figure 2d. 3.2. 3.2. Threat Threat Modelling Modelling 3.2.1. General General Types Types of Threats 3.2.1. No-fly zones, zones, hazardous hazardous weather, weather, high-rise buildings and low-altitude low-altitude control control zones zones need need to to be be No-fly considered when planning routes. considered No-fly zones a UAV cannot fly into. A no-fly zone can be can described as the enclosed No-fly zonesare areregions regionsthat that a UAV cannot fly into. A no-fly zone be described as the area defined the point set point A = set {A1A , A , ·1·, ·A, 2A ,m > 2, as shown in Figure 3a. m }A enclosed areaby defined by the = 2{A , ⋯, m }, m > 2, as shown in Figure 3a.

(a)

(b)

(c)

Figure 3. (a) 2D No-fly zone; (b) 3D slowly hazardous weather model; (c) 3D threat model of a 500-m Figure 3. (a) 2D No-fly zone; (b) 3D slowly hazardous weather model; (c) 3D threat model of a 500-m high TV high TV tower tower described described as as {(114.655225, {(114.655225, 23.269952), 23.269952), (114.658225, (114.658225,23.269952), 23.269952),(114.658225, (114.658225,23.272952), 23.272952), (114.655225, 23.272952), 500}. (114.655225, 23.272952), 500}.

High buildings, slowly hazardous weather and low-altitude control zones are low-altitude High buildings, slowly hazardous weather and low-altitude control zones are low-altitude regions regions that UAVs cannot fly into. Such threats can be described as {{A1 , A2 , A3 , ⋯, Am }, Hmax }, m > 2, that UAVs cannot fly into. Such threats can be described as {{A , A , A3 , · · · , Am }, Hmax }, m > 2, where Hmax is the upper bound. This type of threat is modelled1 as 2observed in Figure 3b,c. These where Hmax 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 threats can be dealt with as terrain data; we transform this type of threat area into terrain before planning. route planning. 3.2.2. Killing Zone of an Air Defense Missile 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 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 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 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 model of a vertical killing zone is shown in Figure 4a. Horizontal lines define the high and low boundaries with heights of Hmax and Hmin , respectively. Dsymin is the minimum slant range of the boundaries with heights of Hmax and Hmin , respectively. Dsymin is the minimum slant range of the far far boundary of the killing zone and Dsymax is the maximum slant range. The minimum slant range and the maximum height angle of the near boundary of the killing zone are Dsjmin and Amax ,

ISPRS Int. J. Geo-Inf. 2016, 5, 184

6 of 30

boundary of the killing zone and Dsymax is the maximum slant range. The minimum slant range6 of and 30 the maximum height angle of the near boundary of the killing zone are Dsjmin and Amax , respectively. The near boundary far boundary the arcs O aswith the centre. shows 3D respectively. The nearand boundary and far are boundary arewith the arcs O as theFigure centre.4b Figure 4b the shows visualization model of the missile killing zone. the 3D visualization model of the missile killing zone. ISPRS Int. J. Geo-Inf. 2016, 5, 184

(a)

(b)

Figure Figure 4. 4. (a) (a) Model Model of of aa missile missile killing killing zone; zone; (b) (b) 3D 3D visualization visualization of of aa missile missile killing killing zone. zone.

3.2.3. Radar Threat Space 3.2.3. Radar Threat Space Radar plays a very important role in modern air defence systems; it is a key threat to the 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 penetration of UAVs that needs to be avoided by route planning. Affected by the terrain masking and and the earth curvature, it is difficult for radar to find a target in low-altitude flight. 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 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 maximally tolerate. Given the detection probability threshold and radar performance parameters, the maximum range of radar Rmax can be estimated by the radar equation [25]. Assuming that the radar maximum range of radar Rmax can be estimated by the radar equation [25]. Assuming that the radar is deployed at Ra (xa ,ya ,Ha ) in Ψ3 ,3 we transform it into Ʊ33 as Ra (Xa ,Ya ,Za ), with the direction of is deployed at Ra (x a , ya , Ha ) in Ψ , we transform it into U as Ra (X a , Ya , Za ), with the direction of pitch andazimuth azimuthangle angleβ,β,and andwe wecan candescribe describeppk((X k ,Yk ,Zk ) on the boundary of the radar pitch angle angle θθ and k Xk , Yk , Zk ) on the boundary of the radar coverage as: coverage as:   f0 (θ) cos max  Xk X=k =XXa a = + RRmax f' (θ)cos (β) (β) (5) Yk = Ya = Rmax' f0 (θ) sin (β) (5)  + Rmax f (θ)sin (β) 0  Z{Y=k =ZYa = R f θ ( ) a max k Zk = Za + Rmax f' (θ) where f0 (θ) describes the beam pattern of the radar antenna in the direction of the pitch angle. where f' (θ) 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: We use the Gaussian function to approximate the beam pattern as: √

2

4ln 2θ2 /θ/θb f0 (f'θ(θ) ) ==ee−−4ln√2θ b

(6) (6)

where thesignal signalbeam beamwidth. width. where θθbb isisthe The radar threat curved bebe obtained. TheThe probability of radar detection on this The radar threat curvedsurface surfacecan canthus thus obtained. probability of radar detection on surface is constant. Outside the curved surface, the radar detection probability is less than the this surface is constant. Outside the curved surface, the radar detection probability is less than the probability probability threshold threshold and and the the UAV UAV is is regarded regarded as as safe. safe. Within Within the the curved curved surface, surface, the the radar radar detection detection probability is higher than the probability threshold and the UAV is not safe. probability is higher than the probability threshold and the UAV is not safe. In low-altitude areas, areas,radar radarsignals signals may masked by terrain. In Figure 5a, radar are In low-altitude may be be masked by terrain. In Figure 5a, radar beamsbeams are partly partly masked by terrain and the shaded area represents the area that signals can reach. The masked by terrain and the shaded area represents the area that signals can reach. The visualization of visualization of can the radar beam canas befollows: constructed as follows: the radar beam be constructed 1. 1. 2.

Calculate the detection range of the radar with different sampled pitch angles by using Equations Calculate the detection range of the radar with different sampled pitch angles by using (5) and (6); this produces a series of sampling points such as A–G, as shown in Figure 5a. Equations (5) and (6); this produces a series of sampling points such as A–G, as shown in Connect the radar centre O to the sampling point and obtain a series of sampling points on this Figure 5a. straight line by sampling with short intervals. Then, we transform these sampling points to Ψ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 OD, the boundary point D is moved to M.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

2.

7 of 30

Connect the radar centre 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 Ψ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 OD, the boundary point D is moved to M.

ISPRS Int. Int. J.J. Geo-Inf. Geo-Inf. 2016, 2016, 5, 5, 184 184 ISPRS

of 30 30 77 of

(a) (a)

(b) (b)

Figure 5. (a) (a) Radar signal signal propagation; (b) (b) Earth curvature's curvature's influence. influence. Figure Figure 5. 5. (a) Radar Radar signal propagation; propagation; (b) Earth Earth curvature's influence.

If we we calculate calculate the the radar radar coverage coverage directly, directly, the the earth earth curvature curvature will will introduce introduce considerable considerable error. error. If If we calculate the radar coverage directly, the earth curvature will introduce considerable error. As shown shown in in Figure Figure 5b, 5b, O O is is the the earth earth centre, centre, BD BD is is the the horizontal horizontal plane, plane, the the circular circular arc arc CE CE is is the As the As shown in Figure 5b, O is the earth centre, BD is the horizontal plane, the circular arc CE is the geoid, geoid, and point B is in the same horizontal plane as point A. When B is very near to A, the geoid, and point B is in the same horizontal plane as point A. When B is very near to A, the and point B is in the same horizontal plane as point A. When B is very near to A, the altitudes of A altitudes of of A A and and B B are are approximately approximately the the same same and and may may be be ignored. ignored. However, However, when when B B is is far far from from altitudes and B are approximately the between same andAmay when is far from A, the elevation A, the the elevation difference andbeB Bignored. due to to However, the influence influence ofBthe the earth curvature, ∆H in in A, elevation difference between A and due the of earth curvature, ∆H difference between A and B due to the influence of the earth curvature, ∆H in Figure 5b, cannot be Figure 5b, cannot be ignored. L represents the length of AB, R is the average curvature radius of Figure 5b, cannot be ignored. L represents the length of AB, R is the average curvature radius of ignored. represents of AB, the earth, earth,Land and ∆H can canthe belength calculated by:R is the average curvature radius of the earth, and ∆H can be the ∆H be calculated by: calculated by: p 2 2− R √R (7) ∆H R22 +++L L22L− (7) −R R (7) ∆H∆H === √ When the the distance distancebetween betweentwo twopoints pointsisis is100 100km, km,the theelevation elevation difference is approximately approximately 785 the distance between two points 100 km, the elevation difference When difference is is approximately 785785 m 3 3 3 3 3 3 m because because of influence the influence influence ofcurvature. the curvature. curvature. First, wetoneed need to convert convert from Ψ to tocalculating Ʊ when when because of the of theof First, we need convert from Ψ from to U when m of the the First, we to Ψ Ʊ calculating the boundary of the the radar radar coverage. Then, we we can calculate theand radar coverage and the boundary of boundary the radar coverage. Then, we can calculate thecan radar coverage perform visibility calculating the of coverage. Then, calculate the radar coverage and 33 and the errors caused by the 3 perform visibility analysis. Finally, the results are converted back to Ψ analysis. Finally, the resultsFinally, are converted back Ψ and the errors by the curvature of the perform visibility analysis. the results aretoconverted back to Ψcaused and the errors caused by curvature of the eliminated. earth have have been been eliminated. eliminated. earth haveof been curvature the earth Figure 66 shows shows the the 3D 3D visualization visualization of of the the radar radar network network coverage coverage under under the the influence influence of of the the Figure earth’s curvature. curvature. earth’s

Figure 6. 6. Visualization Visualization of of the the radar radar networking networking coverage coverage on on the the virtual virtual globe. globe. Figure Figure 6. Visualization of the radar networking coverage on the virtual globe.

3.3. Dynamic Dynamic Threat Threat 3.3. The models models presented presented in in Section Section 3.2 3.2 can can simulate simulate most most of of the the threats threats in in normal normal circumstances. circumstances. The Sometimes the threat movement needs to be considered, for example, the forecast of hazardous hazardous Sometimes the threat movement needs to be considered, for example, the forecast of weather such as a typhoon (hurricane) and other predictable dynamic threat areas. At this time, we we weather such as a typhoon (hurricane) and other predictable dynamic threat areas. At this time, need to to introduce introduce the the dynamic dynamic threat threat model. model. Based Based on on the the general general threat threat model, model, aa time time stamp stamp tt is is need introduced, the property of a dynamic threat at t is described by threat = {{A , A , A ,⋯, A }, m}, introduced, the property of a dynamic threat at t is described by threattt = {{A11, A22, A33,⋯, Am Hmax ,t}, m m >> 2. 2. A A full full path path dynamic dynamic threat threat can can be be described described by by {threat {threattt1,, threat threattt2,, ⋯,threat ⋯,threatttn}, }, where where max,t}, H 1

2

n

ISPRS Int. J. Geo-Inf. 2016, 5, 184

8 of 30

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 t is introduced, the property of a dynamic threat at t is described by threatt = {{A1 , A2 , A3 , · · · , Am }, Hmax , t},  m > 2. A full path dynamic threat can be described by threatt1 , threatt2 , · · · , threattn , where the threat begins at t1 and disappears at tn . As shown in Figure 7a,b, if we want to know the state of a threat at any time tk , we can calculate it as follows: Search the interval [ti , ti+1 ] that satisfies ti ≤ tk < ti+1 to find the threatti and threatti+1 . If not then2016, there is no threat area at this time. ISPRS found, Int. J. Geo-Inf. 5, 184 8 of 30 2. Given the threat tk = 2.4 shown in Figure 7b, use linear interpolation for threatti and threatti+1 2. Given thethreat threatt .tkAs = 2.4 shown in Figure 7b,can useuse linear interpolation forfor threat andvertices threatti of to obtain shown in Figure 7b, we linear interpolation the tfour i +1 k to obtain threat shown in Figure 7b, we can use linear interpolation for the four vertices of threat threat to obtain threat . tk . tAs ti and t i+1 k threatti and threatti + 1 to obtain threattk . 1.

(a)

(b) Figure threat model. (a) There is a dynamic threatthreat described as {{(99.94497, 27.60066), (99.859, Figure 7. 7.Dynamic Dynamic threat model. (a) There is a dynamic described as {{(99.94497, 27.60066), 27.38197), (99.68604, 27.47941), (99.77718, 27.68262), 3000, 0 s}, {(99.54497, 27.40066), (99.459, 27.18197), (99.859, 27.38197), (99.68604, 27.47941), (99.77718, 27.68262), 3000, 0 s}, {(99.54497, 27.40066), (99.459, (99.28604, 27.48262), (99.54497, 27.40066), 3000,27.40066), 3600 s}, {(99.34497, 27.00066), (99.259, 27.18197),27.27941), (99.28604,(99.37718, 27.27941), (99.37718, 27.48262), (99.54497, 3000, 3600 s}, {(99.34497, 26.78197), (99.08604, 26.87941), (99.17718, 27.08262), (99.34497, 27.00066), 3000, 7200 s}}; (b) Perform 27.00066), (99.259, 26.78197), (99.08604, 26.87941), (99.17718, 27.08262), (99.34497, 27.00066), 3000, linear interpolation the dynamic threatfor to the obtain the threat status at anythe time. 7200 s}}; (b) Performfor linear interpolation dynamic threat to obtain threat status at any time.

3.4. 3.4. Build Build the the Grid Grid Planning Planning Space Space Route planning on the Route planning on the virtual virtual globe globe platform platform in in this this study study is is based based on on aa graph graphsearch searchalgorithm. algorithm. We define define aa planning planning graph graph composed composed of of nodes nodes and and edges. edges. The We The graph graph needs needs to to effectively effectively express express threats terrain features. features. threats and and terrain In low-altitude penetration penetration missions, models need need to to be be marked. marked. However, In low-altitude missions, threat threat models However, real-time real-time sampling of a threat body at any point in the planning space will dramatically increase the delay time sampling of a threat body at any point in the planning space will dramatically increase the delay of the planner. This paper introduces a type of grid space with location information, elevation time of the planner. This paper introduces a type of grid space with location information, elevation information, attribute information. information. information, and and other other attribute 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 h by default in this space. Each node stores the information including [lon, lat, alt, H ], where lon and lat represent the longitude and latitude of the node,

ISPRS Int. J. Geo-Inf. 2016, 5, 184

9 of 30

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 h by default in this space. Each node stores the information including [lon, lat, alt, Hmin ], where lon and lat represent the longitude and latitude of the node, respectively, alt represents the terrain altitude superposed by the high buildings and low-altitude control zones, and Hmin is the lowest bound of the radar coverage area or the missile killing zone in this area. If the node h that needs to be extended is higher than Hmin , the node will be eliminated. As the figure shows, nodes A, B, C and D are in a radar coverage area. However, because of the terrain masking effect, points B and C can be used as planning candidate nodes with higher Hmin , whereas nodes A and D with low H will be directly eliminated. As the figure shows, node E, F, G and H ISPRS Int. J. Geo-Inf. 2016, 5, 184 min 9 of 30 areISPRS in a Int. missile threat area and will be eliminated because of their low H . If the current point J. Geo-Inf. 2016, 5, 184 9 is of 30 min located in the no-fly zone, H is set to 0 and the point will be eliminated by the planner. Because of min planner. Because of the planning space generated beforehand, the planner can quickly eliminate some planner. Because of the planning space generated beforehand, theeliminate planner can quickly eliminate some the planning space generated beforehand, the planner can quickly some infeasible nodes. infeasible nodes. infeasible nodes.

Figure 8. Space division. Figure 8. Space division. Figure 8. Space division.

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

(a) (a)

(b) (b)

Figure 9. (a) Route planning with threat; (b) Route planning with threat area expansion. Figure 9. (a) Route planning with threat; Route planning with threat area expansion. Figure 9. (a) Route planning with threat; (b)(b) Route planning with threat area expansion.

3.4.2. Multi-Granularity Planning Space 3.4.2. Multi-Granularity Planning Space 3.4.2.The Multi-Granularity Planning Space finer granularity of the grid space is because the more nodes in planning space, the more The finer granularity of the grid space is because the more nodes in planning space,balance the more precise the planning results will be, although planning will take more time. To obtain The finer granularity of the grid space is the because the more nodes in planning space, athe more preciseperformance the planningand results will be, the planning scenarios will take more time. To obtain a balance between precision foralthough different andtime. navigation performance precise the planning results will be, although the application planning will take more To obtain a balance between performance and precision forneed different scenarios and navigation performance parameters of a UAV, different grid sizes to application beapplication used. scenarios between performance and precision for different and navigation performance parameters of a UAV, different grid sizes need to be used. We need to grid construct different sizes of the grid using DEM data. Common parameters of a resampling UAV, different sizes need to be used. We need resampling to construct different sizes of the using DEMtransform data. Common resampling methods are Box Splines [26], interpolation [27] and the grid discrete wavelet [28]. 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 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 with multi-granularity beforehand and choose different grid sizes according to the requirements of different tasks. different tasks.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

10 of 30

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 m ants. The transfer probability pkij (t) of an ant k (0 ≤ k < m) moving from node i to an adjacent node j in the t round of iteration is expressed as: pkij (t)

=

  

β τα ij (t)ηij (t) β

∑s∈Jk (i) τα is (t)ηis (t)

, if jJk (i)

(8)

0, otherwise

where τij (t) is the residual pheromone of the edge < i, j >. ηij (t) is the heuristic function from node i to node j; when solving the TSP, ηij (t) = 1/cij . cij reflects the cost of the walk from node i to node j. α and β are the weights of the pheromone and the heuristic function, respectively. Adjacent nodes Jk (i) can be chosen by the ant k in node 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: m

τij (t + 1) = ρτij (t) +

∑ ∆τkij (t)

k =1

(9)

ISPRS Int. J. Geo-Inf. 2016, 5, 184

11 of 30

where ρ is the pheromone retention coefficient (0 ≤ ρ < 1). ∆τkij (t) represents the pheromones left at the edge of < i, j > by the ant k at the iteration of round t and is calculated as follows: ( ∆τkij

(t) =

1 , Lk

if < i, j >∈ Qk 0, otherwise

(10)

where Qk is the set of edges traversed by ant k and Lk 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:  n   W = ∑ w(i−1)i i=1 h − hmin l   wij = lijtraj + ε hijtraj max max − h

(11)

min

where wij is the cost of the edge < i, j >, lijtraj is the Euclidean length of the edge < i, j >, which approximately reflects the fuel consumption cost of the UAV, and hijtraj is the average terrain elevation
of node i and j, which is (alt i + altj /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. lmax is the largest Euclidean distance of the two adjacent nodes in the grid. Because planning space is irregular, we set lmax = grid width × 160 KM when the grid width is determined. hmax and hmin are the highest and lowest elevation of the planning space, respectively. ε 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 cij can only reflect the cost of node i to 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: ηij (t) = 2 +

Ti − Tj lmax

(12)

where Tj is the distance between node j and the target point and Ti is the distance between node i and the end point. If Tj > Ti , 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

ISPRS Int. J. Geo-Inf. 2016, 5, 184

12 of 30

proposed accordingly. The eight adjacent nodes of the current point i are A–H, and the implementation method is as follows: 1. 2. 3. 4.

Calculate the average terrain altitude of A–H points. Compare the average altitude with the altitude of point i; if greater than the threshold δ, the UAV at point i is located in the mountainous area. If the UAV is located in the mountainous region, the adjacent A–H points are sorted according to their altitudes, and the 2ε points with highest elevation are set to ηij (t) = 0. If the UAV is located in the non-mountainous region, all the neighbouring nodes are reachable.

After this treatment, η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: ∗ τij (t + 1) = ρτij (t) + ∆τbest (t) ij

(13)

where ∆τijbest∗ (t) is the pheromone increment of global optimal ant: ( ∗ ∆τbest (t) = ij

1 , Wbest∗

if < i, j >∈ Qbest 0, otherwise

(14)

where Qbest is the set of edges traversed by the global optimal ant. Wbest∗ 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 [τ min , τmax ]. ( 1 τmax = 1−1 ρ best W ∗ (15) τmin = τmax n where 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 i to the eight adjacent nodes. Additionally, a two-dimensional array is used to store the distance between 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 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 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.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

13 of 30

Given a UAV flight speed of v, the route point set that the ant walked is {r1 , r2 , · · · , rk } , k > 1; when k = 1, the UAV start time is t0 . If we know the time tk−1 of the ant at node k − 1, when the ant walks to route point k we look up the distance of rk−1 rk , named lrk−1 rk , and the current time tk = tk−1 + lrk−1 rk /v. The time of the ant at any node can thus be determined. The ants excluded the expanding nodes by the following steps: 1. 2. 3. 4.

When the ant k at node 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 Jk (i). By using the method mentioned above, we calculate the time the ant arrives at every expanding node j, named tj . 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 tj . 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, ISPRS Int. J. Geo-Inf. 5, 184 13 of 30 we compile two2016, 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 threat, we introduce theintroduce dynamicthe threat avoidance name the the dynamic threat, we dynamic threat algorithm avoidanceand algorithm and planning name thealgorithm planning MACOD. MACOD. algorithm 4.4. Route Optimization The route planned by the ACO algorithms above is composed of a series series of of broken broken line line segments. segments. To determine determine the the feasible feasibleroute, route,we weneed needto tooptimize optimizethe theroute routeby byaaseries seriesof ofoperations. operations. To 4.4.1. 4.4.1. Route Route Compression Compression To To avoid avoid frequent frequent turning turning of of the the UAV, UAV, itit isis necessary necessary to to compress compress the the waypoints, waypoints, as as shown shown in in Figure the broken Figure 10. 10. When When the broken line line is is compressed, compressed, it it will will be be shorter shorter than than the the original original path path planned planned by by the By compressing compressing the the route, route, errors errors are are introduced; however, the the grid grid search search algorithm. algorithm. By introduced; however, the route route is is optimized to a certain extent. optimized to a certain extent.

Figure 10. Length changed before and after compression. Figure 10. Length changed before and after compression.

When compressing the waypoints, the adjusted route should not cross the threatened area, and When waypoints, the significantly adjusted route should After not cross the threatened the average compressing altitude of thethe route has not been improved. deleting waypoint riarea, , the and the average altitude of the route has not been significantly improved. After deleting waypoint route should meet the criteria: (a) the adjusted route segments {ri−1 , ri+1 } should not be located in the rthreatened should thedifference criteria: (a) the adjusted route segments {ri−1{r, ri+,r1 } }should not be i , the routearea andmeet (b) the value of average altitude between i−1 i+1 and {ri−1 , ri , located in the threatened area and (b) the difference value of average altitude between { r 1 , ri+1 } ri+1 } should not exceed the specified threshold £. The average altitude is calculated asi−follows: 3 and { r , r , r } should not exceed the specified threshold £. The average altitude is calculated as i ithe +1 route with interval ∀ in Ʊ3 , then transform the interval points to Ψ and obtain i−1 for interpolate 3 3 follows: interpolate forthem. the route with interval ∀ in U , then transform the interval points to Ψ and the terrain altitude for obtain the terrain altitude This paper improves for thethem. Douglas-Peucker algorithm [36]; the improved algorithm steps are as This paper improves the Douglas-Peucker algorithm [36]; the improved algorithm steps are follows: as follows: 1. First, the algorithm compares the distance between each waypoint and the link line of two endpoints. If the distance is less than the tolerance ϵ, 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. 2. Select the point that has the farthest distance with the link line of two endpoints as the separation

ISPRS Int. J. Geo-Inf. 2016, 5, 184

1.

2.

14 of 30

First, the algorithm compares the distance between each waypoint and the link line of two endpoints. If the distance is less than the tolerance , 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 ISPRS Int. J. Geo-Inf. 2016, 5, 184 14 of 30 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, DTA is the length of the minimum straight route segment 2 (min[500,V VKTW ]) KTAS +speed, required by the Fly-by Fix; it is determined by the flight turning angle and other parameters R= tan(θ) × 68625.4 and is calculated as follows: { (16) ( αKTAS + VKTW ])2 (min[500,V R == R × tan DTA tan(θ (16) 2) × 68625.4 DTA = R × tan α2 where R is the turning radius, VKTAS is the true airspeed, VKTW is the tailwind, and θ is the slope where angle. R is the turning radius, VKTAS is the true airspeed, VKTW is the tailwind, and θ is the slope angle. UAV UAVflight flighttends tendsto toavoid avoidfrequent frequentrising risingor or descending. descending. If If aa UAV UAV needs needs to to fly fly from from the the waypoint waypoint A to waypoint B with different elevations, the UAV should rise/descend with maximum lifting lifting angle A to waypoint B with different elevations, the UAV should rise/descend with maximum to the horizontal line of line B and maintain level flight B. to B. angle to the horizontal of then B and then maintain leveltoflight

(a)

(b)

Figure 11. 11. (a) (a) Distance Distance of of turn turn anticipation anticipation model; model; (b) (b) Profile Profile view view for for RNP RNP segment segment width. width. Figure

4.4.3. Route Buffer Area 4.4.3. Route Buffer Area The route buffer area provides a vertical and horizontal buffer for the route for UAV safety. For The route buffer area provides a vertical and horizontal buffer for the route for UAV safety. For the the shadow regions shown in Figure 11b, the route buffer area regards the route as the centre line, shadow regions shown in Figure 11b, the route buffer area regards the route as the centre line, with a with a width of 4 × RNP, and the height is a rectangular buffer needed to override the obstacle (ROC). width of 4 × RNP, and the height is a rectangular buffer needed to override the obstacle (ROC). We do We do interpolation for the entire buffer. If the altitude difference between the interpolation point interpolation for the entire buffer. If the altitude difference between the interpolation point and the and the terrain is less than ROC, the height of the waypoint needs to be increased so that the route terrain is less than ROC, the height of the waypoint needs to be increased so that the route buffer area buffer area does not intersect the terrain and the threat zone. does not intersect the terrain and the threat zone. 5. Results 5. Resultsand andDiscussion Discussion In this this chapter, chapter, we we first first introduce introduce the the platform platform and and parameters parameters in in Section Section 5.1. 5.1. Pilot Pilot experiments experiments In with the route planner are reported in Section 5.2 to Section 5.4. In Section 5.5, we compare types with the route planner are reported in Sections 5.2–5.4. In Section 5.5, we compare six six types of of planners. To compare the different algorithms, we do not introduce dynamic threats and dynamic planners. To compare the different algorithms, we do not introduce dynamic threats and dynamic threat avoidance avoidance algorithms algorithms from from Section Section 5.1 5.1 to to Section Section 5.5. 5.5. The The planner planner in in this this paper paper is is currently currently threat named MACO-pl. In Section 5.6, we introduce the dynamic threat and evaluate the performance of named MACO-pl. In Section 5.6, we introduce the dynamic threat and evaluate the performance MACOD-pl. 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

ISPRS Int. J. Geo-Inf. 2016, 5, 184

15 of 30

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. Table 1. The UAV flight parameters. VKTAS (km/h)

VKTW (km/h)

Maximum Lifting Angle

h, 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 hmax and hmin 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: ρ = 0.7; weight coefficient: α = 1, β = 4; ant number: 30; iteration number: 200; node number n between the start point and end point; threshold δ = 10 m; local uplift threshold in route compression: £ = 10 m; and interpolation interval: ∀ = 100 m. This experiment used five coordinate groups, as shown in Table 2. Table 2. Experimental coordinate groups. Serial Number

1

2

3

4

5

Start Point

End Point

Lmin

HavT

(x, y, H)

(x, y, H)

(km)

(m)

99.96463◦

99.43347◦

26.88202◦

27.75616◦

110.613

2339.93

1819.58 m

2530.68 m

90.86051◦

94.46516◦

29.30996◦

29.42218◦

350.051

4492.95

342.512

4.49

1909.435

1186.37

2259.830

203.71

3573.99 m

2951.15 m

121.84255◦

119.05692◦

30.98066◦

32.95193◦

21.51 m

14.17 m

99.44990◦

118.45611◦

27.81014◦

25.09058◦

1723.22 m

397.35 m

109.78029◦

120.41494◦

18.85893◦

37.01404◦

445.07 m

132.31 m

ISPRS Int. J. Geo-Inf. 2016, 5, 184

16 of 30

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 A and BInt. were the start ISPRS J. Geo-Inf. 2016,and 5, 184end points, respectively. 16 of 30

Figure Figure 12. 12. Cartographic Cartographic Map-1. Map-1.

5.2. Get the Initial Route 5.2. Get the Initial Route The parameter ε balances the valley-following ability and the length of the route. Different The parameter ε 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 experimental results can be obtained by setting different coefficients. The initial altitude of the route is is the sum of the node elevation and the ROC value. The experimental results are shown in Figure 13 the sum of the node elevation and the ROC value. The experimental results are shown in Figure 13 and and Table 3. It can be concluded that when 𝜀 was relatively high, the valley-following performance Table 3. It can be concluded that when ε was relatively high, the valley-following performance of the of the algorithm was better. However, a higher 𝜀 will increase the difficulty for an ant to reach the algorithm was better. However, a higher ε will increase the difficulty for an ant to reach the end point, end point, which also increases the running time of the algorithm. When 𝜀 was 0, the algorithm which also increases the running time of the algorithm. When ε was 0, the algorithm calculated the calculated the shortest path; however, the path was longer than the shortest path in the continuous shortest path; however, the path was longer than the shortest path in the continuous world because world because of the impact of irregular topography and the grid space. When 𝜀 ≥ 3, a feasible of the impact of irregular topography and the grid space. When ε ≥ 3, a feasible solution may not solution may not be obtained. Suitable values of ε can be selected according to the application be obtained. Suitable values of ε can be selected according to the application requirements. For the requirements. For the following experiments we chose 𝜀 = 2. following experiments we chose ε = 2. Table 3. Comparison results for different coefficients. Table 3. Comparison results for different coefficients. Coefficients Planning Time (𝜀) (s) Coefficients Planning Time 0 (ε) 1.01(s) 0.5 0 1.09 1.01 1 0.5 1.20 1.09 1.5 1 1.29 1.20 1.29 2 1.5 1.37 1.37 2.5 2 1.55 1.55 3 2.5 N/A 3 N/A

Lroute (km) Lroute (km) 118.428

havT (m)havT (m) 2488.07

119.074 118.428 119.924 119.074 119.995 119.924 119.995 120.617 120.617 123.085 123.085

2351.98 2488.07 2152.15 2351.98 2138.84 2152.15 2138.84 2058.38 2058.38 2028.94 2028.94

Waypoint Number Waypoint Number 86 87 86 88 87 89 88 89 90 90 92 92

Compared with Compared with Lminwith HavT Compared Compared with L+7.07% H+6.33% min avT +7.65% +0.51% +7.07% +6.33% +8.42% -8.03% +7.65% +0.51% +8.48% -8.59% +8.42% −8.03% +8.48% −8.59% +9.04% -12.03% +9.04% −12.03% +11.28% -13.29% +11.28% −13.29%

Coefficients Planning Time (𝜀) (s) 0 1.01 0.5 1.09 1 1.20 1.5 1.29 2 1.37 ISPRS Int. J. Geo-Inf. 2016, 5, 184 2.5 1.55 3 N/A

ISPRS Int. J. Geo-Inf. 2016, 5, 184

Lroute (km)

havT (m)

118.428 119.074 119.924 119.995 120.617 123.085

2488.07 2351.98 2152.15 2138.84 2058.38 2028.94

(a)

Waypoint Number 86 87 88 89 90 92

Compared with Lmin +7.07% +7.65% +8.42% +8.48% +9.04% +11.28%

Compared with HavT +6.33% +0.51% -8.03% -8.59% -12.03% -13.29%

(b)

ISPRS Int. J. Geo-Inf. 2016, 5, 184

17 of 30

17 of 30 17 of 30

(c)

(d)

Figure 13. Experiment results with different height cost coefficients 𝜀. (a) 𝜀 = 0; (b) 𝜀 = 1; (c) 𝜀 = 2;

(c) different height cost coefficients (d)ε. (a) ε = 0; (b) ε = 1; (c) ε = 2; Figure 13. Experiment results with (d) 𝜀 = 2.5. (d) ε = Figure 2.5. 13. Experiment results with different height cost coefficients 𝜀. (a) 𝜀 = 0; (b) 𝜀 = 1; (c) 𝜀 = 2; (d) 𝜀 = 2.5. Figure 14 shows the convergence results of the MACO algorithm when 𝜀 = 2.

Figure Figure 14 shows the convergence results of the MACO algorithm when ε = 2. 14 shows the convergence results of the MACO algorithm when 𝜀 = 2.

(a)

(b)

(a)

(b)

(c)

(d)

Figure 14. Global optimal solutions for different iterations t. (a) t = 1; (c) (d)(b) t = 25; (c) t = 50; (d)

t = 100.

Figure 14. Global optimal solutions for different iterations t. (a) t = 1; (b) t = 25; (c) t = 50; (d)

Figure 14. Global optimal solutions for different iterations t. (a) t = 1; (b) t = 25; (c) t = 50; t = 100. 5.3. Route Optimization (d) t = 100. 5.3. Route Optimization 5.3.1. Route Compression results of Figure 13c were first optimized by compression. As shown in Figure 15 and Table 4, 5.3.1.The Route Compression different results were obtained by setting a different tolerance ϵ. The results of Figure 13clength were first optimized by compression. As shown in Figure 15 and Table 4, Table 4 shows that the of the compressed route may have been shorter than the shortest different results obtained by setting a different tolerance ϵ. path obtained bywere the ACO algorithm in Table 3. 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.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

18 of 30

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, different results were obtained by setting a different tolerance . 4 shows that the length of the compressed route may have been shorter than the shortest ISPRSTable Int. J. Geo-Inf. 2016, 5, 184 18 of 30 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.

(a)

(b)

(c)

(d)

Figure 15. Comparison results for different tolerances ϵ. (a) ϵ = 0.25 nmi; (b) ϵ = 0. 50 nmi; (c) ϵ = 0.75 Figure 15. Comparison results for different tolerances . (a)  = 0.25 nmi; (b)  = 0.50 nmi; nmi; (d) ϵ = 1.00 nmi. (c)  = 0.75 nmi; (d)  = 1.00 nmi. Table 4. Comparison results for different tolerances. Table 4. Comparison results for different tolerances. ϵ (nmi) Lroute (km) havT (m) Waypoint Number Compared with Lmin e (nmi) 120.129 Lroute (km) 2063.33 havT (m) Waypoint with Lmin 0.25 34 Number Compared +8.60% 0.500.25 117.172 16 34 +5.93% 120.129 2060.58 2063.33 +8.60% 117.172 2141.67 2060.58 +5.93% 0.750.50 116.003 9 16 +4.87% 116.003 2120.13 2141.67 +4.87% 1.000.75 115.061 6 9 +4.02% 1.00

115.061

2120.13

6

+4.02%

Compared with HavT Compared with HavT −11.82%

−−11.94% 11.82% −11.94% −8.47% −8.47% −9.39% −9.39%

5.3.2. Smooth Route Turning 5.3.2. Smooth Route Turning As Figure 16 shows, we used Equation (16) and Fly-by Fix to address all UAV turning corners. 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 With respect to height the results of Figure 16,turning the route length was 115.495 an average terrain of 2142.80 m. The curve is a after set ofadjustment route control points withkm an with interval average terrain height of 2142.80 m. The turning curve is a set of route control points with an interval of 50 m. of 50 m.

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

19 of 30 19 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184

19 of 30

(a)

(b)

Figure 16. Turning, Turning,falling fallingand andrising rising corner adjustment. Before turning corner adjustment; (b) Figure 16. corner adjustment. (a) (a) Before turning corner adjustment; (b) After After turning corner adjustment. turning corner adjustment.

5.3.3. Route Buffer Areas 5.3.3. Route Buffer Areas Figure 17a is the result of the introduction of route buffer zones. A route segment that does not Figure 17a is the result of the introduction of route buffer zones. A route segment that does not (a) is marked in red. Figure 17b is (b) meet the flight safety requirements the result after height adjustment. meet the flight safety requirements is marked in red. Figure 17b is the result after height adjustment. AfterFigure the introduction the route buffercorner zone, the route length wasturning adjusted to 115.609 km with 16. Turning,of falling and rising (a) Before corner adjustment; (b) an After the introduction of the route buffer zone,adjustment. 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 After turning corner adjustment. 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 figure of route and terrain. The Terrain Altitude-curve is the terrain altitude that the planning route passes through, andAreas the 2RNP Terrain Altitude-curve is the highest terrain altitude within the 2RNP 5.3.3. Buffer passesRoute 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 range. The initial was unsafe; we obtainedofa route safe route that couldAmeet safety requirements Figure 17a is route the result of the introduction buffer zones. routethe segment that does not after adjusting the elevation. after adjusting the elevation. meet the flight safety requirements is marked in red. Figure 17b is the result after height adjustment. route optimization, LUSbuffer = 0; however, average of the route changed AfterAfter the introduction of the route zone, thethe route lengthheight was adjusted to 115.609 km and withthe an average terrain passed the UAV also height changed. Errors introduced by 18 theisoptimization average height height of 3228.90 m, through and the by terrain average was 2142.50 m. Figure the vertical processofmay cause route to pass through a threat area, mentioned Section 3.4.1, and the figure route and aterrain. The Terrain Altitude-curve is as thewas terrain altitudein that the planning route threat area should thus expanded when it is mappedisinto planning space. The within expansion width passes through, and thebe 2RNP Terrain Altitude-curve the the highest terrain altitude the 2RNP of the threat zone should be greater than the navigation performance error (0.5 nmi); the width of the range. The initial route was unsafe; we obtained a safe route that could meet the safety requirements ◦ grid width planning space was 1 node. corresponding expansion in the 0.01 after adjusting the elevation.

(a)

(b)

Figure 17. Route optimization with route buffer areas. (a) Before height adjustment; (b) After

height adjustment.

(a)

(b)

Figure optimization with route buffer areas. (a) Before height adjustment; (b) After Figure 17. 17. Route Route optimization with route buffer areas. (a) Before height adjustment; (b) After height adjustment. height adjustment.

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

20 of 30 20 of 30

Figure Figure 18. 18. Vertical Vertical flight flight path. path.

After route optimization, 5.4. Multi-Granularity Results LUS = 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 In Section 3.4.2,may we cause introduced planning space forindifferent optimization process a routeatomulti-granularity pass through a threat area, as was suitable mentioned Section applications. Based on different grids, the experimental results are shown in Table 5. 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); Table 5. Comparison of results with different granularity. the width of the corresponding expansion in the 0.01° grid width planning space was 1 node. Grid Width

e (nmi)

Processing Time (s)

5.4. Multi-Granularity Results ◦

Lroute (km)

havT (m)

Compared with Lmin

Compared with HavT

0.007 0.75 2.15 115.117 1975.70 +4.07% −15.57% ◦ 0.007 0.50 2.15 116.899 1947.84 planning +5.68% −for 16.76% In Section 3.4.2, we introduced a multi-granularity space suitable different 0.006◦ 0.75 2.73 114.404 1979.44 +3.43% −15.41% ◦ applications. Based on different grids, the experimental results are shown in Table 5. 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% Table 5. Comparison of results with different granularity. 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% Grid Width ϵ (nmi) Processing Time (s) L (km) h (m) Compared with L Compared with HavT route avT min 0.003◦ 0.75 6.18 114.498 1961.64 +3.51% −16.17% 0.007° 0.75 2.15 115.117 1975.70 +4.07% −15.57% ◦ 0.003 0.50 6.18 115.300 1925.10 +4.24% −17.73% 0.007° 0.50 2.15 116.899 1947.84 +5.68% −16.76% 0.002◦ 0.75 7.37 115.033 1973.27 +4.00% − 15.67% 0.006° 0.75 2.73 114.404 1979.44 +3.43% −15.41% 0.002◦ 0.50 7.37 115.711 1961.85 +4.61% − 16.16% 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% As the grid reduced and 115.603 the planning space became closer to the real world, the 0.005° 0.50width was 3.21 1905.73 +4.51% −18.56% 0.004° 0.75 was more4.24 114.403 −16.41% route we obtained consistent with reality.1955.89 The algorithm+3.43% ran slower; however, when the 0.004° 0.50 4.24 115.523 1923.16 +4.44% −17.81% grid width was reduced to a certain value limited by factors such as the input parameters and the 0.003° 0.75 6.18 114.498 1961.64 +3.51% −16.17% experimental terrain data source, the results no longer improved. In+4.24% the experiment scene shown in 0.003° 0.50 6.18 115.300 1925.10 −17.73% ◦ and  was 0.50 nmi. The following Figure 19, a better effect was achieved when grid width was 0.005 0.002° 0.75 7.37 115.033 1973.27 +4.00% −15.67% 0.002° 0.50 7.37 parameters. 115.711 1961.85 +4.61% −16.16% experiments are based on these

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 ϵ was 0.50 nmi. The following experiments are based on these parameters.

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

21 of 30 21 of 30

(b)

(a)

(c)

Figure 19. 19. (a) (a) The Theoptimal optimalresult resultfor forgrid gridwidth width = 0.005°, ◦ ,  ϵ== 0.50 Figure = 0.005 0.50 nmi, nmi, εε ==2, 2,terrain terrain average average altitude = 1,905.73 m, and average altitude of the route after introducing the route buffer = 2698.37 altitude = 1905.73 m, and average altitude of the route after introducing the route buffer = 2698.37 m; m; (b) are details. local route details. (b,c) areand local(c) route

5.5. Comparison 5.5. Comparison Experiments Experiments On our ourvirtual virtual globe platform we implemented also implemented theplanner route of planner ofACO the pheromone basic ACO On globe platform we also the route the basic pheromone model named ACO-pl and the Rank-based pheromone model named RAS-pl. The model named ACO-pl and the Rank-based pheromone model named RAS-pl. The number of retaining number of retaining RAS-pl was valley-following 5. We used the local valley-following described in ants of RAS-pl was 5.ants We of used the local method described in method Section 4.2.2 and the Section 4.2.2 and the shortest distance from the target as the heuristic function to realize an A* based shortest distance from the target as the heuristic function to realize an A* based planner, named A*-pl. planner, named A*-pl. Additionally, we implemented the planners based on PSO and EA, named PSO-pl and EA-pl, Additionally, we these implemented the planners on PSO and EA, encode named method PSO-pl and EA-pl, respectively. Because two methods use thebased dimension-reduction to achieve 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 evolutionary computation, the local valley-following method in Section 4.2.2 wasaltitude not applicable. This study used the minimal flight altitude method in [12] to update the average of the route. study used theinterval minimaldifference flight altitude method in [12] to updateinterval the average altitude for of the Because of the between adjacent waypoints, interpolation the route. route Because of the interval difference between adjacent waypoints, interval interpolation for theweight route was needed when calculating the average altitude. In the following experiments, the inertial was needed calculating average the following the inertial weight coefficient ofwhen PSO-pl was linearthe from 0.9 toaltitude. 0.4. The In iteration count ofexperiments, PSO-pl and EA-pl was 200 and 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. the population size was 30. When implementing the six types of planners, this study used data structures to store When implementing the sixthe types of planners, thispre-calculated study used data structures to each storenode prepre-calculated distances between adjacent nodes and distances between calculated distances between the adjacent nodes and pre-calculated distances between each node and and the end node. The adjacent waypoints in MACO-pl, ACO-pl, RAS-pl and A*-pl are adjacent the end node. The adjacent waypoints in MACO-pl, ACO-pl, RAS-pl and A*-pl are adjacent nodes in nodes in the grid planning space, as shown in Figure 20a. The algorithm proposed and used in this the grid planning space, as shown in FigureWhile 20a. The algorithm proposedinand usedand in this study study achieved computational acceleration. the adjacent waypoints PSO-pl EA-pl are 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 movenot in necessarily adjacent nodes, as shown in Figure 20b, the waypoints P and Q can only move in two two directions. At this point, P and Q are not adjacent nodes and thus distances must be calculated. directions.the Atwaypoints this point,QPand andHQ not adjacent nodes thus must be Qcalculated. Although areare adjacent nodes, we can and obtain thedistances distance between and H by Although the waypoints Q and H are adjacent nodes, we can obtain the distance between Q and H a table lookup. by a table lookup.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

22 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184

22 of 30

(a)

(b)

Figure MACO, ACO, Figure 20. 20. The The relationship relationshipbetween betweenwaypoint waypointand andnode nodefor fordifferent differentalgorithms. algorithms.(a)(a) MACO, ACO, RAS, RAS, A*; A*; (b) (b) PSO, PSO, EA. EA.

The comparison experiments used the same grid size and UAV parameters and were conducted 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 and Figure 21. for different coordinate groups. The results are shown in Table 6 and Figure 21. Table 6. Planning results of 1–5 coordinate groups. Table 6. Planning results of 1–5 coordinate groups. Serial Method Processing Time (s) Number Serial Processing Method 1NumberMACO-pl 3.21 Time (s) 1 ACO-pl 10.07 1 MACO-pl 1 RAS-pl 5.66 3.21 1 ACO-pl 1 A*-pl 2.16 10.07 1 1 PSO-pl RAS-pl 9.78 5.66 1 1 EA-pl A*-pl 10.84 2.16 1 MACO-plPSO-pl 2 12.39 9.78 1 2 ACO-pl EA-pl 38.4810.84 2 MACO-pl 2 RAS-pl 20.4012.39 2 2 A*-pl ACO-pl 148.1438.48 2 2 PSO-pl RAS-pl 40.0220.40 2 2 EA-pl A*-pl 41.28148.14 2 MACO-plPSO-pl 3 10.9640.02 2 3 ACO-pl EA-pl 36.4241.28 3 RAS-pl 18.6710.96 3 MACO-pl 3 A*-pl ACO-pl 116.2336.42 3 3 PSO-pl RAS-pl 39.5718.67 3 3 EA-pl A*-pl 38.95116.23 3 4 80.1239.57 3 MACO-plPSO-pl 4 ACO-pl EA-pl 201.1538.95 3 4 RAS-pl 4 MACO-pl 149.9880.12 4 A*-pl ACO-pl N/A201.15 4 4 PSO-pl RAS-pl 278.26 4 149.98 4 EA-pl A*-pl 299.38N/A 4 5 78.84278.26 4 MACO-plPSO-pl 5 ACO-pl 197.43 4 EA-pl 299.38 5 RAS-pl 5 MACO-pl 136.0178.84 5 A*-pl ACO-pl N/A197.43 5 5 PSO-pl 301.44 5 RAS-pl 136.01 5 EA-pl A*-pl 279.64N/A 5

5 5

PSO-pl EA-pl

301.44 279.64

𝐋𝐫𝐨𝐮𝐭𝐞 (km)

L route (km) 115.603 118.831 115.603 117.267 118.831 116.875 117.267 118.679 116.875 119.441 118.679 398.782 119.441 389.907 398.782 399.823 389.907 396.934 399.823 386.679 396.934 385.041 386.679 362.433 385.041 368.345 365.184 362.433 363.245 368.345 366.176 365.184 364.504 363.245 2115.333 366.176 2135.830 364.504 2136.106 2115.333

Compared with Lmin Compared havT (m) 1905.73 +4.51% with Lmin 2262.32 +7.43% 1905.73 +4.51% 1914.28 +6.02% 2262.32 +7.43% 1916.29 +5.66% 1914.28 +6.02% 2177.42 +7.29% 1916.29 +5.66% 2190.32 +7.98% 2177.42 +7.29% 3471.04 +13.92% 2190.32 +7.98% 4074.93 +11.39% 3471.04 +13.92% 3623.61 +14.22% 4074.93 +11.39% 3681.73 +13.39% 3623.61 +14.22% 4081.56 +10.46% 3681.73 +13.39% 4094.39 +10.00% 4081.56 +10.46% 2.85 +5.82% 4094.39 +10.00% 2.98 +7.54% 2.932.85 +6.62% +5.82% 2.942.98 +6.05% +7.54% 4.012.93 +6.91% +6.62% 3.952.94 +6.42% +6.05% 1043.90 +10.78% 4.01 +6.91% 1097.04 +11.85% 3.95 +6.42% 1063.77 +11.87% 1043.90 +10.78%

havT (m)

2135.830

1097.04

+11.85%

2104.857 2136.106 2109.891 2371.012 2104.857 2401.023 2109.891 2390.174 2371.012

1137.25 1063.77 1131.32 84.59 1137.25 116.56 1131.32 85.99 84.59

+10.23% +11.87% +10.50% +4.92% +10.23% +6.25% +10.50% +5.76% +4.92%

2401.023

116.56

+6.25%

2398.322 2390.174 2388.459

127.76 85.99 134.31

+6.13% +5.76% +5.69%

2398.322 2388.459

127.76 134.31

+6.13% +5.69%

Compared with HavT Compared −18.56% with HavT −3.32% −−18.19% 18.56% − 3.32% −18.10% −−6.94% 18.19% −−6.39% 18.10% − 6.94% −22.74% − 6.39% −9.30% −−19.35% 22.74% − 9.30% −18.05% −−9.15% 19.35% −−8.87% 18.05% − 9.15% −36.53% − 8.87% −33.63% −−34.74% 36.53% −−34.52% 33.63% −−10.69% 34.74% −−12.03% 34.52% −−12.01% 10.69% −−7.53% 12.03% −−10.33% 12.01%

−7.53% −−4.14% 10.33% −4.64% −58.48% − 4.14% −42.78% − 4.64% −−57.79% 58.48% −42.78% −−37.28% 57.79% −34.07%

−37.28% −34.07%

ISPRS Int. J.J. Geo-Inf. Geo-Inf. 2016, 2016, 5, 5, 184 184 ISPRS

23 of of 30 30 23

(a)

(b)

(c) Figure Comparisonwith Figure 21. 21. Planning Planningresults resultsofof1–5 1–5coordinate coordinategroups. groups.(a)(a)Processing Processingtime time (s); (s); (b) (b) Comparison with L ; (c) Comparison with H . avT Lmin ; (c)min Comparison with HavT .

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

ISPRS Int. J. Geo-Inf. 2016, 5, 184

24 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

24 of 30 24 of 30

(a) (a)

(b) (b)

(c) (d) (c) (d) Figure 22. Planning results of 2–5 coordinate groups using MACO-pl. (a) Coordinate group 2; (b) Figure 22. Planning results of 2–5 coordinate groups using MACO-pl. (a) Coordinate group 2; Figure 22. Planning of 2–5 coordinate using MACO-pl. Coordinate group 3;results (c) Coordinate group 4; groups (d) Coordinate group 5. (a) Coordinate group 2; (b) (b) Coordinate group 3; (c) Coordinate group 4; (d) Coordinate group 5. Coordinate group 3; (c) Coordinate group 4; (d) Coordinate group 5.

(a) (a)

(b) (b)

(c) (c)

(d) (d)

Figure 23. Local details of the routes of 2–5 coordinate groups. (a) Coordinate group 2; (b) Figure 23. Local details of the routesgroup of 2–54coordinate groups. (a) Coordinate group 2; (b) Coordinate group 3; (c) Coordinate ; (d) Coordinate group 5. Figure 23. Local details of the routes of 2–5 coordinate groups. (a) Coordinate group 2; (b) Coordinate Coordinate group 3; (c) Coordinate group 4; (d) Coordinate group 5. group 3; (c) Coordinate group 4; (d) Coordinate group 5.

ISPRS Int. J. Geo-Inf. 2016, 5, 184

25 of 30

After marking the radar zone, no-fly zone, and missile threat zone, we did experiments on coordinate 1, 5, 2,184 4, and 5. Table 7 and Figure 24 show the results of different planners. ISPRS Int. J. groups Geo-Inf. 2016, 25 of 30 We demonstrated that the MACO-pl had a better solution. Figure 25 shows local results by MACO-pl marking the radar zone, no-fly zone, and missile threat zone, we could did experiments onthe of four After coordinate groups after plotting threatened areas. The planned routes make use of coordinate groups 1, 2, 4, and 5. Table 7area. and Figure 24 show the results of different planners. We terrain condition to avoid the threatened demonstrated that the MACO-pl had a better solution. Figure 25 shows local results by MACO-pl of Table 7. Planning resultsareas. in threat four coordinate groups after plotting threatened Theenvironments. planned routes could make use of the terrain condition to avoid the threatened area. Serial Number 1 Serial 1 Number 11 11 11 11 12 12 22 2 2 2 2 2 2 2 4 2 4 4 44 44 44 44 45 55 55 55 55 5 5 5

Method MACO-pl

Processing Lroute havT (m) LFA (km) Time (s)7. Planning (km) Table results in threat environments. 5.63

Method Time (s) ACO-plProcessing 17.02

RAS-pl MACO-pl A*-pl ACO-pl RAS-pl PSO-pl A*-pl EA-pl PSO-pl MACO-pl EA-pl ACO-pl MACO-pl RAS-pl ACO-pl A*-pl RAS-pl PSO-pl A*-pl EA-pl PSO-pl MACO-pl EA-pl ACO-pl MACO-pl RAS-pl ACO-pl A*-pl RAS-pl PSO-pl A*-pl EA-pl PSO-pl EA-pl MACO-pl MACO-pl ACO-pl ACO-pl RAS-pl RAS-pl A*-pl A*-pl PSO-pl PSO-pl EA-pl EA-pl

9.39 5.63 5.27 17.02 9.39 29.84 5.27 31.15 29.84 20.85 31.15 75.35 20.85 37.28 75.35 326.92 37.28 91.49 326.92 99.02 91.49 127.39 99.02 391.50 127.39 271.98 391.50 N/A 271.98 401.32 N/A 384.56 401.32 384.56 84.26 84.26 209.92 209.92 140.15 140.15 N/A N/A 329.48 329.48 286.77 286.77

136.480

2725.64

0

L138.737 route (km)

h2870.48 avT (m)

LFA (km) 0

137.158 136.480 136.497 138.737 137.158 138.528 136.497 137.917 138.528 421.922 137.917 431.943 421.922 423.871 431.943 422.775 423.871 430.563 422.775 428.719 430.563 2250.836 428.719 2301.760 2250.836 2284.541 2301.760

2757.80 2725.64 2789.72 2870.48 2757.80 2955.42 2789.72 2918.81 2955.42 4379.01 2918.81 4687.66 4379.01 4470.49 4687.66 4450.11 4470.49 4716.47 4450.11 4671.99 4716.47 1086.37 4671.99 1143.74 1086.37 1109.44 1143.74 1109.44

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

2284.541

2304.283 2300.575 2304.283 2300.575 2438.912 2438.912 2494.026 2494.026 2477.678

1217.07 1203.64 1217.07 1203.64 86.22 86.22 119.14 119.14 86.85 86.85

2477.678

2489.205 2489.205 2503.418

129.49

129.49 136.61 136.61

2503.418

(a)

(b) Figure 24. Cont.

0 0 0 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

0 0 0

Compared with Lmin +23.39% Compared +25.43% with Lmin +24.00% +23.39% +23.40% +25.43% +24.00% +25.24% +23.40% +24.68% +25.24% +20.53% +24.68% +23.39% +20.53% +21.09% +23.39% +20.78% +21.09% +23.00% +20.78% +22.47% +23.00% +17.88% +22.47% +20.55% +17.88% +19.64% +20.55% +19.64% +20.68%

Compared with HavT +16.48% Compared +22.67% with HavT +17.86% +16.48% +19.22% +22.67% +17.86% +26.30% +19.22% +24.74% +26.30% −2.54% +24.74% +4.33% −2.54% − 0.50% +4.33% − 0.95% −0.50% +4.97% −0.95% +3.98% +4.97% − 8.43% +3.98% − 3.59% −8.43% − 6.49% −3.59% −6.49% +2.59%

+20.48% +20.68% +20.48% +7.92% +7.92% +10.36% +10.36% +9.64% +9.64%

+1.46% +2.59% −+1.46% 57.68% −57.68% − 41.51% −41.51% − 57.37% −57.37%

+10.15%

−36.43% −36.43% − 32.94%

+10.15% +10.78% +10.78%

−32.94%

ISPRS Int. J. Geo-Inf. 2016, 5, 184

26 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184

26 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184

26 of 30

(c) (c) Figure 24. Planning results under threat environments. (a) Processing time (s); (b) Comparison Figure 24. Planning results under threat environments. (a) Processing time (s); (b) Comparison with Figure Planning results under (a) Processing time (s); (b) Comparison with Lmin24. ; (c) Comparison with threat HavTenvironments. . Lmin ; (c) Comparison with HavT . with Lmin ; (c) Comparison with HavT .

(a) (a)

(c)

(b) (b)

(d)

(c) Figure 25. Low altitude penetration routes using MACO-pl of different(d) cases. (a) Case 1; (b) Case 2; Figure 25. Low altitude penetration routes using MACO-pl of different cases. (a) Case 1; (b) Case 2; (c) Case 4; (d) Case 5.penetration routes using MACO-pl of different cases. (a) Case 1; (b) Case 2; Figure 25. Low altitude (c) Case 4; (d) Case 5. (c) Case 4; (d) Case 5. 5.6. Experiments on Dynamic MACO

5.6. Experiments on Dynamic MACO This section introduces the dynamic threat area discussed in Section 3.3. Based on MACOD-pl, 5.6. Experiments on Dynamic MACO we used coordinate groupthe 1 todynamic test the algorithm. Considering thatin the route optimization This section introduces threat area discussed Section 3.3. Based algorithms on MACOD-pl, This section introduces the dynamic threat area discussed in Section 3.3. Based on MACOD-pl, would bring approximately 1%–3% error to the initial route length, which is approximately km we used coordinate group 1 to test the algorithm. Considering that the route optimization1–4 algorithms we used coordinate group 1 tothe test the algorithm. Considering that the route optimization algorithms under the current scenario, maximum time error was approximately 40 s, and the moving speed would bring approximately 1%–3% error to the initial route length, which is approximately 1–4 km would bring approximately errorwas to the initial route length, which is approximately 1–4 km of the dynamic threat area 1%–3% in the figure approximately 12.6 m/s. The maximum error range was under thethe current scenario, the maximum time 40 s,s,and speed under current scenario, the maximum timeerror errorwas wasapproximately approximately andthe the moving moving speed of thus approximately 510 m. To ensure route safety, the boundary of the 40 dynamic threat area was theofdynamic threat area in the figure was approximately 12.6 m/s. The maximum error range was thus the dynamic area in theinfigure approximately 12.6 expanded bythreat 0.01°. As shown Figurewas 26 and Table 8, when them/s. startThe timemaximum of the UAVerror was range not thewas approximately 510 m. To meet ensure the boundary of theand dynamic threatroutes area wasarea expanded the UAV could theroute threatsafety, atroute different timesthe and places, were thuswas thussame, approximately 510 m. To ensure safety, boundary ofthe theplanned dynamic threat ◦ . As shown in Figure 26 and Table 8, when the start time of the UAV was not the same, the UAV by expanded 0.01different. When the start time of the UAV was between 0 and 2000 s, the dynamic threat area would by 0.01°. As shown in Figure 26 and Table 8, when the start time of the UAV was not the block UAV from flyingthe in the valley area. When the start of the UAV 3000 s, the dynamic could meet the threat atmeet different times places, and the time planned routes were thus different. When same, thethe UAV could threat atand different times and places, and the was planned routes were thus threat area would pass through the valley area and the result of route planning was similar to the thedifferent. start timeWhen of thethe UAV was between 0 andwas 2000between s, the dynamic threat area would threat block the from start time of the UAV 0 and 2000 s, the dynamic areaUAV would flying inthe theUAV valley area. When thevalley start area. time When of thethe UAV was 3000 s, the dynamic threat would block from flying in the start time of the UAV was 3000 s, thearea dynamic threat area would passarea through and the result of route planning similar to thethe pass through the valley and the the valley result area of route planning was similar to thewas result without threat. Because of the introduction of dynamic threats, the dynamic threat avoidance algorithm needed

ISPRS Int. J. Geo-Inf. 2016, 5, 184

27 of 30

ISPRS Int. J. Geo-Inf. 2016, 5, 184

27 of 30

to consume many computing operating of efficiency the MACOD-pl algorithm result without the threat. resources, Because ofand thethe introduction dynamicofthreats, the dynamic threat was muchavoidance lower than that of the MACO-pl algorithm. 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. Table 8. The results of different start times under dynamic threats. Start Time (s) Start time 0 (s) 1000 0 2000 1000 3000 2000 3000

Table 8. The results of different start times under dynamic threats. Planning Time (s) Lroute (km) havT (m) Compared with Lmin Compared with HavT Planning time Lroute (km) havT (m) Compared with Lmin Compared − with HavT 10.01 133.850 2112.65 +21.01% 9.71% (s) 12.14 135.633 2724.18 +22.62% +16.42% 10.01 133.850 2112.65 +21.01% −9.71% 12.89 136.204 2725.49 +23.14% +16.48% 12.14 135.633 2724.18 +22.62% +16.42% 9.49 115.377 1906.71 +4.31% −18.51% 12.89 136.204 2725.49 +23.14% +16.48% 9.49 115.377 1906.71 +4.31% −18.51%

(a)

(c)

(b)

Figure 26. The results of different start times under dynamic threats. (a) Start time = 0 s; (b) Start Figure 26. The results of different start times under dynamic threats. (a) Start time = 0 s; (b) Start time = 1000 s; (c) Start time = 3000 s. time = 1000 s; (c) Start time = 3000 s.

6. Conclusions and Future Work

6. Conclusions and Future Work

This paper presented an improved ACO-based planner based on the virtual globe platform to

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

ISPRS Int. J. Geo-Inf. 2016, 5, 184 ISPRS Int. J. Geo-Inf. 2016, 5, 184

28 of 30

28 of 30

Many other types of planning algorithms exist; for example, some algorithms use a rotated

other algorithms to thetovirtual platformofand obtainFurther betterresearch results would coordinate frame reduceglobe the dimension the efficiently problem [12,18]. on howbe to worthwhile, apply as wellother as further research on the model of different dynamic threats and the efficient use of dynamic algorithms to the virtual globe platform and efficiently obtain better results would be as well as further research on the model of different dynamic threats and the efficient threatworthwhile, avoidance algorithms. use of dynamic threat avoidance algorithms. Acknowledgments: The authors would like to thank the support of Shenzhen Government, Project No. GRCK2015092515503432. We would likelike to thank NASA for providing data. Government, Project No. Acknowledgments: The authors also would to thank the support of Shenzhen GRCK2015092515503432. We would also like to thank NASA for providing data.

Author Contributions: Ming Zhang and Yuesheng Zhu designed the project; Ming Zhang, Chen Su, Yuan Liu and Mingyuan Hu preparedMing the manuscript. All the Zhu authors havethe participated in Zhang, the project. Author Contributions: Zhang and Yuesheng designed project; Ming Chen Su, Yuan Liu and Mingyuan Hu prepared the manuscript. All the authors have participated in the project.

Conflicts of Interest: The authors declare no conflicts of interest. Conflicts of Interest: The authors declare no conflicts of interest.

Appendix A Appendix A

Figure A1. Threat modelling in the route planning system in which a radar threat zone is marked

Figure A1. Threat modelling in the route planning system in which a radar threat zone is interactively. marked interactively. ISPRS Int. J. Geo-Inf. 2016, 5, 184

29 of 30

Figure A2. Screenshot of a flight simulation in the route planning system in which a simulated flight

Figure A2. Screenshot of a flight simulation in the route planning system in which a simulated flight experiment is performed for the saved route. experiment is performed for the saved route. References 1. 2.

3.

Kimon, P.V.; George, J.V. Handbook of Unmanned Aerial Vehicles; Springer: Heidelberg, Germany, 2015. Neto, A.A.; Macharet, D.G.; Campos, M.F.M. Feasible RRT-based path planning using seventh order Bézier curves. In Proceedings of the 2010 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Taipei, Taiwan, 18–22 October 2010. Goerzen, C.; Kong, Z.; Mettler, B. A survey of motion planning algorithms from the perspective of

ISPRS Int. J. Geo-Inf. 2016, 5, 184

29 of 30

References 1. 2.

3. 4. 5. 6.

7. 8. 9. 10. 11.

12. 13.

14.

15. 16. 17. 18. 19. 20. 21.

22. 23.

Kimon, P.V.; George, J.V. Handbook of Unmanned Aerial Vehicles; Springer: Heidelberg, Germany, 2015. Neto, A.A.; Macharet, D.G.; Campos, M.F.M. Feasible RRT-based path planning using seventh order Bézier curves. In Proceedings of the 2010 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Taipei, Taiwan, 18–22 October 2010. Goerzen, C.; Kong, Z.; Mettler, B. A survey of motion planning algorithms from the perspective of autonomous UAV guidance. J. Intell. Robot. Syst. 2010, 57, 65–100. [CrossRef] Nelson, D.R.; Barber, D.B.; McLain, T.W.; Beard, R.W. Vector field path following for miniature air vehicles. IEEE Trans. Robot. 2007, 23, 519–529. [CrossRef] Chen, H.; Chang, K.; Agate, C.S. UAV path planning with Tangent-plus-Lyapunov vector field guidance and obstacle avoidance. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 840–856. [CrossRef] Lee, J.; Pippin, C.; Balch, T. Cost based planning with RRT in outdoor environments. In Proceedings of the IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS 2008), Nice, France, 22–26 September 2008. Wen, N.; Zhao, L.; Su, X.; Ma, P. UAV online path planning algorithm in a low altitude dangerous environment. IEEE/CAA J. Autom. Sin. 2015, 2, 173–185. Wen, N.; Zhao, L.; Su, X.; Ma, P. Online UAV path planning in uncertain and hostile environments. Int. J. Mach. Learn. Cyber. 2015, 1–19. [CrossRef] Ergezer, H.; Leblebicioglu, K. Path planning for UAVs for maximum information collection. IEEE Trans. Aerosp. Electron. Syst. 2013, 49, 502–520. [CrossRef] Zheng, C.; Li, L.; Xu, F.; Sun, F.; Ding, M. Evolutionary route planner for unmanned air vehicles. IEEE Trans. Robot. 2005, 21, 609–620. [CrossRef] Mittal, S.; Deb, K. Three-dimensional offline path planning for UAVs using multiobjective evolutionary algorithms. In Proceedings of the IEEE Congress on Evolutionary Computation (CEC 2007), Singapore, 25–28 September 2007. Yang, P.; Tang, K.; Lozano, J.A.; Cao, X. Path planning for single unmanned aerial vehicle by separately evolving waypoints. IEEE Trans. Robot. 2015, 31, 1130–1146. [CrossRef] Fu, Y.; Ding, M.; Zhou, C.; Hu, H. Route planning for Unmanned Aerial Vehicle (UAV) on the sea using hybrid differential evolution and quantum-behaved particle swarm optimization. IEEE Trans. Syst. Man Cybern. Syst. 2013, 43, 1451–1465. [CrossRef] Ling, X.; Hao, Y. Effective 3-D path planning for UAV in presence of threat netting. In Proceedings of the 2015 Fifth International Conference on Communication Systems and Network Technologies (CSNT), Gwalior, India, 4–6 April 2015. Duan, H.B.; Yu, Y.X.; Zhang, X.Y.; Shao, S. Three-dimension path planning for UCAV using hybrid meta-heuristic ACO-DE algorithm. Simul. Model. Pract. Theory 2010, 18, 1104–1115. [CrossRef] Garcia, M.; Viguria, A.; Ollero, A. Dynamic graph-search algorithm for global path planning in presence of hazardous weather. J. Intell. Robot. Syst. 2013, 69, 285–295. [CrossRef] Zhou, Y.; Bao, Z.; Wang, R.; Qiao, S.; Zhou, Y. Quantum wind driven optimization for unmanned combat air vehicle path planning. Appl. Sci. 2015, 5, 1457–1483. [CrossRef] Yao, P.; Wang, H. Dynamic Adaptive Ant Lion Optimizer applied to route planning for unmanned aerial vehicle. Soft Comput. 2016. [CrossRef] Le, Y.; Peng, G. Google Earth as a virtual globe tool for Earth science applications at the global scale: progress and perspectives. Int. J. Remote Sens. 2012, 33, 3966–3986. Zheng, G.; Zheng, Y. Radar netting technology & its development. In Proceedings of the 2011 IEEE CIE International Conference on Radar, Chengdu, China, 24–27 October 2011. Bell, D.G.; Kuehnel, F.; Maxwell, C.; Kim, R.; Kasraie, K.; Gaskins, T.; Hogan, P.; Coughlan, J. NASA world wind: Opensource GIS for mission operations. In Proceedings of the 2007 IEEE Aerospace Conference, Big Sky, MT, USA, 3–10 March 2007. United States Standard for Performance Based Navigation (PBN) Instrument Procedure Design Document Information. Available online: http://www.faa.gov/regulations_policies/ (accessed on 10 October 2015). Mingan, Z. Killing zone boundary and defense efficiency of surface-air missile. Tactical Missile Technol. 1992, 1, 1–8. (In Chinese)

ISPRS Int. J. Geo-Inf. 2016, 5, 184

24. 25. 26. 27. 28. 29. 30. 31.

32. 33. 34. 35. 36. 37.

30 of 30

Haifeng, L. A Study on Spatial Modeling Method of Threats in Mission Rehearsal of Tactical Aviation. Master’s Thesis, National University of Defense Technology, Changsha, China, 2006. (In Chinese) Skolnik, M.I. Introduction to Radar Systems, 3rd ed.; McGraw-Hill Book Company: New York, NY, USA, 2002. Hebert, D.J. A box spline subdivision pyramid algorithm. Appl. Math. Lett. 1999, 12, 57–62. [CrossRef] Donoho, D.L.; Yu, P.Y. Nonlinear pyramid transforms based on median-interpolation. SIAM J. Math. Anal. 2000, 31, 1030–1061. [CrossRef] Vishwanath, M. The recursive pyramid algorithm for the discrete wavelet transform. IEEE Trans. Signal Process. 1994, 42, 673–676. [CrossRef] Dorigo, M.; Maniezzo, V.; Colorni, A. Ant system: Optimization by a colony of cooperating agents. IEEE Trans. Syst. Man Cybern. Cybern. 1996, 26, 29–41. [CrossRef] [PubMed] Dorigo, M.; Gambardella, L.M. Ant colony system: A cooperative learning approach to the traveling salesman problem. IEEE Trans. Evol. Comput. 1997, 1, 53–66. [CrossRef] Adubi, S.A.; Misra, S. A comparative study on the ant colony optimization algorithms. In Proceedings of the 11th International Conference on Electronics, Computer and Computation (ICECCO), Abuja, Nigeria, 29 September–1 October 2014. Albinati, J.; Oliveira, S.E.; Otero, F.E.; Pappa, G.L. An ant colony-based semi-supervised approach for learning classification rules. Swarm Intell. 2015, 9, 315–341. [CrossRef] Farahnakian, F.; Ashraf, A.; Pahikkala, T.; Liljeberg, P.; Plosila, J.; Porres, I.; Tenhunen, H. Using ant colony system to consolidate VMs for green cloud computing. IEEE Trans. Serv. Comput. 2015, 8, 187–198. Wang, Z.; Xing, H.; Li, T.; Yang, Y.; Qu, R.; Pan, Y. A modified ant colony optimization algorithm for network coding resource minimization. IEEE Trans. Evol. Comput. 2015, 1, 1–18. [CrossRef] Adel, A.; Akbar, H. Autonomously implemented versatile path planning for mobile robots based on cellular automata and ant colony. Int. J. Comput. Intell. Syst. 2012, 5, 39–52. Saalfeld, A. Topologically consistent line simplification with the douglas-peucker algorithm. Cartogr. Geogr. Inf. Sci. 1999, 26, 7–18. [CrossRef] Stützle, T.; López-Ibánez, M.; Pellegrini, P.; Maur, M.; De Oca, M.M.; Birattari, M.; Dorigo, M. Parameter Adaptation in Ant Colony Optimization. Autonomous Search; Springer: Heidelberg, Germany, 2011. © 2016 by the authors; licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC-BY) license (http://creativecommons.org/licenses/by/4.0/).