# Improved F-RRT∗ Algorithm for Flight-Path Optimization in Hazardous Weather

Authors: Xue Qiu; Yaohui Li; Rui Jin; Zhi Zhao; Jiajun Li; Donglin Lu; Linhui Ma
Journal: International Journal of Aerospace Engineering
 (2022)

---

## Abstract

Hazardous weather has become a major cause of flight delays in recent years. With the development of satellite navigation systems, the study of flight-path optimization under hazardous weather conditions has become especially important. In this study, radar data were used as the basis for the initial flight-restricted area under hazardous weather conditions, and the Graham algorithm was used to delineate the dynamic flight-restricted area by comprehensively considering the hazardous weather boundary changes along with the speed and direction. Then, under the grid environment model, the range of influence, size, and distribution characteristics of the flight-restricted area was examined, and the path optimization model was created according to constraints related to the path distance, corner size, and number of turning points. An improved F-RRT∗ algorithm was developed to solve the model. The algorithm can overcome the problems of traditional path planning algorithms, such as strong randomness, poor guidance, slow convergence speed, unsmooth paths, and poor tracing smoothness. Finally, a simulation analysis was conducted on the Guiyang–Guangzhou route in China as an example. This study can address the drawbacks of existing research on route change and provide sufficient theoretical support and reference for the implementation of specific route change plans in the future.


---

## Body

## 1. Introduction

Hazardous weather, which primarily refers to small and medium-scale weather elements such as thunderstorms, squall lines, icing, and wind shear [1], has become a major threat to air transport safety and causes flight delays. In recent years, it has accounted for approximately 50% of the total flight accidents. It is necessary to identify hazardous weather accurately and provide the best route during aircraft operation, which is of great significance to ensure the safety of aircraft navigation [2].At present, diversion models are generally used to determine the optimal route under hazardous weather, and several studies have been reported. Ahlstrom and Jaggard evaluated the implications of hazardous weather on flight operations and proposed a method to identify dangerous weather during flight [3]. Liu et al. comprehensively studied hazardous weather, FFA, aircraft, and the flight coordination of aircraft and established an aircraft diversion decision model in extreme weather [4]. Chen et al. studied the coverage path planning problem of autonomous heterogeneous unmanned aerial vehicles (UAVs) over a limited number of areas and realized a flight path covering all the areas of interest [5]. Ding and Sui comprehensively analyzed the damage of route-based systems on meteorological data, identified the flight segments that pilots should avoid in route selection, and then used an A∗ algorithm to divert the route [6]. Minyang et al., to solve the path planning problem of multi-degree-of-freedom robots in complex environments, proposed an RRT∗ algorithm that can address the frequent demand for UAVs [7]. Song et al. introduced several representative methods of path planning in the field of UAVs and provided several parameters for the 3D flight path planning of UAVs [8]. Chen et al. focused on the coverage path planning problem of heterogeneous UAVs and proposed an ant-colony-system-based algorithm to select appropriate UAV paths and effectively cover all areas [9].However, most current studies focus on the diversion of UAVs [10] and a few on the diversion of flights in hazardous weather, as found that most researchers focused on the route planning of flight diversions [11] with not fully consider the implications of dangerous weather levels and movement trends on the routes. To address these limitations, this study delineates and predicts dangerous weather areas, utilizes these areas as constraints, evaluates other limiting factors such as aircraft performance and flight safety, and establishes a diversion path. Further, a mathematical model of route planning is realized. In addition, this study proposes an improved F-RRT∗ algorithm based on the goal-biased strategy [12]. Most researchers used traditional path planning algorithms, such as A∗ [13], genetic [14], simulated annealing [15], and ant colony optimization algorithms, to realize flight diversion models [16]. However, the computational complexity of these algorithms is exponentially related to the degrees of freedom of the aircraft. Therefore, based on a previous study [17], this study proposes an improved F-RRT∗ algorithm based on a target bias strategy combined with the angle limit. It has demonstrated the ability to determine safe trajectories between two states in high-dimensional problems using system dynamics and constraints. Unlike other algorithms, the RRT∗ algorithm does not require an initial guess because it can determine a feasible path in a few iterations and then optimize the path in subsequent steps [18]. Initially, the traditional RRT∗ algorithm was used for robot path planning. Recently, it has been employed to study UAVs [19]. However, it has not been successfully applied for flight path planning. To address this limitation, this study modifies an RRT∗ algorithm according to the characteristics of the flight and obtains an F-RRT∗ algorithm, which is suitable for flight path optimization. This approach solves the problems of strong randomness, poor guidance, slow convergence, tortuous paths, and the poor tracking smoothness of traditional path planning algorithms [20]. It can more effectively solve the problem of flight path planning in dangerous weather.


## 2. Delineation of Flight-Restricted Areas in Hazardous Weather

Hazardous weather during flights mainly refers to hail, heavy rain, and thunderstorms. These hazardous weather conditions have adverse effects on flight safety and normal operation of the airspace system [21]. Reasonable and accurate delineation of the flight-restricted area is essential to reduce these adverse effects. In this study, we mainly used the monitored radar echo data to filter and divide the data in order to determine the initial range of the flight-restricted area [22]. The Graham algorithm was then used for the dynamic prediction of the flight-restricted area.


### 2.1. Initial Designation of Flight-Restricted Area in Hazardous Weather

In this study, the flight-restricted area was determined based on the radar echo map. First, the area ofradarweatherechointensity≥41dBZ was determined as the dangerous flight area according to the characteristics of heavy rainfall [23]. Then, image processing techniques were used to identify the hazardous weather areas in grayscale. Finally, the outline of the flight-restricted area was circled using convex polygons. The steps to delineate the flight-restricted area are as follows.
(1)
Obtain the relevant hazardous weather data, extract the areas withradarweatherechointensity≥41dBZ, and calculate the radar reflectance. Each radar reflectance corresponds to an RGB color value, and the area constituted by the pixel points combining these color values is the hazardous weather area(2)
Binarize the images of severe weather areas, thus highlighting them and making the background grayer(3)
Delineate the initial outline of the severe weather area. The binary image obtained in the second step is further recognized to improve the rough outline of the severe weather area(4)
Approximate the contour lines using irregular polygons. In this step, the resulting contour set is approximated using irregular polygons of a given accuracy, and a simple polygonal contour set is created after approximation(5)
Create a convex polygon. The convex packets of each simple polygon are located; the area enclosed by these convex polygons is the flight-restricted area. Subsequently, the delineated danger zone is processed(a)
Denoising process: according to the flight performance database, a danger zone of less than seven miles can be neglected when the aircraft is flying at a speed of no greater than 420 kt(b)
Combined processing: according to civil aviation regulations, the minimum separation between clouds should be at least 20 km for the aircraft to pass through them; clouds that do not meet the above condition are considered as a whole [24]


### 2.2. Dynamic Prediction of Flight-Restricted Area in Hazardous Weather

To perform the dynamic prediction of the flight-restricted area under hazardous weather, we first need to determine the initial convex polygon formed by the hazardous weather boundary area. Then, we introduce the gray prediction model, combine the data before the hazardous weather boundary point, and use theGraham algorithm to dynamically predict the changes in the hazardous weather boundary point within a short period of time in the future. Finally, the predicted area is extended outward by 25 km to determine the extent of the flight-restricted area at a certain time in the future.The algorithm first scans the planar point set to find the extreme points in the leftmost, rightmost, topmost, and bottommost directions to construct an initial convex packet and delete all the points inside the initial convex packet. Then, it groups the remaining point sets and generates a new convex packet for each group by using theGraham algorithm. Finally, the vertices of all subsets of the convex packets are considered as a new point set, and the Graham algorithm is applied again to generate the final convex packet.Based on the above analysis, theGraham algorithm is divided into the following three steps. (1) The individual reference points on the graph are determined. (2) The polar angles are extended outward from each reference point, and the resulting polar angles are sorted to form a simple polygon. (3) The concave points on the simple polygon are removed, and the individual convex points are then joined to form a convex polygon. The time complexity of the convex packet algorithm is largely determined by that of the sorting algorithm. Because the time complexity of (1) and (3) is On, that of (2) is Onlogn. Therefore, the total time complexity is Onlogn. At this point, the speed of computing the convex packets is relatively slow. Therefore, if the time complexity of the sorting algorithm can be reduced, the speed of computing the convex packets can be improved significantly [25]. The process is illustrated in Figure 1 below.Figure 1
Construction of a convex packet using theGraham algorithm.


