An Air Route Network Planning Model of Logistics UAV Terminal Distribution in Urban Low Altitude Airspace

Abstract

Traditional terminal logistics distribution in urban areas is mainly concentrated on the ground, which leads to increasingly serious air pollution and traffic congestion. With the popularization of unmanned aerial vehicle (UAV) techniques and the reform of low altitude airspace, terminal logistics distribution is expected to be carried out by drones. Therefore, it is of great significance to construct a reasonable air route network for logistics UAV to ensure the safety and efficiency of operations. In this paper, a single route planning model and an air route network planning model for UAV were constructed by fully considering the complex urban low altitude environment, the flight performance of UAV and the characteristics of logistics tasks to regulate the flights of drones. Then, taking Jiangjun Road Campus of Nanjing University of Aeronautics and Astronautics as an example, the improved cellular automata (CA) was adopted to search for the optimal route between different waypoints, and the optimal spanning tree algorithm was used to construct the route network. The experimental results demonstrated that the improved CA could greatly reduce search time and obtain the optimal route while enhancing safety. With the satisfaction of the voyage, the needs of logistics and distribution constraints, a network that had smaller intersection points and redundancy was generated. The models and core ideas proposed in this paper can not only regulate operation of drones but also provide a solid foundation for the distribution of logistics UAV in the future.

1. Introduction

In the information age, UAV logistics, as a rapidly developing new strategic industry, has shown great virtues in both military and civilian fields. In the civilian sector, with the development of e-commerce and the aggravation of distribution demand, the traditional distribution mode of supply centers and couriers gradually has shown defects. Issues such as road congestion, air pollution and low distribution efficiency brought by the last mile community terminal distribution cannot be ignored. Therefore, UAV logistics has become a new distribution mode with broad prospects and vast audiences, driving the iterative update of logistics UAV technology.

UAV logistics belongs to the domain of Urban Air Mobility (UAM) [1]. It is defined by NASA as a safe and efficient urban air traffic system, which is a significant support for social and economic activities and the essence of future urban construction and development [2]. With the vigorous expansion of technologies such as autonomous driving, 5G communications, and electric propulsion [3], Uber [4], Airbus [5], NASA [6], Ehang [7] and other companies have proposed UAM concepts successively. In addition, they are expanded with Electric Vertical Take-off and Landing (eVTOL) as the main vehicle. With regard to the logistics field in UAM, Amazon [8], DHL [9], JD [10] and other enterprises have successively attempted to deliver parcels by drones in order to seize opportunities in this emerging field. In the outbreak of COVID-19 in 2020, the unique strengths of UAV, such as non-contact service and flexibility, forged its core competitiveness in the logistics field. RAND Corporation predicts that by 2030, drones will replace more than 20% of ground logistics distribution [11]. By then, it will not be surprising to see drones delivering packages over cities, and logistics drone operations will become normal.

Based on previous research, air route network planning can be divided into two categories: the local air route network planning and the global air route network planning. The former is the local optimization and/or adjustment of segments and nodes based on the optimization theory; the latter is to abandon the existing network and construct a new layout in a certain region [12]. At present, the air route network of UAV terminal distribution has not been formed, so it belongs to the latter. It can generally be transformed into a path search problem in state space, and a search algorithm is adopted to solve the route [13]. A wide array of methods have been used on this issue. In [14,15], the improved A* algorithm was established to solve the route planning model of low-altitude logistics UAV with multiple restriction conditions. The particle swarm optimization algorithm was proposed for the first time in [16]. Then, multi-intelligent particle filter on this basis was created to deal the path planning issue in unknown environment in [17]. Based on the improved ant colony algorithm, the route of single drone was constructed in [18,19]. It shows that the improved algorithm can greatly improve search efficiency and reduce path redundancy. In [20], the same algorithm was adopted to plan routes for cluster drones. Another method called CA aims to settle the network node optimization problem of avoiding the three areas. In [21], three air routes in Beijing flight information area were taken as examples to verify the effectiveness of the CA. However, the previously mentioned researches only planning a single route, and there are few studies on the air route network. In addition, the influencing factors considered in the planning are relatively simple, and the constraints of the UAV’s performance and external operating environment are not fully considered.