### 2.3. Movement Prediction Model of Flight-Restricted Area

Previous studies proposed dynamic predictions for light-prohibited areas under the assumption that the moving speed and direction of extreme weather are known. However, in practical applications, these parameters cannot easily be determined. Moreover, the relationship between the movement of a flight-restricted area along thex and y directions and time can be studied based only on the constantly updated meteorological radar data. Subsequently, the movement prediction model of this area can be established.As the update frequency of meteorological radar data used in this study is 5 min, the starting timetn of the forecast is used as the benchmark. By converting the forecast time tii=1,2,3,⋯ into the forecast time slice TiTi=1,2,3,⋯, the conversion formula can be expressed as
(1)Ti=ti−tn5.


#### 2.3.1. Prediction Model of Uniform Movement

The motion of a thunderstorm can be modeled as a uniform motion along a certain direction in a short time as its component speeds along thex and y directions are unchanged. Therefore, the centroid positions of a thunderstorm, xi and yi, correspond to n historical time slices TiTi=1,2,3,⋯. Using the relationships among the speed, distance, and time, the moving speeds vx and vy along the x and y directions, respectively, are obtained as follows:
(2)vx=∑i=0n−1xi+1−xin−1,vy=∑i=0n−1yi+1−yin−1.According to the moving speeds,vx and vy, and the thunderstorm positions, xn and yn, corresponding to the n time slices Tn, we can predict the position of the thunderstorm after ii=1,2,3,⋯n time slices as follows:
(3)xn+i=xn+vx∗i,yn+i=yn+vy∗i.


#### 2.3.2. Prediction Model of Uniform Speed Movement

In Section2.3.1, the thunderstorm motion is regarded as a uniform motion. However, in summer, thunderstorms exhibit uniform acceleration or uniform deceleration during the development and dissipation phases. Therefore, according to the centroid positions of the thunderstorm, xi and yi, corresponding to n historical time slices TiTi=1,2,3,⋯n and using the relationship between acceleration, speed, and time, the accelerations ax and ay along the x and y directions, respectively, are obtained as follows:
(4)ax=∑i=0n−2xi+2−xi+1−xi+1−xin−2,(5)ay=∑i=0n−2yi+2−yi+1−yi+1−yin−2.In Equations (4) and (5), xi+2−xi+1−xi+1−xiandyi+2−yi+1−yi+1−yi represent the velocity change rates of two adjacent time slices. By calculating the average accelerations ax and ay and using the thunderstorm positions xn and yn corresponding to the time slice Tn, the position of the thunderstorm after ii=1,2,3,⋯n time slices can be expressed as follows:
(6)xn+i=xn+xn−xn−1∗i+12axi2,yn+i=yn+yn−yn−1∗i+12ayi2.


## 2.1. Initial Designation of Flight-Restricted Area in Hazardous Weather

In this study, the flight-restricted area was determined based on the radar echo map. First, the area ofradarweatherechointensity≥41dBZ was determined as the dangerous flight area according to the characteristics of heavy rainfall [23]. Then, image processing techniques were used to identify the hazardous weather areas in grayscale. Finally, the outline of the flight-restricted area was circled using convex polygons. The steps to delineate the flight-restricted area are as follows.
(1)
Obtain the relevant hazardous weather data, extract the areas withradarweatherechointensity≥41dBZ, and calculate the radar reflectance. Each radar reflectance corresponds to an RGB color value, and the area constituted by the pixel points combining these color values is the hazardous weather area(2)
Binarize the images of severe weather areas, thus highlighting them and making the background grayer(3)
Delineate the initial outline of the severe weather area. The binary image obtained in the second step is further recognized to improve the rough outline of the severe weather area(4)
Approximate the contour lines using irregular polygons. In this step, the resulting contour set is approximated using irregular polygons of a given accuracy, and a simple polygonal contour set is created after approximation(5)
Create a convex polygon. The convex packets of each simple polygon are located; the area enclosed by these convex polygons is the flight-restricted area. Subsequently, the delineated danger zone is processed(a)
Denoising process: according to the flight performance database, a danger zone of less than seven miles can be neglected when the aircraft is flying at a speed of no greater than 420 kt(b)
Combined processing: according to civil aviation regulations, the minimum separation between clouds should be at least 20 km for the aircraft to pass through them; clouds that do not meet the above condition are considered as a whole [24]


## 2.2. Dynamic Prediction of Flight-Restricted Area in Hazardous Weather

To perform the dynamic prediction of the flight-restricted area under hazardous weather, we first need to determine the initial convex polygon formed by the hazardous weather boundary area. Then, we introduce the gray prediction model, combine the data before the hazardous weather boundary point, and use theGraham algorithm to dynamically predict the changes in the hazardous weather boundary point within a short period of time in the future. Finally, the predicted area is extended outward by 25 km to determine the extent of the flight-restricted area at a certain time in the future.The algorithm first scans the planar point set to find the extreme points in the leftmost, rightmost, topmost, and bottommost directions to construct an initial convex packet and delete all the points inside the initial convex packet. Then, it groups the remaining point sets and generates a new convex packet for each group by using theGraham algorithm. Finally, the vertices of all subsets of the convex packets are considered as a new point set, and the Graham algorithm is applied again to generate the final convex packet.Based on the above analysis, theGraham algorithm is divided into the following three steps. (1) The individual reference points on the graph are determined. (2) The polar angles are extended outward from each reference point, and the resulting polar angles are sorted to form a simple polygon. (3) The concave points on the simple polygon are removed, and the individual convex points are then joined to form a convex polygon. The time complexity of the convex packet algorithm is largely determined by that of the sorting algorithm. Because the time complexity of (1) and (3) is On, that of (2) is Onlogn. Therefore, the total time complexity is Onlogn. At this point, the speed of computing the convex packets is relatively slow. Therefore, if the time complexity of the sorting algorithm can be reduced, the speed of computing the convex packets can be improved significantly [25]. The process is illustrated in Figure 1 below.Figure 1
Construction of a convex packet using theGraham algorithm.


## 2.3. Movement Prediction Model of Flight-Restricted Area

Previous studies proposed dynamic predictions for light-prohibited areas under the assumption that the moving speed and direction of extreme weather are known. However, in practical applications, these parameters cannot easily be determined. Moreover, the relationship between the movement of a flight-restricted area along thex and y directions and time can be studied based only on the constantly updated meteorological radar data. Subsequently, the movement prediction model of this area can be established.As the update frequency of meteorological radar data used in this study is 5 min, the starting timetn of the forecast is used as the benchmark. By converting the forecast time tii=1,2,3,⋯ into the forecast time slice TiTi=1,2,3,⋯, the conversion formula can be expressed as
(1)Ti=ti−tn5.


### 2.3.1. Prediction Model of Uniform Movement

The motion of a thunderstorm can be modeled as a uniform motion along a certain direction in a short time as its component speeds along thex and y directions are unchanged. Therefore, the centroid positions of a thunderstorm, xi and yi, correspond to n historical time slices TiTi=1,2,3,⋯. Using the relationships among the speed, distance, and time, the moving speeds vx and vy along the x and y directions, respectively, are obtained as follows:
(2)vx=∑i=0n−1xi+1−xin−1,vy=∑i=0n−1yi+1−yin−1.According to the moving speeds,vx and vy, and the thunderstorm positions, xn and yn, corresponding to the n time slices Tn, we can predict the position of the thunderstorm after ii=1,2,3,⋯n time slices as follows:
(3)xn+i=xn+vx∗i,yn+i=yn+vy∗i.


### 2.3.2. Prediction Model of Uniform Speed Movement