Therefore, this paper proposes an air route network planning method based on improved CA and optimal spanning tree after considering the characteristics of urban low altitude environment and the operational limitations of logistics UAV. Firstly, the three-dimensional modeling of the planned region is carried out by using the grid method to construct the low altitude grid environment. Then, under the constraints of complex low altitude environment and UAV performance, an improved CA based on cost function is raised to search for the optimal route between target points, thus generating the optimal route set. After that, according to the regional logistics needs and operational features of service providers, the optimal spanning tree is adopted to select the routes from the set. Finally, the air route network for trunk and branch is constructed. Multiple indicators are selected to conduct comprehensive evaluation and analysis to ensure the rationality and scientificity of the network.

2. Models

With the support of detailed traffic flow data, the planning method for the air route network in high altitude has been very mature,. However, since China has not adopted UAV to carry out logistics transportation on a large scale, the research on the route network of urban low-altitude logistics UAV is still in the initial stage, and detailed flow data cannot be obtained. The environment at low altitude is more complex, with more obstacles and protection zones than at high altitude. Traditional route network planning is mainly aimed at high-altitude areas without obstacles. The main restrictions at high-altitude are bad weather, military restricted zones and so on. Low-altitude route planning has not only the basic factors of high-altitude areas, but also the avoidance of dense buildings, trees and other factors in low-altitude areas. Therefore, the traditional route network planning method is not suitable for the urban low-altitude logistics UAV route network planning. Without the support of traffic flow data, the route network planning was divided into two steps: single air route planning and air route network planning.

2.1. Problem Description

In a certain region, there are numerous demand points and supply centers belonging to different service providers. The demand points have distribution requirements for each service provider. They need to be within the service range of each supply center. To ensure the safe, stable and smooth operation of drones and reduce the difficulty of the network management, this paper constructs a unified air route network of terminal distribution. Through it, various supply centers can distribute cargoes and mails to demand points. This network has two layers, and each of them is at the same flight level. The upper layer is the distribution network; the lower layer is the return network. After completing the distribution task at the demand point, the drones of different service providers can only choose to return to their own service provider’ s supply center.

2.2. Airspace Grid

The grid method is adopted to model the airspace and divide it into grid cells. At the same flight level, the space where the drones perform logistics tasks is a rectangle with length $L$, width $W$ and ABCD as vertexes. It is divided into grid cells whose length and width are respectively $l_{\mathrm{g}}$ and $v_{\mathrm{g}}$. Then, the task space can be divided into $m\times n$ grid cells, where $m=\mathrm{int}\left(L/l_{\mathrm{g}}\right),n=\mathrm{int}\left(W/v_{\mathrm{g}}\right)$, $\mathrm{int}\left(\right)$ is a upward rounding function. The center of each grid cell is considered an alternative waypoint.

Multi-rotor UAVs that can take off and land vertically are often used in terminal distribution. In the airspace of UAV operation, there are special areas such as obstacles, prohibited areas and restricted areas, which need to be avoided when planning an air route. Given that, protection zones were set up to enhance the safety margin of the route. The scope of the protection zone is the smallest convex polygon formed by special area expanding $z/2$ outward, where $z$ is the wheelbase of the drone.

Based on the discriminant function $P\left(x,y\right)$, whether the grid cell with center point $\left(x,y\right)$ has protection zone is identified.