In Section2.3.1, the thunderstorm motion is regarded as a uniform motion. However, in summer, thunderstorms exhibit uniform acceleration or uniform deceleration during the development and dissipation phases. Therefore, according to the centroid positions of the thunderstorm, xi and yi, corresponding to n historical time slices TiTi=1,2,3,⋯n and using the relationship between acceleration, speed, and time, the accelerations ax and ay along the x and y directions, respectively, are obtained as follows:
(4)ax=∑i=0n−2xi+2−xi+1−xi+1−xin−2,(5)ay=∑i=0n−2yi+2−yi+1−yi+1−yin−2.In Equations (4) and (5), xi+2−xi+1−xi+1−xiandyi+2−yi+1−yi+1−yi represent the velocity change rates of two adjacent time slices. By calculating the average accelerations ax and ay and using the thunderstorm positions xn and yn corresponding to the time slice Tn, the position of the thunderstorm after ii=1,2,3,⋯n time slices can be expressed as follows:
(6)xn+i=xn+xn−xn−1∗i+12axi2,yn+i=yn+yn−yn−1∗i+12ayi2.


## 2.3.1. Prediction Model of Uniform Movement

The motion of a thunderstorm can be modeled as a uniform motion along a certain direction in a short time as its component speeds along thex and y directions are unchanged. Therefore, the centroid positions of a thunderstorm, xi and yi, correspond to n historical time slices TiTi=1,2,3,⋯. Using the relationships among the speed, distance, and time, the moving speeds vx and vy along the x and y directions, respectively, are obtained as follows:
(2)vx=∑i=0n−1xi+1−xin−1,vy=∑i=0n−1yi+1−yin−1.According to the moving speeds,vx and vy, and the thunderstorm positions, xn and yn, corresponding to the n time slices Tn, we can predict the position of the thunderstorm after ii=1,2,3,⋯n time slices as follows:
(3)xn+i=xn+vx∗i,yn+i=yn+vy∗i.


## 2.3.2. Prediction Model of Uniform Speed Movement

In Section2.3.1, the thunderstorm motion is regarded as a uniform motion. However, in summer, thunderstorms exhibit uniform acceleration or uniform deceleration during the development and dissipation phases. Therefore, according to the centroid positions of the thunderstorm, xi and yi, corresponding to n historical time slices TiTi=1,2,3,⋯n and using the relationship between acceleration, speed, and time, the accelerations ax and ay along the x and y directions, respectively, are obtained as follows:
(4)ax=∑i=0n−2xi+2−xi+1−xi+1−xin−2,(5)ay=∑i=0n−2yi+2−yi+1−yi+1−yin−2.In Equations (4) and (5), xi+2−xi+1−xi+1−xiandyi+2−yi+1−yi+1−yi represent the velocity change rates of two adjacent time slices. By calculating the average accelerations ax and ay and using the thunderstorm positions xn and yn corresponding to the time slice Tn, the position of the thunderstorm after ii=1,2,3,⋯n time slices can be expressed as follows:
(6)xn+i=xn+xn−xn−1∗i+12axi2,yn+i=yn+yn−yn−1∗i+12ayi2.


## 3. Multiobjective Redirected Path Optimization Model

### 3.1. Three Types of Flight-Restricted Areas

Considering the shape, size, and distribution characteristics of the flight-restricted areas, this study classified them into three main categories—block-shaped, strip-shaped, and scattered. A block-shaped flight-restricted area is a convex polygon with the horizontal scale of the flight-restricted area usually larger than100×100km2; the affected flight area is usually large [26]. Both dynamic flight-restricted areas caused by block clouds and those formed over a period of time are categorized as block-shaped flight-restricted areas.There are generally two types of strip-shaped flight-restricted areas. The first is a single-strip flight-restricted area with a long and narrow convex polygon. The second is a series of small, banded flight-restricted areas with long and narrow strips distributed along both sides of the route. The strip-shaped flight restriction zone covers almost 2–3 points. The ratio between the long and short sides of the outer rectangle is greater than two.Scattered flight-restricted areas refer to multiple small-scale flight-restricted areas along the route on both sides of the arbitrary distribution, which have the characteristics of small scale, large number, and small impact range. The main difference between the scattered and block-shaped flight-restricted areas is that the latter is limited to a single flight-restricted area, whereas the former consists of multiple flight-restricted areas. The geometric scale of the latter is much smaller than that of the former. Unlike the strip-shaped distributed flight-restricted area, the scattered flight-restricted area is distributed among a small number of route points. The strip distribution is not obvious. Moreover, the distribution is irregular, and the scale of the flight-restricted area is small [27]. The distribution of banded flight-restricted areas covers almost the entire route, involving at least five waypoints with distinct banding characteristics. The dynamic flight-restricted areas caused by the scattered small-scale clouds and the flight-restricted areas formed during the fixed period belong to this category. The distribution of the scattered flight-restricted areas has no specific pattern, and the number of restricted areas may be large or small. Therefore, the grid approach has advantages of simplicity, intuitiveness, and ease of implementation.(Abridged here)


### 3.2. Description of Related Concepts

#### 3.2.1. Turning Points of the Change of Course

The number of path-optimized turning points determines the workload of the controllers and pilots when performing path planning. The number of turning points for redirecting an aircraft in clear weather is usually 2–4, and the maximum number of turns for an aircraft in hazardous weather is 7–8.


#### 3.2.2. Turning Angle of the Aircraft

The turning angle of an aircraft is the amount of constant change in its heading during flight. As shown in Figure2, θ indicates the angular magnitude of the turn. When radar control is used, the aircraft’s turn angle usually does not exceed 60°. The turning angle when flying in restricted areas under hazardous weather conditions usually does not exceed 90°.Figure 2
Turning angle of the aircraft.


#### 3.2.3. Segment Distance

The segment distance is the straight-line distance between two turning points. This distance should be long enough to ensure that an aircraft completes a turn smoothly [28]. Therefore, for a normal flight, the distance between two turning points must exceed 7.4km, and the turning radius must exceed 3.7km.


### 3.3. Modification of Navigation Model and Algorithm

This study focused on the path optimization of a single flight; therefore, in addition to avoiding the flight-restricted area, other factors such as path length, magnitude of turning angle, and number of turns [29] should be considered during the actual modeling process. Under actual circumstances, the turning angle of the flight around the flight-restricted area should not exceed the maximum angle allowed for the aircraft performance, and the turn radius should not be lower than the minimum radius of turn of the aircraft. In the path optimization process, the number of turns should be minimized to avoid aircraft congestion [30]. Therefore, this model adopts the shortest path, minimum turn angle, and minimum number of turns as the objectives, considering the above factors. Simultaneously, to achieve a unified metric for the mathematical modeling of multiple objectives, it was decided to perform dimensionless processing of multiple metrics. The model built based on the above requirements is as follows:Objective function was as follows:(7)minC=ω1Cl+ω2Ct+ω3Ca,(8)Cl=1−dps,pDdps,pl+∑i=1Nrdpi,pi+1+dpNr,pD,(9)Ct=NtNmax,(10)Ca=∑i=1NrcosθNt.Binding conditions were as follows:(11)ω1+ω2+ω3=1,(12)ω1≥0,ω2≥0,ω3≥0,(13)Dpi,pi+1Fi≥Dmin∀i=1⋯Nt,j=1⋯n,(14)dpi,pi+1≥dmin∀i=1⋯Nt,(15)θ<θmax,(16)Nt≤Nmax.Here are the following definition of parameters:Cl,Ct,Ca: path length cost, corner number cost, and turning angle change costω1,ω2,ω3: path length cost weighting, corner number cost weighting, and turning angle change cost weightingPs,Pi,PD: path starting point, turning point, and path ending point in the path optimization processDpi,pi+1Fi: minimum straight-line distance of the path between two turning points from the flight-restricted areadps,pD: initial distance between the starting and ending points of the change of coursedpi,pi+1: straight-line distance between two turning pointsDmin: minimum distance to deviate from the restricted area, usually 10kmdmin: shortest distance allowed between two turning points, usually 7.4kmNt: actual number of cornersNmax: maximum number of corners allowedθ: corner sizeθmax: maximum turning angle allowedFi: i-th flight-restricted areaAccording to the above formula of the objective function, the best path can be obtained when the minimum value of the objective function is considered [31]. Equations (7)–(10) normalize the costs for covering the segment distance, given number of turns, and angle change, and ensure that the smaller the value of Cl, the shorter is the path length after path optimization; the smaller the value of Ct, the smaller is the number of turns, and the smaller the value of Ca, the smaller is the angle change. Equation Equation (11) ensures that the sum of the weight coefficients is 1 and positive; Equation (12) ensures the safety of path optimization, and Equations (13)–(16) limit the segment distance, maximum turning angle, and number of turning points.


## 3.1. Three Types of Flight-Restricted Areas

Considering the shape, size, and distribution characteristics of the flight-restricted areas, this study classified them into three main categories—block-shaped, strip-shaped, and scattered. A block-shaped flight-restricted area is a convex polygon with the horizontal scale of the flight-restricted area usually larger than100×100km2; the affected flight area is usually large [26]. Both dynamic flight-restricted areas caused by block clouds and those formed over a period of time are categorized as block-shaped flight-restricted areas.There are generally two types of strip-shaped flight-restricted areas. The first is a single-strip flight-restricted area with a long and narrow convex polygon. The second is a series of small, banded flight-restricted areas with long and narrow strips distributed along both sides of the route. The strip-shaped flight restriction zone covers almost 2–3 points. The ratio between the long and short sides of the outer rectangle is greater than two.Scattered flight-restricted areas refer to multiple small-scale flight-restricted areas along the route on both sides of the arbitrary distribution, which have the characteristics of small scale, large number, and small impact range. The main difference between the scattered and block-shaped flight-restricted areas is that the latter is limited to a single flight-restricted area, whereas the former consists of multiple flight-restricted areas. The geometric scale of the latter is much smaller than that of the former. Unlike the strip-shaped distributed flight-restricted area, the scattered flight-restricted area is distributed among a small number of route points. The strip distribution is not obvious. Moreover, the distribution is irregular, and the scale of the flight-restricted area is small [27]. The distribution of banded flight-restricted areas covers almost the entire route, involving at least five waypoints with distinct banding characteristics. The dynamic flight-restricted areas caused by the scattered small-scale clouds and the flight-restricted areas formed during the fixed period belong to this category. The distribution of the scattered flight-restricted areas has no specific pattern, and the number of restricted areas may be large or small. Therefore, the grid approach has advantages of simplicity, intuitiveness, and ease of implementation.(Abridged here)


## 3.2. Description of Related Concepts

### 3.2.1. Turning Points of the Change of Course

The number of path-optimized turning points determines the workload of the controllers and pilots when performing path planning. The number of turning points for redirecting an aircraft in clear weather is usually 2–4, and the maximum number of turns for an aircraft in hazardous weather is 7–8.


### 3.2.2. Turning Angle of the Aircraft

The turning angle of an aircraft is the amount of constant change in its heading during flight. As shown in Figure2, θ indicates the angular magnitude of the turn. When radar control is used, the aircraft’s turn angle usually does not exceed 60°. The turning angle when flying in restricted areas under hazardous weather conditions usually does not exceed 90°.Figure 2
Turning angle of the aircraft.


### 3.2.3. Segment Distance

The segment distance is the straight-line distance between two turning points. This distance should be long enough to ensure that an aircraft completes a turn smoothly [28]. Therefore, for a normal flight, the distance between two turning points must exceed 7.4km, and the turning radius must exceed 3.7km.


## 3.2.1. Turning Points of the Change of Course

The number of path-optimized turning points determines the workload of the controllers and pilots when performing path planning. The number of turning points for redirecting an aircraft in clear weather is usually 2–4, and the maximum number of turns for an aircraft in hazardous weather is 7–8.


## 3.2.2. Turning Angle of the Aircraft

The turning angle of an aircraft is the amount of constant change in its heading during flight. As shown in Figure2, θ indicates the angular magnitude of the turn. When radar control is used, the aircraft’s turn angle usually does not exceed 60°. The turning angle when flying in restricted areas under hazardous weather conditions usually does not exceed 90°.Figure 2
Turning angle of the aircraft.


## 3.2.3. Segment Distance

The segment distance is the straight-line distance between two turning points. This distance should be long enough to ensure that an aircraft completes a turn smoothly [28]. Therefore, for a normal flight, the distance between two turning points must exceed 7.4km, and the turning radius must exceed 3.7km.


## 3.3. Modification of Navigation Model and Algorithm

This study focused on the path optimization of a single flight; therefore, in addition to avoiding the flight-restricted area, other factors such as path length, magnitude of turning angle, and number of turns [29] should be considered during the actual modeling process. Under actual circumstances, the turning angle of the flight around the flight-restricted area should not exceed the maximum angle allowed for the aircraft performance, and the turn radius should not be lower than the minimum radius of turn of the aircraft. In the path optimization process, the number of turns should be minimized to avoid aircraft congestion [30]. Therefore, this model adopts the shortest path, minimum turn angle, and minimum number of turns as the objectives, considering the above factors. Simultaneously, to achieve a unified metric for the mathematical modeling of multiple objectives, it was decided to perform dimensionless processing of multiple metrics. The model built based on the above requirements is as follows:Objective function was as follows:(7)minC=ω1Cl+ω2Ct+ω3Ca,(8)Cl=1−dps,pDdps,pl+∑i=1Nrdpi,pi+1+dpNr,pD,(9)Ct=NtNmax,(10)Ca=∑i=1NrcosθNt.Binding conditions were as follows:(11)ω1+ω2+ω3=1,(12)ω1≥0,ω2≥0,ω3≥0,(13)Dpi,pi+1Fi≥Dmin∀i=1⋯Nt,j=1⋯n,(14)dpi,pi+1≥dmin∀i=1⋯Nt,(15)θ<θmax,(16)Nt≤Nmax.Here are the following definition of parameters:Cl,Ct,Ca: path length cost, corner number cost, and turning angle change costω1,ω2,ω3: path length cost weighting, corner number cost weighting, and turning angle change cost weightingPs,Pi,PD: path starting point, turning point, and path ending point in the path optimization processDpi,pi+1Fi: minimum straight-line distance of the path between two turning points from the flight-restricted areadps,pD: initial distance between the starting and ending points of the change of coursedpi,pi+1: straight-line distance between two turning pointsDmin: minimum distance to deviate from the restricted area, usually 10kmdmin: shortest distance allowed between two turning points, usually 7.4kmNt: actual number of cornersNmax: maximum number of corners allowedθ: corner sizeθmax: maximum turning angle allowedFi: i-th flight-restricted areaAccording to the above formula of the objective function, the best path can be obtained when the minimum value of the objective function is considered [31]. Equations (7)–(10) normalize the costs for covering the segment distance, given number of turns, and angle change, and ensure that the smaller the value of Cl, the shorter is the path length after path optimization; the smaller the value of Ct, the smaller is the number of turns, and the smaller the value of Ca, the smaller is the angle change. Equation Equation (11) ensures that the sum of the weight coefficients is 1 and positive; Equation (12) ensures the safety of path optimization, and Equations (13)–(16) limit the segment distance, maximum turning angle, and number of turning points.


## 4. Model Solution Using Algorithm for Change of Navigation

### 4.1. RRT∗ Algorithm

Traditional path planning algorithms such asA∗, shortest path, and genetic algorithms, which require modeling the obstacles in a defined space, are not suitable for solving the path planning of a multidegree-of-freedom aircraft in a complex environment [31]. The fast extended random tree algorithm, that is, the RRT∗ algorithm, can solve path planning problems in high-dimensional spaces and under complex constraints by performing collision detection on sampled points in the state space. The intention is to generate random points and then search for a step forward from the starting point to the target point. In this manner, it can effectively avoid obstacles, prevent the path from falling into the local minima, and converge faster [32].According to Figure3, the steps of the RRT∗ algorithm are as follows.Figure 3
Expansion diagram of the RRT∗ algorithm.Initialize the random tree, adding nodes starting with an arbitrary random tree with onlyQinit.
(1)
Execute thesample function to obtain a random point Qnearest in the mapIterate through all the nodes in the tree and find the pointQnearest with the least cost to reach Qrand.
(2)
Execute theextend function to obtain the extension point Qnew of the specified length in the direction of Qnearest to Qrand and perform collision detection on Qnew. If false, Qnearest is designated as the parent node of Qnew, and the line connecting the two points is drawn(3)
Determine whetherQnew has reached the specified target range, if it has, end the loop; otherwise, continue to execute the loop until the target point Qrand is found(4)
To increase the practicality of the algorithm, a running time limit and a limit on the number of nodes searched are added to the above process. In other words, the algorithm is judged to have failed after searching for a finite number of nodes within the specified time and not finding the target point