$$P\left(x,y\right)=\{\begin{array}{ccc}\hfill 0,& & \mathrm{there} \mathrm{is} \mathrm{no} \mathrm{protection} \mathrm{zone} \mathrm{in} \mathrm{grid} \mathrm{cell}\hfill \\ \hfill -1,& & \mathrm{otherwise}\hfill \end{array}$$

(1)

In the light of neighbor structures, cellular neighbors are divided into three forms, such as Von. Neumann neighborhood, Moore neighborhood and extended Moore neighborhood in Figure 1.

For grid cell $\left(x,y\right)$, the density of the protection zone in the neighbors of it is calculated by analysing $u$ grid cells adjacent to it. It is defined as the risk degree.

$$r\left(x,y\right)=\frac{1}{u}\left|{\displaystyle \sum _{k=1}^{u}P}\left(x_k,y_k\right)\right|$$

(2)

where, $u$ is the number of cellular neighbors.

2.3. Single Air Route Planning Model of Logistics UAV

The center of the grid cell is considered as a waypoint, and the route consists of a series of adjacent waypoints. Set the coordinates of the start point and the end point of the route as $\left(x_0,y_0\right)$ and $\left(x_n,y_n\right)$ respectively. A series of adjacent waypoints of the planning route are $\left(x_i,y_i\right),i=1,2,3\cdots n-1$.

2.3.1. Objective Functions

(1) Length of single route.

$$D={\displaystyle \sum _{i=1}^{n-1}\sqrt{{\left(x_{i+1}-x_i\right)}^2+{\left(y_{i+1}-y_i\right)}^2}}$$

(3)

(2) Risk degree of single route.

$$R={\displaystyle \sum _{i=1}^{n-1}r\left(x_i,y_i\right)}={\displaystyle \sum _{i=1}^{n-1}\left[\frac{1}{u_i}\left|{\displaystyle \sum _{k=1}^{u_i}{P}_i\left(x_k,y_k\right)}\right|\right]}$$

(4)

To sum up, the objective function of a single route can be described by the following formular:

$$\min Q={\alpha }_1\cdot D+{\alpha }_2\cdot \left(\beta R\right),\begin{array}{cc}& {\alpha }_1+{\alpha }_2=1\end{array}$$

(5)

where $Q$ is the cost of route, and ${\alpha }_1,{\alpha }_2$ are weight coefficients, and $\beta$ is scaling factor.

2.3.2. Constraint Conditions

(1) Turning angle. Due to the limitation of the UAV’s maximum turning radius, the turning angle of each single turn within the route needs to be restricted. If the coordinates of the three adjacent waypoints are $\left(x_{i-1},y_{i-1}\right),\left(x_i,y_i\right),\left(x_{i+1},y_{i+1}\right)$, the turning angle of the route can be obtained by the vector formula based on cosine law, and it should be satisfied with the following formula:

$$0\le {\phi }_i=\arccos \left[\frac{\left(x_i-x_{i-1}\right)\left(x_{i+1}-x_i\right)+\left(y_i-y_{i-1}\right)\left(y_{i+1}-y_i\right)}{\sqrt{{\left(x_i-x_{i-1}\right)}^2+{\left(y_i-y_{i-1}\right)}^2}\sqrt{{\left(x_{i+1}-x_i\right)}^2+{\left(y_{i+1}-y_i\right)}^2}}\right]\le {\phi }_{\max }$$

(6)

where ${\phi }_{\max }$ refers to the maximum turning angle of the route.

(2) Turning buffer distance. For the sake of a safe turn, the drone needs to have a buffer distance before turning. If the coordinates of two adjacent turning points are $\left(x_i,y_i\right),\left(x_j,y_j\right)$, the default is to fly in a straight line between the turning points.

$$\sqrt{{\left(x_i-x_j\right)}^2+{\left(y_i-y_j\right)}^2}>d_{\min }$$

(7)

where $d_{\min }$ means the minimum turning buffer distance.

2.4. Air Route Network Planning Model of Logistics UAV

After all the optimal routes of any two points in the region have been obtained through the algorithm in Section 2.3, the appropriate nodes and the corresponding optimal route between nodes would be selected to construct the air route network according to the characteristics and constraints of the logistics task to model the logistics UAV network in the region.

2.4.1. Related Conception

(1) Adjacent matrix. For any two nodes, $A_{ij}$ is used to express whether there is a direct route between $i$ and $j$. Direct route means a route without transit points from the start point to the end point.

$$A_{ij}=\{\begin{array}{l}\begin{array}{cc}0,& \begin{array}{rr}\hfill & \hfill \mathrm{there} \mathrm{is} \mathrm{no} \mathrm{direct} \mathrm{route} \mathrm{between} i \mathrm{and} j\end{array}\end{array}\hfill \\ \begin{array}{cc}1,& \begin{array}{cc}& \mathrm{otherwise}\end{array}\end{array}\hfill \end{array}$$

(8)

(2) Connection matrix. For any two nodes, $B_{ij}$ is applied to indicate whether there is at least one direct or transit route between $i$ and $j$. Transit route refers to the route from the start point to the end point after passing through the transit points.

$$B_{ij}=\{\begin{array}{l}\begin{array}{cc}0,& \begin{array}{cc}& \mathrm{there} \mathrm{is} \mathrm{no} \mathrm{connected} \mathrm{route} \mathrm{between} i \mathrm{and} j\end{array}\end{array}\hfill \\ \begin{array}{cc}1,& \begin{array}{cc}& \mathrm{otherwise}\end{array}\end{array}\hfill \end{array}$$

(9)

(3) Distance matrix. For any two nodes, $D_{ij}$ is utilized to show the distance of direct route between $i$ and $j$.

2.4.2. Objective Functions

There are $k$ different service providers in the region, and the total number of supply centers of the service provider $w$ is $e_w$, $w=1,2,3\cdots ,k$; $u_j^w$ is the $j$th supply center of $w$, $j\in [1,e_w]$; $s_i$ is the demand point, $i=1,2,3\cdots ,m$; and $N$ represents the set of all points in the region.

(1) Node importance. In order to reflect the positions of demand point and supply center in the air route network, the importance index $I$ is established.

$$\begin{array}{l}{I}_{s_i}={\eta }_1\frac{N_{s_i}}{\overline{N}}+{\eta }_2\frac{G_{s_i}}{\overline{G}}+{\eta }_3\frac{Mo{n}_{s_i}}{\overline{Mon}}\\ I_{u_j^w}={\eta }_1^{’}\frac{S{t}_{u_j^w}}{\overline{St}}+{\eta }_2^{’}\frac{S{q}_{u_j^w}}{\overline{Sq}}+{\eta }_3^{’}\frac{N{u}_{u_j^w}}{\overline{Nu}}\end{array}$$

(10)

where $N_{s_i},G_{s_i},Mo{n}_{s_i}$ express demand point’s population, the number of deliveries and average income respectively; $S{t}_{u_j^w},S{q}_{u_j^w},N{u}_{u_j^w}$ represent supply center’s employees, site area and daily delivery volume, respectively; $\overline{N},\overline{G},\overline{Mon},\overline{St},\overline{Sq},\overline{Nu}$ are means of the mentioned factors; and ${\eta }_i,{\eta }_i^{’}$ are weights.

(2) Segment weight. The segment weight $q$ between point $a$ and point $b$ is composed of segment gravity $F$ and distance $D$, which can be expressed as follows.

$$F_{ab}=M\frac{I_a{I}_b}{D_{ab}}{A}_{ab},\begin{array}{cc}& a,b\in N\end{array}$$

(11)

$$q_{ab}={\lambda }_1\frac{\mu }{F_{ab}}+{\lambda }_2{D}_{ab}{A}_{ab},\begin{array}{cc}& a,b\in N\end{array}$$

(12)

where ${\lambda }_1,{\lambda }_2$ are weight coefficients and $\mu$ is the scaling factor.

In conclusion, the objective function of the air route network can be described by the following formula:

$$\min Q^{’}=\min \left(\frac{1}{2}{\displaystyle \sum _{j=1}^m{\displaystyle \sum _{i=1}^m{q}_{s_i}{}_{s_j}+{\displaystyle \sum _{i=1}^m{\displaystyle \sum _{w=1}^k{\displaystyle \sum _{j=1}^{e_w}{q}_{s_i}{}_{u_j^w}}}}+\frac{1}{2}{\displaystyle \sum _{w=1}^k{\displaystyle \sum _{i=1}^{e_w}{\displaystyle \sum _{j=1}^{e_w}{q}_{u_i^w}{}_{u_j^w}}}}}}\right),\begin{array}{cc}& i\ne j\end{array}$$

(13)

where $Q^{’}$ is the cost of the network.

2.4.3. Constraint Conditions

(1) Voyage distance. The route length between demand point and supply center ought to be within the maximum range of UAV $L_{\max }$

$$\begin{array}{l}2D_{s_i,u_j^w}\le L_{\max }\\ 2D_{s_i,s_j}\le L_{\max }\begin{array}{cc},& i\ne j\end{array}\\ \begin{array}{cc}2D_{u_i^w,u_j^w}\le L_{\max },& i\ne j\end{array}\end{array}$$

(14)

For a route with transit, if the drone starts from the supply center $u$, and passes through $n$ transit points $x$ before arriving at the demand point. Then, it goes straight from the demand point to the supply center as shown in Figure 2. The total voyage distance needs to meet the following formula:

$$D_{u,x_1}+{\displaystyle \sum _{i=1}^{n-1}{D}_{x_i,x_{i+1}}}+D_{x_n,s}+D_{s,u}^{’}\le L_{\max }$$

(15)

where $D_{s,u}^{’}$ represents the distance of return route.

(2) Demand point. Under the premise that there is no isolated demand point in the network, arbitrary demand point $s_i$ needs to be delivered by all service providers.

$${\displaystyle \prod _{w=1}^k{\displaystyle \sum _{j=1}^{e_w}{B}_{s_i,u_j^w}}}\ge 1$$

(16)

(3) Supply center. Under the premise that there is no isolated supply center in the network, arbitrary supply center $u_j^w$ ought to be linked with at least one demand point.

$${\displaystyle \sum _{i=1}^m{A}_{s_i,u_j^w}}\ge 1$$

(17)

(4) Service provider. Based on the requirements of logistics operation, there is no direct route between supply centers of different service providers.

$${\displaystyle \sum _{a=1}^{w-1}{\displaystyle \sum _{j=1}^{e_w}{\displaystyle \sum _{b=1}^{e_a}{A}_{u_j^w,u_b^a}}}}+{\displaystyle \sum _{c=w+1}^k{\displaystyle \sum _{j=1}^{e_w}{\displaystyle \sum _{b=1}^{e_c}{A}_{u_j^w,u_d^c}}}}=0$$

(18)

3. Algorithms

3.1. Single Air Route Generation Based on the Improved CA

3.1.1. Basic Principle

CA was first created by J. von Neumann, the father of computers, in order to simulate the self-replication function of living systems. It is composed of cells, space, state, neighbors and evolution rule [22], which can be expressed as the following quaternion:

$$A\left(C_q,S,N,R\right)$$

(19)

where $A$ represents CA system, $C_q$ expresses $q$-dimentional space, $S$ indicates finite discrete state of cell, $N$ is the set of cellular neighbors, and $R$ is the evolution rule.

3.1.2. Parameters

In the layered air route network, the routes at the same flight level are regarded as a two-dimensional plane network; thus, $q=2$. Based on the low altitude grid environment, the Moore neighborhood is adopted, so $u=8$. Then, the start point cell state value is 2, the end point cell state value is 3, the protection zone cell state value is −1, and the passable cell state value is 0.

3.1.3. Cost Function

In traditional CA, the route planning only depends on the step values. However, there are many points with the same step values, so the traditional CA randomly selects one of the points with the same value, which leads to winding, rather than straight routes. To overcome the defects of traditional CA, the cost function is introduced as a signal for route searching, so a better next step can be quickly found as long as all the costs of the potential next points are calculated. The cost function can optimize the route between start point and end point, and it can be expressed as follows:

$$P_{now}=d_{start,now}+d_{now,end}+\beta r_{now}$$

(20)

where $d_{start,now}$ indicates the Euclidean distance from the start point to the current cell, $d_{now,end}$ represents the Euclidean distance from the current cell to the end point, $r_{now}$ expresses the risk degree of current cell, and $\beta$ is the scaling factor.

3.1.4. Pseudo-Code

The first stage of this algorithm was to calculate all the step values from start point to end point. The matrix cell_status was to record the state value of each point. Neighbour_status takes the mean of the state value set of the current cell’s neighbours. The status of start point and end point were set as 2 and 3 respectively. Geo_info was the matrix which recorded the information of obstacle or non-obstacle, if this point was obstacle, the value of geo_info[i][j] equaled −1 and it would be 0 for point without obstacle.

All the values of cell_status except start point, end point and obstacles minus 2 are the step values of each cell from start point to themselves.

The second stage is to find the best route from end point to start point. The route is searched in descending order of cellular steps, which means the search orientation of the grid cell. A cell in the neighborhood of the end point was selected, which had the smallest state value among all neighbors. Then according to Equation (20), the costs p of those cells were obtained, and the operations shown in pseudo code were performed. Step[i][j] represents the step values of current cell to the start point. Neighbour_step indicates the step values of current cell’s neighbour. Neighbour_p meaned the costs of neighbour cells. Route_point was the matrix of the selected waypoint.

3.2. Air Route Network Planning Based on the Optimal Spanning Tree

3.2.1. Basic Principle

The optimal spanning tree algorithm was adopted to generate a connected tree without cycles among all nodes, which minimizes the sum of weights of all segments. All points are set as $a_i,i=1,2,3,\cdots ,N$, including the supply center and the demand point. Based on the Equation (12), the weight of each branch is defined as follows:

$$w_{ij}=\{\begin{array}{ll}{\lambda }_1\frac{\mu }{F_{ij}}+{\lambda }_2{D}_{ij},\hspace{1em}\hfill & \mathrm{there} \mathrm{is} \mathrm{a} \mathrm{connected} \mathrm{route} \mathrm{between} \mathrm{point} i \mathrm{and} \mathrm{point} j\\ 0,\hfill & \mathrm{otherwise}\hfill \end{array}$$

(21)

3.2.2. Algorithm Flow

On the grounds of above analyses, the processes can be described as follows:

Step 1. A $N\times N$ route weight matrix $W$ is constructed.

$$\begin{array}{c}\begin{array}{ccc}\begin{array}{cc}& \end{array}{w}_1& \cdots & w{}_N\end{array}\\ \begin{array}{c}w{}_1\\ ⋮\\ w{}_N\end{array}\left\{\begin{array}{ccc}w{}_{11}& \cdots & w{}_{1N}\\ ⋮& \ddots & ⋮\\ w{}_{N1}& \cdots & w{}_{NN}\end{array}\right\}\end{array}$$

(22)

Step 2. Random line k in the matrix is selected, and it is recorded in the branch list. That means the tree grows from this node. Then, the corresponding column of the line is deleted to avoid forming a circle.

Step 3. The smallest element $w_{ij}$ in the branch list is chosen, and it is put in the route list. Then, row i is placed in the branch list and column j is deleted. The elements and columns in the route list and the branch list are not chosen again.

Step 4. Examine whether all columns have been added to the branch list. If all of them are added, the algorithm terminates, and the adjustment and expansion stage is entered. Otherwise, return to Step 3.

Step 5. Check whether all demand points and supply centers meet (18) and (19). If the constraints are met, proceed to the next step; otherwise, after connecting the demand points that do not satisfy the constraints with the adjacent need points that can satisfy (16) and (17), put them in the route list and proceed to the next step.

Step 6. For a number of demand points with greater importance, in order to meet their distribution needs, they are directly connected with the supply centers. Then, these routes are added into the route list as a branch so that the network is gradually transformed from tree to network.

Step 7. The repetitive routes from the route list are deleted to get the final network structure.

4. Simulation and Analysis

4.1. Simulation Environment

To prove the validity of the method proposed in this paper, Jiangjun Road Campus of Nanjing University of Aeronautics and Astronautics was selected as the environment of the air route network which can be seen in Figure 3. According to the actual situation, there were three supply centers belonging to two different service providers in this scenario. The demand points were fixed and concentrated in student dormitories and office buildings. Based on relevant data collection and analysis, experimental parameters and node information are shown in

Table 1. Parameter setting for air route network.

Parameter	Value	Parameter	Value
Length of grid cell $l_{\mathrm{g}}$	5 m	Width of grid cell $v_{\mathrm{g}}$	5 m
Planned area	1.0865 km2	Supply centers’ number	3
Distance weight ${\alpha }_1$	0.6	Demand points’ number $m$	25
Risk degree weight ${\alpha }_2$	0.4	Weight of segment gravity ${\lambda }_1$	0.5
Scaling factor of risk degree $\beta$	100	Weight of segment distance ${\lambda }_2$	0.5
Maximum turning Angle ${\phi }_{\max }$	$\pi /2$	Scaling factor of segment $\mu$	100
Turning buffer distance $d_{\min }$	5 m	Maximum range of UAV $L_{\max }$	2800 m
Service providers’ number $k$	2	UAV speed	5 m/s

and

Table 2. Index of demand point (D) and supply center (S).

No.	Type	Coordinate	Node Importance	No.	Type	Coordinate	Node Importance
1	D	(14, 7)	0.059638	15	D	(26, 76)	0.036818
2	D	(13, 27)	0.0601	16	D	(37, 125)	0.115057
3	D	(35, 21)	0.024546	17	D	(18, 125)	0.050571
4	D	(32, 36)	0.012273	18	D	(29, 162)	0.050571
5	D	(37, 8)	0.012273	19	D	(29, 197)	0.019176
6	S of A	(92, 39)	0.7	20	D	(56, 191)	0.019176
7	D	(129, 24)	0.024546	21	D	(89, 193)	0.019176
8	D	(145, 25)	0.024546	22	D	(125, 193)	0.019176
9	D	(172, 28)	0.024546	23	D	(147, 193)	0.019176
10	D	(179, 55)	0.024546	24	D	(170, 189)	0.019176
11	D	(86, 77)	0.024546	25	D	(169, 149)	0.019176
12	D	(8, 49)	0.02654	26	D	(195, 148)	0.115057
13	D	(8, 61)	0.02654	27	S of A	(113, 139)	0.3
14	D	(8, 76)	0.024546	28	S of B	(113, 125)	0.5

Notes: A denotes the first service provider, and B represents the second service provider.

respectively. All of the elevation data was obtained by field measurements. After airspace rasterization, a matrix of 215 by 205 was obtained. The altitude of logistics network was assumed at 50 m and the return layer was set at 30 m. The height of obstacles was mainly concentrated in the range of 30 to 50 m and the highest was 112 m.

4.2. Performance of Algorithm

According to the above descriptions, simulation was conducted. The contrasts of air route between improved CA and traditional CA can be seen in

Table 3. Comparison between traditional CA and improved CA.

Index	Traditional CA	Improved CA
Total time consuming/s	1826	1410.1
Average time consuming/s	4.7	3.5
Average number of search steps	135	108
Average single route length/m	716	625.3
Average flight time/s	143.2	125.1
Average number of turns	17	6

, Figure 4 and Figure 5, comprehensive factors of the air route network based on improved algorithm are shown in

Table 4. Air route network evaluation index.

Index	Optimized Network
Total length $D/m$	7233.3
Nonlinear coefficient $\xi$	1.25
Segment number $M$	30
Intersection number	5
Network connectivity $\frac{2M}{N}$	2.14
Network density $\frac{M}{S}$	6.657

. Based on the improved CA algorithm, the route could effectively avoid obstacles and reduced the risk degree, especially in the area with dense obstacles, which effectively guarantees the flight safety. From the perspective of configuration, it was apparent that the route was relatively straight compared with the route solved by traditional CA, which decreased unnecessary turns and enabled the UAV to operate more quickly. The average distance of single route was 20% less than the traditional CA so that the flight time per route decreased by 18.2 s. Furthermore, through orientation of cellular step values and localization of cost, the improved algorithm was more intelligent in the search process. The average single path search time was only 3.5 s, which was 1.2 s quicker than that of the traditional search time.

In this case, the 28 waypoints in the region constituted 378 origin-destination (OD) pairs. The route set was established according to the optimal routes which found by the improved CA among all OD pairs. As shown in Figure 6, the optimal spanning tree algorithm was used to select 30 routes from it to build the distribution network, and the routes of the UAV directly returning to each supply center after completing the distribution task at the demand point constituted the return route network. Therefore, the double layer air route network was composed of logistics route network and direct return route, which is shown in Figure 7.

In the optimal spanning tree algorithm, from the view of the rationality of the distribution route network, in the absence of traffic flow data, the idea of gravity was innovatively introduced, which could effectively connect the important nodes and reduce the waste of resources. At the same time, the distances of all routes were within the range of 2800 m of the UAV. From the perspective of route configuration, the optimal spanning tree algorithm could reduce the intersection points, and the total number of intersection points was as low as 5, which reduced the conflict among different routes. The non-linear coefficient was 1.25, indicating that there was an approximate linear connection between waypoints and transportation was relatively convenient in this network. In the view of logistics distribution, the route network could not only meet the distribution requirements of service providers but also meet the service needs of customers. In terms of the route network connectivity, the global network connectivity reached 2.14, denoting that the network connectivity was mature and capable of rapid distribution.

5. Conclusions

An air route network planning algorithm based on the improved CA and the optimal spanning tree for logistics UAV in urban low altitude environment was proposed in this paper. The cost function was introduced in CA, which could effectively reduce route searching time, decrease unnecessary turns, and greatly shorten route length. In addition, the safety factor of the route was effectively improved, including the area of dense obstacles, which could also effectively avoid obstacles and would not be trapped in a local search. Furthermore, by introducing gravity in optimal spanning tree, the important nodes could be connected as far as possible to ensure mature network connectivity and convenience. Based on this air route network planning method, a network with less redundancy, minor airway intersections, and reasonable structure could be established. With the large-scale operation of drones, the probability of conflicts between different flights increases. In the future, we will further consider conflict resolution between drones. Another promising direction is to research the dynamic route networks to guarantee the safe and stable operation of drones.