### 4.2. Improved F-RRT∗ Algorithm

To solve the problems of strong randomness, poor orientation, slow convergence speed, and poor smoothness of path twisting and tracing of thetraditionalrapidlyexploringrandomtree∗algorithm, an improved RRT∗ algorithm based on a target-biasing strategy combined with angle restriction is proposed—the F−RRT∗ algorithm.First, a target bias threshold is set based on the traditional RRT∗ algorithm, and the randomness of tree growth is reduced by expanding the internal nodes of the tree through the target bias strategy [33]. Simultaneously, a local expansion mechanism is introduced to avoid the local optimum problem caused by changing the sampling structure. Then, a bidirectional tree growth strategy is used to improve the convergence speed of the algorithm. Finally, redundant nodes are removed, and cubic spline interpolation is used to increase the smoothness of the path.


#### 4.2.1. Extended Step

In the RRT∗ algorithm, Qinit is extended to Qnearest by a step to obtain the sampling point Qnew, expressed as
(17)Qnew=Qinit+φQnearest−QinitLQinit,Qnearest,where φ is the extended step size, and LQinit,Qnearest is the Euclidean distance between the points Qinit and Qnearest.


#### 4.2.2. Target Bias Strategy

The traditional RRT∗ algorithm has a large range of sampling points and is prone to randomness; therefore, the algorithm search time is longer [34]. By adopting the target bias strategy, the target points are used as sampling points for random tree expansion with a certain probability P, which reduces the number of sampling points and accelerates the efficiency of the random tree expansion [35].


#### 4.2.3. Metric Function

The objective of the metric function is to find the point closest to the selected pointQinit in the randomly expanded tree, namely, Qnearest. In the traditional RRT∗ algorithm, the metric function of the Euclidean distance is used [36]. The improved F-RRT∗ algorithm introduces an angle constraint while considering two-point spacing, which makes the entire path smoother.The distance and angle are unified such thatLd represents the distance function, and Mθ represents the declination function after planning. The node that minimizes the function Mxi,xj is near the node formed in the search extension tree, which can be used as Qnew, as shown in Figure 4.
(18)Ld=dmax−ddmax,Mθ=θmax−θθmax,Li,j=xj−xi,θi,j=θj−θi,Mxi,xj=LdLi,j+Mθθi,j,where the coordinates of the point Qinit are xi,yi, and those of Qnearest are xj,yj; dmax and θmax are the maximum distance and maximum angle of the search, respectively.Figure 4
Extended diagram of the F-RRT∗ algorithm.


#### 4.2.4. Path Pruning and Smoothing

The sampling points of the traditional RRT algorithm are highly random, and the generated paths often have numerous redundant nodes [37]. The efficiency in complex environments with many obstacles is low, making path pruning an essential step. The newly generated path nodes are filtered based on the traditional paths to remove redundant nodes, and the pruned paths are then simplified using a greedy algorithm [38].The specific pruning steps are as follows:Begin from the starting pointQinit and find all its Qnearest nodes. If the connection between Qinit and one of its Qnearest nodes does not contain obstacles, the Qnearest node is taken as the new node Qnew, and all other Qnearest nodes are deleted. Conversely, if the connection between Qinit and one of its Qnearest nodes encounters obstacles, then the above process is repeated starting from the parent node of the point Qinit.
(1)
Find the connection betweenQinit and the next Qnearest node until the Qrand node is added to the new path. Finally, the identified nodes are connected by a smooth curve [39]TheK-times B-spline curve equation is as follows:
(19)Pu=∑i=onBi,kuPii=1,2⋯n,where Pi is the node of the curve control, Bi,ku is the K-times B-spline basis function, and the node vector u=u0,u1,⋯,un+k+1 determines the K-times segmented curve [40].The basis functions of the tripleB-spline curve are as follows:
(20)b0=16−u3+3u2−3u+1,b1=163u3−6u2+4,b2=16−3u3+6u2+3u+1,b3=16u3.


### 4.3. Steps to Improve the RRT∗ Algorithm

The main difference between the improved RRT∗ algorithm and RRT∗ algorithm is that the former has two additional recomputation processes for the new node Qnew, which are as follows:
(1)
The process of reselecting the parent node forQnew contains an additional reselection process(2)
An additional rewiring random tree processBeginning from the starting pointQinit, we search for the “nearest neighbors” in the selected range; that is, we search for the suitable Qnearest node among all the Qnearest nodes from Qinit as an alternative to the parent node Qnew. The distance from the starting point Qinit to the nearest Qnearest node and the distance from the nearest Qnearest node to the parent node Qnew are calculated [41]. The specific steps are illustrated in Figure 5. Figure 5(a) shows the expansion process of the random tree at a certain moment, where node 0 represents the starting point Qinit, node 9 is the newly generated node Qnearest, node 6 is the parent node Qnew of node 9, and the generated path sequence number is 0-4-6-9 with a distance cost of 10+5+1=16. Figure 5(b) shows the process of identifying the parent node once again based on Figure 5(a). The newly generated Qnearest node, that is, node 9, is considered as the center, and its “near neighbors” are detected within the selected range [42], that is, nodes 5, 6, and 8; they are Qnearest composed of paths 0-1-5-9, 0-4-6-9, and 0-1-5-8-9, and the distances are 3+5+3=11, 10+5+1=16, and 3+5+1+3=12, respectively. It can be seen that node 5 is the closest to node 9; therefore, the parent node of node 9 is changed from the original node 6 to node 5, and the results are shown in Figure 5(b).Figure 5
Process of reselecting the parent node.
(a)(b)After reselecting the parent node, the rewiring process is performed on the random tree to further reduce the distance cost of the random extended tree [43]. The specific steps are illustrated in Figure 6. Figure 6(a) shows the rewiring process of the random tree at a certain moment, where node 0 is the starting point Qinit and node 9 is the Qnew node; that is, node 9 is the center, and its “nearest neighbors” are found in the selected range, which are nodes 5, 6, and 8. The distances from the starting point to each “nearest neighbor” node are 8, 15, and 9, respectively. Figure 6(b) shows the rewiring process based on Figure 6(a). Suppose that the parent node of node 5 is not node 1 but node 9, and the path is 0-4-9-5 with distance 17, which is greater than the original path distance 8, then the parent node of node 5 is not changed. Similarly, find the new parent node of node 6. Assuming that it is node 9 with path 0-1-5-9-6 and distance 12, which is less than the original path distance of 15, the parent node of node 6 is changed from node 4 to node 9. Finally, we find the new parent node of node 8. Assuming it is node 9, the path is 0-1-5-9-8, and the distance is 14, which is greater than the original path distance 9; therefore, the parent node of node 8 is not changed [44]. In summary, only the parent node of node 6 must be changed, and the newly generated paths are shown in Figure 6(b). From the above, it can be seen that the significance of rewiring is to reduce the distance cost of some newly generated node paths, thus reducing the distance of the entire optimized path.Figure 6
Process of rewiring a random tree.
(a)(b)The core steps of the improvedRRT∗ algorithm lie in the above two processes, by which the redundant paths can be reduced and optimal paths can be found.


## 4.1. RRT∗ Algorithm

Traditional path planning algorithms such asA∗, shortest path, and genetic algorithms, which require modeling the obstacles in a defined space, are not suitable for solving the path planning of a multidegree-of-freedom aircraft in a complex environment [31]. The fast extended random tree algorithm, that is, the RRT∗ algorithm, can solve path planning problems in high-dimensional spaces and under complex constraints by performing collision detection on sampled points in the state space. The intention is to generate random points and then search for a step forward from the starting point to the target point. In this manner, it can effectively avoid obstacles, prevent the path from falling into the local minima, and converge faster [32].According to Figure3, the steps of the RRT∗ algorithm are as follows.Figure 3
Expansion diagram of the RRT∗ algorithm.Initialize the random tree, adding nodes starting with an arbitrary random tree with onlyQinit.
(1)
Execute thesample function to obtain a random point Qnearest in the mapIterate through all the nodes in the tree and find the pointQnearest with the least cost to reach Qrand.
(2)
Execute theextend function to obtain the extension point Qnew of the specified length in the direction of Qnearest to Qrand and perform collision detection on Qnew. If false, Qnearest is designated as the parent node of Qnew, and the line connecting the two points is drawn(3)
Determine whetherQnew has reached the specified target range, if it has, end the loop; otherwise, continue to execute the loop until the target point Qrand is found(4)
To increase the practicality of the algorithm, a running time limit and a limit on the number of nodes searched are added to the above process. In other words, the algorithm is judged to have failed after searching for a finite number of nodes within the specified time and not finding the target point


## 4.2. Improved F-RRT∗ Algorithm

To solve the problems of strong randomness, poor orientation, slow convergence speed, and poor smoothness of path twisting and tracing of thetraditionalrapidlyexploringrandomtree∗algorithm, an improved RRT∗ algorithm based on a target-biasing strategy combined with angle restriction is proposed—the F−RRT∗ algorithm.First, a target bias threshold is set based on the traditional RRT∗ algorithm, and the randomness of tree growth is reduced by expanding the internal nodes of the tree through the target bias strategy [33]. Simultaneously, a local expansion mechanism is introduced to avoid the local optimum problem caused by changing the sampling structure. Then, a bidirectional tree growth strategy is used to improve the convergence speed of the algorithm. Finally, redundant nodes are removed, and cubic spline interpolation is used to increase the smoothness of the path.


### 4.2.1. Extended Step

In the RRT∗ algorithm, Qinit is extended to Qnearest by a step to obtain the sampling point Qnew, expressed as
(17)Qnew=Qinit+φQnearest−QinitLQinit,Qnearest,where φ is the extended step size, and LQinit,Qnearest is the Euclidean distance between the points Qinit and Qnearest.


### 4.2.2. Target Bias Strategy

The traditional RRT∗ algorithm has a large range of sampling points and is prone to randomness; therefore, the algorithm search time is longer [34]. By adopting the target bias strategy, the target points are used as sampling points for random tree expansion with a certain probability P, which reduces the number of sampling points and accelerates the efficiency of the random tree expansion [35].


### 4.2.3. Metric Function

The objective of the metric function is to find the point closest to the selected pointQinit in the randomly expanded tree, namely, Qnearest. In the traditional RRT∗ algorithm, the metric function of the Euclidean distance is used [36]. The improved F-RRT∗ algorithm introduces an angle constraint while considering two-point spacing, which makes the entire path smoother.The distance and angle are unified such thatLd represents the distance function, and Mθ represents the declination function after planning. The node that minimizes the function Mxi,xj is near the node formed in the search extension tree, which can be used as Qnew, as shown in Figure 4.
(18)Ld=dmax−ddmax,Mθ=θmax−θθmax,Li,j=xj−xi,θi,j=θj−θi,Mxi,xj=LdLi,j+Mθθi,j,where the coordinates of the point Qinit are xi,yi, and those of Qnearest are xj,yj; dmax and θmax are the maximum distance and maximum angle of the search, respectively.Figure 4
Extended diagram of the F-RRT∗ algorithm.


### 4.2.4. Path Pruning and Smoothing

The sampling points of the traditional RRT algorithm are highly random, and the generated paths often have numerous redundant nodes [37]. The efficiency in complex environments with many obstacles is low, making path pruning an essential step. The newly generated path nodes are filtered based on the traditional paths to remove redundant nodes, and the pruned paths are then simplified using a greedy algorithm [38].The specific pruning steps are as follows:Begin from the starting pointQinit and find all its Qnearest nodes. If the connection between Qinit and one of its Qnearest nodes does not contain obstacles, the Qnearest node is taken as the new node Qnew, and all other Qnearest nodes are deleted. Conversely, if the connection between Qinit and one of its Qnearest nodes encounters obstacles, then the above process is repeated starting from the parent node of the point Qinit.
(1)
Find the connection betweenQinit and the next Qnearest node until the Qrand node is added to the new path. Finally, the identified nodes are connected by a smooth curve [39]TheK-times B-spline curve equation is as follows:
(19)Pu=∑i=onBi,kuPii=1,2⋯n,where Pi is the node of the curve control, Bi,ku is the K-times B-spline basis function, and the node vector u=u0,u1,⋯,un+k+1 determines the K-times segmented curve [40].The basis functions of the tripleB-spline curve are as follows:
(20)b0=16−u3+3u2−3u+1,b1=163u3−6u2+4,b2=16−3u3+6u2+3u+1,b3=16u3.


## 4.2.1. Extended Step

In the RRT∗ algorithm, Qinit is extended to Qnearest by a step to obtain the sampling point Qnew, expressed as
(17)Qnew=Qinit+φQnearest−QinitLQinit,Qnearest,where φ is the extended step size, and LQinit,Qnearest is the Euclidean distance between the points Qinit and Qnearest.


## 4.2.2. Target Bias Strategy

The traditional RRT∗ algorithm has a large range of sampling points and is prone to randomness; therefore, the algorithm search time is longer [34]. By adopting the target bias strategy, the target points are used as sampling points for random tree expansion with a certain probability P, which reduces the number of sampling points and accelerates the efficiency of the random tree expansion [35].


## 4.2.3. Metric Function

The objective of the metric function is to find the point closest to the selected pointQinit in the randomly expanded tree, namely, Qnearest. In the traditional RRT∗ algorithm, the metric function of the Euclidean distance is used [36]. The improved F-RRT∗ algorithm introduces an angle constraint while considering two-point spacing, which makes the entire path smoother.The distance and angle are unified such thatLd represents the distance function, and Mθ represents the declination function after planning. The node that minimizes the function Mxi,xj is near the node formed in the search extension tree, which can be used as Qnew, as shown in Figure 4.
(18)Ld=dmax−ddmax,Mθ=θmax−θθmax,Li,j=xj−xi,θi,j=θj−θi,Mxi,xj=LdLi,j+Mθθi,j,where the coordinates of the point Qinit are xi,yi, and those of Qnearest are xj,yj; dmax and θmax are the maximum distance and maximum angle of the search, respectively.Figure 4
Extended diagram of the F-RRT∗ algorithm.


## 4.2.4. Path Pruning and Smoothing

The sampling points of the traditional RRT algorithm are highly random, and the generated paths often have numerous redundant nodes [37]. The efficiency in complex environments with many obstacles is low, making path pruning an essential step. The newly generated path nodes are filtered based on the traditional paths to remove redundant nodes, and the pruned paths are then simplified using a greedy algorithm [38].The specific pruning steps are as follows:Begin from the starting pointQinit and find all its Qnearest nodes. If the connection between Qinit and one of its Qnearest nodes does not contain obstacles, the Qnearest node is taken as the new node Qnew, and all other Qnearest nodes are deleted. Conversely, if the connection between Qinit and one of its Qnearest nodes encounters obstacles, then the above process is repeated starting from the parent node of the point Qinit.
(1)
Find the connection betweenQinit and the next Qnearest node until the Qrand node is added to the new path. Finally, the identified nodes are connected by a smooth curve [39]TheK-times B-spline curve equation is as follows:
(19)Pu=∑i=onBi,kuPii=1,2⋯n,where Pi is the node of the curve control, Bi,ku is the K-times B-spline basis function, and the node vector u=u0,u1,⋯,un+k+1 determines the K-times segmented curve [40].The basis functions of the tripleB-spline curve are as follows:
(20)b0=16−u3+3u2−3u+1,b1=163u3−6u2+4,b2=16−3u3+6u2+3u+1,b3=16u3.


## 4.3. Steps to Improve the RRT∗ Algorithm

The main difference between the improved RRT∗ algorithm and RRT∗ algorithm is that the former has two additional recomputation processes for the new node Qnew, which are as follows:
(1)
The process of reselecting the parent node forQnew contains an additional reselection process(2)
An additional rewiring random tree processBeginning from the starting pointQinit, we search for the “nearest neighbors” in the selected range; that is, we search for the suitable Qnearest node among all the Qnearest nodes from Qinit as an alternative to the parent node Qnew. The distance from the starting point Qinit to the nearest Qnearest node and the distance from the nearest Qnearest node to the parent node Qnew are calculated [41]. The specific steps are illustrated in Figure 5. Figure 5(a) shows the expansion process of the random tree at a certain moment, where node 0 represents the starting point Qinit, node 9 is the newly generated node Qnearest, node 6 is the parent node Qnew of node 9, and the generated path sequence number is 0-4-6-9 with a distance cost of 10+5+1=16. Figure 5(b) shows the process of identifying the parent node once again based on Figure 5(a). The newly generated Qnearest node, that is, node 9, is considered as the center, and its “near neighbors” are detected within the selected range [42], that is, nodes 5, 6, and 8; they are Qnearest composed of paths 0-1-5-9, 0-4-6-9, and 0-1-5-8-9, and the distances are 3+5+3=11, 10+5+1=16, and 3+5+1+3=12, respectively. It can be seen that node 5 is the closest to node 9; therefore, the parent node of node 9 is changed from the original node 6 to node 5, and the results are shown in Figure 5(b).Figure 5
Process of reselecting the parent node.
(a)(b)After reselecting the parent node, the rewiring process is performed on the random tree to further reduce the distance cost of the random extended tree [43]. The specific steps are illustrated in Figure 6. Figure 6(a) shows the rewiring process of the random tree at a certain moment, where node 0 is the starting point Qinit and node 9 is the Qnew node; that is, node 9 is the center, and its “nearest neighbors” are found in the selected range, which are nodes 5, 6, and 8. The distances from the starting point to each “nearest neighbor” node are 8, 15, and 9, respectively. Figure 6(b) shows the rewiring process based on Figure 6(a). Suppose that the parent node of node 5 is not node 1 but node 9, and the path is 0-4-9-5 with distance 17, which is greater than the original path distance 8, then the parent node of node 5 is not changed. Similarly, find the new parent node of node 6. Assuming that it is node 9 with path 0-1-5-9-6 and distance 12, which is less than the original path distance of 15, the parent node of node 6 is changed from node 4 to node 9. Finally, we find the new parent node of node 8. Assuming it is node 9, the path is 0-1-5-9-8, and the distance is 14, which is greater than the original path distance 9; therefore, the parent node of node 8 is not changed [44]. In summary, only the parent node of node 6 must be changed, and the newly generated paths are shown in Figure 6(b). From the above, it can be seen that the significance of rewiring is to reduce the distance cost of some newly generated node paths, thus reducing the distance of the entire optimized path.Figure 6
Process of rewiring a random tree.
(a)(b)The core steps of the improvedRRT∗ algorithm lie in the above two processes, by which the redundant paths can be reduced and optimal paths can be found.


## 5. Simulation Verification

Southern China is a thunderstorm-prone region. Hence, this study considered the typical route of Guiyang–Guangzhou in China as an example and studied the delineation of the flight-restricted area and optimization of the route affected by thunderstorms and traffic disruption. The Guiyang–Guangzhou route comprises the route points Guiyang–Guilin–Gaoyao–Pingzhou–Guangzhou, as shown in Figure7.Figure 7
High-altitude aerial map of Guiyang–Guangzhou route in China.


### 5.1. Flight-Restricted Area at the Start of the Change of Flight

Through meteorological satellite and Doppler meteorological radar observation images, the Chinese aviation meteorological department observed a band of thunderstorm clouds in the Guiyang–Guilin flight segment and some thunderstorm clouds in the Guilin–Gaoyao flight segment. According to the method introduced in Section1, the initial flight-restricted area is set up; the area has a radarecho≥41dBZ at a certain time t, as shown in Figure 8 [45]. To facilitate the analysis, the flight restricted area obtained following the processing in Figure 8 is analyzed separately in Figure 9. To simplify the graphical representation, Guiyang is denoted as A, Guilin as B, Gaoyao as C, Pingzhou as D, and Guangzhou as E in Figure 9. The boundary points are s1–s10, with the location of Guiyang Airport in China as the origin, magnetic north as the positive direction of the y-axis, and magnetic west as the positive direction of the x-axis [46]. The coordinates of Guiyang, Guilin, and Gaoyao in this map are (0, 0), (152.50, −53.26), and (394.30, −256.72), respectively. The coordinates of S1–S10 are (213.64, −24.78), (265.64, −153.32), (299.64, −97.82), (304.14, −219.32), (249.64, −182.82), (202.14, −187.20), (174.14, −5.32), (174.14, −32.82), (97.26, −34.14), and (97.26, −6.33), respectively. The lumpy thunderstorm cloud AFFA was observed to be moving in the direction of 40° east of north; the rates of movement were v1=40km/h and v2=10km/h, and the safety margin was L1=25km; the expansion speed of the thunderstorm boundary point was Vs2=20km/h [47]. The direction of movement of the banded thunderstorm cloud BFFA was 20° south by west; the rates of movement were v3=30km/h and v4=10km/h, and the safety margin was L2=25km; the expansion speed of the thunderstorm boundary point was Vs1=20km/h, and the acceleration was a=5km/h2. The flight from Guiyang to Guangzhou was expected to start at a certain time t [48], t0=0.5h after arriving in Guilin. The cruising speed was vs=800km/h.Figure 8
The initial flight-restricted areas.Figure 9
Initial boundary of hazardous weather on the Guiyang–Guangzhou route in China.Based on the initial flight-restricted area, the dynamic flight-restricted areas AFFA2 and BFFA2 are determined according to the prediction method introduced in Sections2 and 2.1, as shown in Figure 10. Owing to the uncertainty of the thunderstorm movement, the boundary of the initial flight restricted area is shifted outward by 25 km [49]. Then, it is expanded outward with a speed v=20km/h and an acceleration a=3km/h using the uniform speed prediction method to obtain the finalized dynamic flight-restricted areas AFFA2 and BFFA2. First, the movement of the flight-restricted area was considered according to the boundary point Si of the affected area by using the Graham algorithm. The initial dynamic flight-restricted areas AFFA1 and BFFA1 were determined after an elapsed time t=0.5h, as shown in Figure 10. Because of the uncertainty of the thunderstorm movement, the flight-restricted area boundary was shifted outward by 25 km [49]. Considering the variation in the thunderstorm boundary, the finalized dynamic flight-restricted areas AFFA2 and BFFA2 were obtained by expanding outward with velocity v=20km/h and acceleration a=3km/h2 [50]. The coordinates of the vertices are (210.44, −4.70), (255.64, −53.32), (319.64, −96.82), (314.14, −239.32), (233.64, −221.82), (197.14, −197.20), (144.14, 4.33), (144.14, −42.82), (76.26, −42.82), and (76.26, 4.33) assuming that the actual flight-restricted areas are AFFA2 and BFFA2, as shown in Figure 11, and path planning is performed based on them.Figure 10
Dynamic flight-restricted area on the Guiyang–Guangzhou route in China.Figure 11
The finalized dynamic flight-restricted areas.The above figure depicts the delineation of the flight-restricted area. The next step is to set the model parameters and solve the model. First, we set up the simulation environment with the location of Guiyang Airport in China as the origin, magnetic south as the positivey-axis direction, magnetic west as the positive x-axis direction, starting point as [0,0], ending point as [860,480], and extended step as 20 m; the standard deviation σ=5.6,kp=420000. The standard threshold of the path length was 586 m. The search was terminated when the number of sampled nodes was greater than 18000 or the path length was less than the standard threshold.The results obtained using the improved F-RRT∗ algorithm based on the above parameters are shown in Figure 12. From the following graphs, it can be seen that this algorithm enables the flight to bypass the flight-restricted area affected by dangerous weather more smoothly and finds the least costly path more quickly when compared with the other algorithms [51].Figure 12
Path planning effect.The paths obtained by the following four algorithms were compared and analyzed in terms of the path planning time, number of turning points, and path distance. The analysis results are shown in Figure13. In terms of the path planning time, the F-RRT∗ algorithm uses the shortest time to identify a path. The effect is significant as the number of sampling points increases. In terms of the number of turning points, it can identify more valid nodes, whereas in terms of the path distance, it also determines the shortest path. The specific parameters are shown in Table 1. The simulation results show that the F−RRT∗ algorithm has a shorter path planning time, fewer turning points, and shorter and smoother paths than the other RRT∗ algorithms. It can be better used for flight path planning.Figure 13
Comparison and analysis of the capability of four algorithms.Table 1
Comparative performance analysis of various RRT∗ algorithms.
AlgorithmAverage planning time (s)Average number of sampled nodes in spaceAverage number of effective nodesAverage path length (km)Traditional RRT∗ algorithm45.6613559641847.49Bi-RRT∗ algorithm43.5611688971693.65PAVS-RRT∗ algorithm50.6912679121280.77F-RRT∗ algorithm39.6810658431183.84


## 5.1. Flight-Restricted Area at the Start of the Change of Flight

Through meteorological satellite and Doppler meteorological radar observation images, the Chinese aviation meteorological department observed a band of thunderstorm clouds in the Guiyang–Guilin flight segment and some thunderstorm clouds in the Guilin–Gaoyao flight segment. According to the method introduced in Section1, the initial flight-restricted area is set up; the area has a radarecho≥41dBZ at a certain time t, as shown in Figure 8 [45]. To facilitate the analysis, the flight restricted area obtained following the processing in Figure 8 is analyzed separately in Figure 9. To simplify the graphical representation, Guiyang is denoted as A, Guilin as B, Gaoyao as C, Pingzhou as D, and Guangzhou as E in Figure 9. The boundary points are s1–s10, with the location of Guiyang Airport in China as the origin, magnetic north as the positive direction of the y-axis, and magnetic west as the positive direction of the x-axis [46]. The coordinates of Guiyang, Guilin, and Gaoyao in this map are (0, 0), (152.50, −53.26), and (394.30, −256.72), respectively. The coordinates of S1–S10 are (213.64, −24.78), (265.64, −153.32), (299.64, −97.82), (304.14, −219.32), (249.64, −182.82), (202.14, −187.20), (174.14, −5.32), (174.14, −32.82), (97.26, −34.14), and (97.26, −6.33), respectively. The lumpy thunderstorm cloud AFFA was observed to be moving in the direction of 40° east of north; the rates of movement were v1=40km/h and v2=10km/h, and the safety margin was L1=25km; the expansion speed of the thunderstorm boundary point was Vs2=20km/h [47]. The direction of movement of the banded thunderstorm cloud BFFA was 20° south by west; the rates of movement were v3=30km/h and v4=10km/h, and the safety margin was L2=25km; the expansion speed of the thunderstorm boundary point was Vs1=20km/h, and the acceleration was a=5km/h2. The flight from Guiyang to Guangzhou was expected to start at a certain time t [48], t0=0.5h after arriving in Guilin. The cruising speed was vs=800km/h.Figure 8
The initial flight-restricted areas.Figure 9
Initial boundary of hazardous weather on the Guiyang–Guangzhou route in China.Based on the initial flight-restricted area, the dynamic flight-restricted areas AFFA2 and BFFA2 are determined according to the prediction method introduced in Sections2 and 2.1, as shown in Figure 10. Owing to the uncertainty of the thunderstorm movement, the boundary of the initial flight restricted area is shifted outward by 25 km [49]. Then, it is expanded outward with a speed v=20km/h and an acceleration a=3km/h using the uniform speed prediction method to obtain the finalized dynamic flight-restricted areas AFFA2 and BFFA2. First, the movement of the flight-restricted area was considered according to the boundary point Si of the affected area by using the Graham algorithm. The initial dynamic flight-restricted areas AFFA1 and BFFA1 were determined after an elapsed time t=0.5h, as shown in Figure 10. Because of the uncertainty of the thunderstorm movement, the flight-restricted area boundary was shifted outward by 25 km [49]. Considering the variation in the thunderstorm boundary, the finalized dynamic flight-restricted areas AFFA2 and BFFA2 were obtained by expanding outward with velocity v=20km/h and acceleration a=3km/h2 [50]. The coordinates of the vertices are (210.44, −4.70), (255.64, −53.32), (319.64, −96.82), (314.14, −239.32), (233.64, −221.82), (197.14, −197.20), (144.14, 4.33), (144.14, −42.82), (76.26, −42.82), and (76.26, 4.33) assuming that the actual flight-restricted areas are AFFA2 and BFFA2, as shown in Figure 11, and path planning is performed based on them.Figure 10
Dynamic flight-restricted area on the Guiyang–Guangzhou route in China.Figure 11
The finalized dynamic flight-restricted areas.The above figure depicts the delineation of the flight-restricted area. The next step is to set the model parameters and solve the model. First, we set up the simulation environment with the location of Guiyang Airport in China as the origin, magnetic south as the positivey-axis direction, magnetic west as the positive x-axis direction, starting point as [0,0], ending point as [860,480], and extended step as 20 m; the standard deviation σ=5.6,kp=420000. The standard threshold of the path length was 586 m. The search was terminated when the number of sampled nodes was greater than 18000 or the path length was less than the standard threshold.The results obtained using the improved F-RRT∗ algorithm based on the above parameters are shown in Figure 12. From the following graphs, it can be seen that this algorithm enables the flight to bypass the flight-restricted area affected by dangerous weather more smoothly and finds the least costly path more quickly when compared with the other algorithms [51].Figure 12
Path planning effect.The paths obtained by the following four algorithms were compared and analyzed in terms of the path planning time, number of turning points, and path distance. The analysis results are shown in Figure13. In terms of the path planning time, the F-RRT∗ algorithm uses the shortest time to identify a path. The effect is significant as the number of sampling points increases. In terms of the number of turning points, it can identify more valid nodes, whereas in terms of the path distance, it also determines the shortest path. The specific parameters are shown in Table 1. The simulation results show that the F−RRT∗ algorithm has a shorter path planning time, fewer turning points, and shorter and smoother paths than the other RRT∗ algorithms. It can be better used for flight path planning.Figure 13
Comparison and analysis of the capability of four algorithms.Table 1
Comparative performance analysis of various RRT∗ algorithms.
AlgorithmAverage planning time (s)Average number of sampled nodes in spaceAverage number of effective nodesAverage path length (km)Traditional RRT∗ algorithm45.6613559641847.49Bi-RRT∗ algorithm43.5611688971693.65PAVS-RRT∗ algorithm50.6912679121280.77F-RRT∗ algorithm39.6810658431183.84


## 6. Conclusion

This paper presented an in-depth study on the redirecting strategy of flight paths under hazardous weather conditions. An improved F-RRT∗ algorithm based on the target-biasing strategy combined with angle restriction was proposed. The traditional RRT∗ algorithm has been conventionally used for robot path planning; it has never been used for flight-path planning. In this study, we improved the traditional RRT∗ algorithm based on the characteristics of flight and derived an algorithm applicable to flight-path optimization, that is, the F-RRT∗ algorithm.The main aspects of this study and its findings are as follows.(1)
The radar data for the Guiyang–Guangzhou route in China were obtained, processed and filtered to determine the initial hazardous weather area on the route based on the radar echo values. Then, theGraham algorithm was proposed to predict the flight-restricted area according to the boundary change of hazardous weather to obtain the final flight-restricted area on the route. The flight-restricted area was divided into three categories based on its shape and distribution characteristics: block-shaped, strip-shaped, and scattered. The concepts related to changes in flight were briefly explained(2)
Based on the above analysis, the model was constructed as a grid environment model with the flight route distance, route turn size, and number of turning points as the objectives, and the model was explained(3)
An improved F-RRT∗ algorithm based on a target-biasing strategy combined with angle restriction was proposed based on the flight characteristics. It solved the problems of traditional path planning algorithms, such as strong randomness, poor guidance, slow convergence speed, poor path curvature, and tracing smoothness(4)
The model obtained under this algorithm was simulated and analyzed for the Guiyang–Guangzhou route in China as an example. The obtained rerouting paths were compared in terms of the path planning time, path distance, and number of turning points. The simulation results showed that the F-RRT∗ algorithm had shorter planning time, fewer turning points, smaller turning angle, and smoother path than the other RRT∗ algorithms. This method is effective in reducing flight delays while ensuring safety interval standardsAlthough the achievements of this study are significant, several aspects can be further improved in future research, for example, the establishment of a more detailed dynamic environment model, improvement of the metric for bad weather, and consideration of three-dimensional modeling of the flight paths. It is also possible to conduct in-depth research on multiflight and multiroute flight plans and expand from diversion path planning for one route to that for route networks.



---

*Source: 1166968-2022-08-22.xml*