Abstract
Under convective weather conditions, aircraft rerouting in terminal airspace is essential to ensure flight safety and reduce air traffic delays. Traditional air traffic rerouting approaches do not combine convective weather information using radar image recognition with controller workload, additional fuel consumption, delay loss, time cost, and route length at terminal airspace in tactical ATFM phase for exact rerouting decision. In accordance with the safety and economic principles of the route network in terminal airspace and in consideration of the changes in the speed and height of aircraft in terminal airspace, a multiobjective rerouting planning model is established in this study for terminal airspace under convective weather conditions in tactical air traffic flow management phase. Then traffic simulation is conducted to analyze the capacity and delays of the rerouting in the approach path in Shanghai terminal area. Experimental results indicate that this model can increase airspace capacity and operational efficiency of air traffic compared with the traditional air traffic rerouting approaches.
1. Introduction
With the rapid development of air transportation, the problem of air traffic congestion is exacerbated. The contradiction between limited airspace resources and traffic demand becomes increasingly prominent. As a result, flight delays, economic losses, and the risk of security incidents increase in an air transportation system. Convective weather is an important factor detrimental to the safety of aviation operation. Data show that weather factors primarily cause flight accidents, accounting for approximately 50% of all flight accidents [1]. Furthermore, weather has become the main causal factor in flight delays [2, 3]. According to an in-depth research of the Federal Aviation Administration (FAA) of the United States, approximately 65%–75% of flight delays longer than 15 minutes are directly caused by weather conditions [2]. Estimating airport capacity to formulate appropriate traffic management strategies is remarkably difficult because of the uncertainty of weather forecast for terminal airspace. In the 12 airports included in the National Airspace System (NAS) Plan released by the FAA in 2008 and 2009, the end time of the weather-related ground delay program was averagely 95 minutes earlier than that of the initially planned ground delay program, which was nearly 2 hours earlier than the end time of the revised ground delay program [4, 5]. For the Next Generation Air Transportation System, the FAA and National Aeronautics and Space Administration study the reduction of the negative impact of weather conditions on NAS. They also intend to explore improving the Collaborative Air Traffic Management by increasing the degree of automation and improving decision support tools [6].
The FAA determined the weather-affected area of the air traffic management unit according to the weather radar echo and the 30-minute earlier short-term forecast to investigate the impact of inclement weather conditions on route network analysis in tactical air traffic flow management (ATFM) phase. The following process was adopted: a grid pattern was used to describe the entire environment, and the area influenced by severe weather condition was divided into 0–6 grades [7–9], which were represented by different colors. A study of the Lincoln Laboratory at the Massachusetts Institute of Technology shows that the third- or higher-grade weather condition can cause aircraft bumpiness and consequently reduce flight performance [10]. When the intensity of the echo is higher than 41 dBZ, the flight is unauthorized to proceed. The application of this criterion can be used to identify the restricted flight area clearly according to the route; however, the impact of the thunderstorm area on route coverage should be considered because of the structural constraints in the terminal area to improve and perfect the method applicable to it. If the flight phase is affected by the weather condition, the maximum flow/minimum cut theorem can be used to analyze the capacity of the affected sector [7, 11]. In the NAS of the United States, the impact of convective weather condition on the route of jet aircraft should be assessed. For this impact assessment, a convective weather avoidance model (CWAM), weather avoidance field, maximum flow/minimum cut sets, and route-blocking rate are used [8, 12]. The storm area is composed of cumulonimbus clouds and is a closed curve with inward increasing precipitation grades. A geometric method can be used to draw the shape of thunderstorm clouds according to the characteristics of the cumulonimbus clouds. This method can also effectively provide reference for rerouting according to the vertex of the convex curve. The maximum flow/minimum cut theorem can be used in a route to analyze the influence of a thunderstorm area in a terminal area on the degrees of route coverage. The coverage area of a thunderstorm and the intensity of precipitation in a thunderstorm area can be combined in accordance with statistical theory to analyze the traffic condition in a terminal area when considering the intensity of precipitation and topological structure of the terminal area.
Many air traffic management experts have studied route network in terminal airspace and established rerouting path planning models for aircraft to avoid the effect of complex weather conditions on flights. The existing studies on this topic are conducted from five different aspects. First, a stochastic dynamic planning model and a probabilistic network model are established, and the decision of rerouting in terminal airspace is made based on the probability and random distribution of weather and traffic events. For example, McCrea et al. (2008) [13] proposed a model for large-scale airspace planning and collaborative decision-making to generate the rerouting path under the constraints of flight safety, operation quality of airline, and controller workload in a sector under severe weather conditions. Second, a geometric model is established for rerouting in airspace to avoid severe weather conditions, and a dynamic rerouting path is established for aircraft in dynamic weather conditions. For example, Rubnich and Delaura (2010) [14] proposed a CWAM in en route airspace. Matthews and DeLaura (2010) [15] quantified the performance in a weather-affected area to predict the influence of convective weather on the operation on en route airspace. A geometric compensation model was proposed by Yoon et al. (2011) [16]. Third, an integer programming model for an aircraft in terminal airspace is established for the planning of rerouting path of a dynamic route network in terminal airspace. For example, a dynamic airspace was proposed by Michalek and Balakrishnan (2010) [17] for a stochastic model of route availability. Zou (2010) [18] used a flexible tree based on the approach and departure routes in transition airspace to establish a flow control model for aircraft along the approach path under hazardous weather, special use airspace, and operational constraints on merge points. Bertsimas et al. (2011) [19] proposed an integer programming method to solve the problems in large-scale air traffic management. In this method, all the flight phases of aircraft are covered to obtain the prioritization scheme of traffic management. Agustin et al. (2012) [20] analyzed the route network under severe weather conditions by constructing the scene tree of severe weather conditions and establishing a stochastic mixed 0-1 optimization model under uncertainty in the airport arrival and departure capacity, the sector capacity, and the flight demand. Dell’Olmo and Lulli (2003) [21] proposed a two-level structure for air traffic management by establishing a mathematical model; the structure, which considered the aircraft traffic flow and the capacity limit of the operational route, was based on the static route network of the aircraft traffic flow and the distribution of traffic flow. Bertsimas and Patterson (2000) [22] proposed the rerouting of aircraft based on a dynamic network flow approach in consideration of the dynamically changing weather conditions. Janić (2015) [23] proposed an evaluation method for the resilience, friability, and costs of air transportation network. Fourth, some scholars have attempted to improve the algorithm of the traditional rerouting planning model to establish dynamic rerouting path planning models suitable for complex weather conditions and consequently enhance the accuracy of optimization results. One of these algorithms is the Lagrange method proposed by Patterson (1997) [24]. Samà et al. (2012, 2014, and 2016) [25–27] proposed mixed-integer linear programming for traffic rerouting in terminal airspace, implemented centralized and rolling horizon traffic control paradigms, and presented a number of algorithmic improvements implemented in the AGLIBRARY solver in order to improve the possibility of finding good quality solutions quickly. Gangi et al. (2016) [28] developed an optimization method for signal settings consistent with traffic flow dynamics and used a heuristic algorithm and multicriteria optimization. Fifth, some scholars have analyzed a large amount of simulation and traffic survey data in a terminal area using big data statistical analysis and machine learning techniques to explore the rerouting path planning method for air traffic in a terminal area under complex weather conditions [29–32].
The above-mentioned study is conducted from five different aspects based on methodology. Meanwhile, air traffic rerouting problem can be divided into strategic, pretactical, and tactical phases in accordance with ATFM. In the tactical ATFM phase, the effective rerouting management of traffic flows during convective weather events in terminal airspace requires decision support tools that can translate weather information into anticipated air traffic operational impact. Some related studies focus on a geometric model [14, 17, 33], big data statistical analysis and machine learning techniques [9, 34], and a probabilistic network model [10, 13, 35] in order to provide pilot with rerouting suggestion from convective weather forecast systems.
Although the studies presented have solved the air traffic flow management rerouting problem in the en route flight phase under dynamic weather conditions, they mainly focused on the economic costs of aircraft delays to airlines as a target; however, they did not consider the controller workload. In the literature, most of the covered rerouting problems focus on the en route flight phase. The restricted area is only extrapolated as a geometric shape. Moreover, these literatures do not combine convective weather information using radar image recognition with controller workload, additional fuel consumption, delay loss, time cost, and route length at terminal airspace in tactical ATFM phase for exact rerouting decision. Thus, the rerouting solutions offered in the literature cannot completely meet the operational requirements because of the dynamic changes in the topological network structure of the terminal area.
To overcome the drawbacks of previous methods, a new method is proposed in this study. In this method, the factors influencing the normal operation of aircraft are first identified and classified based on the precipitation grade; the digital images of radar echoes are processed; and the probability model for traffic rerouting is established. Based on stochastic dynamic programming, traffic volume of each flight path is considered as load weights, each of which is assigned to each arc segment of the approach path. A rerouting model is established for a dynamic multiobjective route network in a terminal area under short-term convective weather conditions. Furthermore, the optimization problem of the dynamic network in the terminal area is solved. Thus, the present study presents the foundation for further research on the management of air traffic flow.
The remainder of this paper is organized as follows. In Section 2, the probability model for traffic rerouting is established by obtaining the contour lines of the restricted thunderstorm areas with different precipitation levels, the dynamic airspace network in the terminal area, and radar echoes. In Section 3, a multiobjective programming model is established for additional fuel consumption, delay loss, time cost, and route length; this model is based on the topological structure of the dynamic network in the terminal airspace. Section 4 validates the proposed model by applying it to the approach path in Shanghai terminal areas. Finally, Section 5 presents a summary of the current study and research prospects.
2. Establishment of the Probability Model for the Traffic Rerouting in Thunderstorm Areas
2.1. Detection of Rerouting according to the Color Index
A thunderstorm area has different precipitation levels; thus, the probabilities of rerouting aircraft in the terminal sector vary across regions. The region with high precipitation has a relatively high traffic rerouting probability. Therefore, in this study, the precipitation levels in the thunderstorm area are categorized according to the regional rainfall intensity represented by different color values in the cloud chart. Subsequently, the probability of traffic rerouting in different regions can be analyzed.
In this study, represents the number of aircraft that pass the region in the terminal sector under -level precipitation; represents the number of aircraft that pass this region and need to be rerouted. We can calculate based on geographic location of aircraft under -level precipitation using contour extraction of a region using the binary image. For the region in the terminal sector under -level precipitation, the rerouting probability is calculated by Precipitation is categorized into four levels: , , , and , as shown in Tables 1–4.
O-level precipitation is minimal; thus, in a thunderstorm area with this precipitation level, an aircraft does not need to consider the weather condition, and it can directly fly through the restricted flight zone, in which the probability of rerouting is close to 0.
2.2. Identification of Thunderstorms by Digital-Image Processing
2.2.1. Acquisition of the Binary Image of an Area according to the Different Colors of Precipitation Levels
A binary image can be acquired by obtaining the color attributes of , , and . First, the RGB values of all the colors in the picture are extracted and then saved as the values of variables R, G, and B. Second, the relationship among the three colors (i.e., R, G, and B), the relational functions of the color values, and the threshold relationship are determined. Subsequently, the desired color is extracted. Third, the pixels that cannot meet the conditions are colored black. After the color extraction is finished, the fixed threshold method is used to acquire the binary image, as shown in Figure 1.
2.2.2. Contour Extraction of a Region Using the Binary Image
If is regarded as the gray value of any point in the image, then the gradient of the point is expressed asThe Sobel edge detection operators are used to detect the pixel matrix of the edge image. The Sobel operators are the matrices—A and B. represents the gradient along the -direction in the image, whereas represents the gradient along the -direction in the image. Subsequently, formulas (8), (9), and (10) are substituted into formulas (6) and (7) to determine the gradient size of the pixel point. The maximum gradient is detected, the pixel point is given the white property, and the contour line of the specified color region is determined. The contour maps of the regions with the , , and levels in the thunderstorm area are obtained using MATLAB (please see appendix which is codes of green contour, and the codes of other color contours are similar).
2.3. Influence of Thunderstorm Coverage on Rerouting Probability
In the terminal area of the route network, the thunderstorm coverage of the route affects the normal flight of the aircraft. In a previous study [11], the maximum flow/minimum cut theorem was applied to establish a thunderstorm area, adapting to the terminal area of the route network according to the following one-dimensional linear model of the thunderstorm area [11]:where is the width of the route, is the length of zones restricted by severe weather conditions, and is the flight width of Required Navigation Performance, that is, a single vehicle for a one-way traffic flow.
2.3.1. Condition 1: Unilateral Thunderstorm Area
In the maximum flow/minimum cut theorem, the minimum cut is the shortest distance from the protected zone, which is unaffected by the thunderstorm, to the thunderstorm area. The minimum cut is denoted by , where represents the single side width of the route, is the vertical distance between the most distant point in the thunderstorm area and the intersecting route, and represents the safety redundancy of the thunderstorm area, as illustrated in Figure 2.
The aircraft has equal probabilities of occupying the positions of the transversal of a specified height in the route. Thus, equidistribution can be used to describe the proportion of aircraft in the areas affected by thunderstorms. The probability density is expressed as where and are minimum and maximum. If and , then the probability density in formula (12) can be modified towhere represents the sphere of influence in the thunderstorm area when safety redundancy is considered, that is, the starting point of the minimum cut . By conducting field research and consulting pilots, the following conditions are assumed.
If , the probability of an aircraft flying in an unaffected area can be expressed as The proportion of the affected aircraft is expressed as follows:
If , which is in accordance with the actual operation of pilots, the entire flight path of the aircraft is affected by the thunderstorm area when the nominal route is covered with a track line. Under this condition, .
2.3.2. Condition 2: Bilateral Thunderstorm Area and the Minimum Cut Perpendicular to the Route
This scenario can be illustrated as a linear type, as shown in Figure 3.
Under this condition, suggests thatwhere represents the length of the minimum cut and is the shortest distance covering the opposite protected area. If both sides are covered, both sides of form the radial distance .
2.3.3. Condition 3: Bilateral Thunderstorm Area and the Minimum Cut Not Perpendicular to the Route
In the terminal area, the aircraft flies along the flight path; thus, deviating from the original nominal track line is impossible. The following relation can be used to judge whether the minimum distance of the two closed curves, in which the minimal cut stays, is satisfied:where and represent the vertical distance between the two sides of the distance from the thunderstorm area corresponding to the minimum cut and and denote the safety redundancy of the thunderstorm area.
If this condition is satisfied, then the route is covered, and ; otherwise, (see Figure 4).
2.3.4. Combination of the First Three Conditions
From what is mentioned above, we can see that the feasible aircraft paths are contained in the area defined by closed curves. The shape of the radar echo can be composed of a number of thunderstorm areas. Furthermore, the closed curves of each thunderstorm area have mutual independence. Thus, the effects of the curves on the segment of the route can be considered separately.
Therefore, based on the classification of precipitation levels and the given information, three precipitation grades (three different color values) exist. The influence on the coverage of the green area en route is considered. Each color belongs to the nested relation. Therefore, the width of the green area can be determined by assuming that the width of the green area is the width of the route. Furthermore, the influence of the yellow area on the route can be determined as follows: (a) when the yellow area covers the nominal track, then all the aircraft are affected; (b) when the yellow area does not cover the nominal track, then the influence ratio can be obtained according to the uniform distribution.
The influence ratio of the red area can be determined using the same procedure.
2.4. Determining the Flight Path in the Approach Segment Affected by Thunderstorms
2.4.1. Verification if the Route in the Terminal Area Passes through the Thunderstorm Area
The gray value of each pixel in the contour map can be obtained simultaneously based on the regional contours of the , , and precipitation levels. According to the implied conditions in the MATLAB program, the background of the image is black, whereas the contour line is white; that is, the gray value of the pixels in the background area is 0, whereas that of the pixels along the contour line is 255. Therefore, the coordinates of the pixels along the regional contour lines of the , , and precipitation levels are acquired by obtaining the coordinates of the pixels whose gray value is 255.
Any section of the route can be composed of several straight legs. Each linear route can be assumed as a monotonically increasing (decreasing) function within a certain range of the defined domain. When a point along the contour line satisfies the corresponding function condition, the straight leg can be judged as passing through the thunderstorm area.
Therefore, the contour map of the thunderstorm area with three different levels of precipitation can be obtained from the given information to determine whether the straight leg passes through the area with , , or precipitation level. The straight legs are merged, and the terminal area route that passes though the areas with , , and precipitation levels is determined.
2.4.2. Use of the Radar Echo Image to Analyze the Rerouting Probability in the Terminal Sector
The degree of coverage has an effect on the execution probability of a flight path, if the rerouting probability of segment , which is affected by the thunderstorm, is .
Therefore, the influence of the coverage of the precipitation area should be considered to judge the rerouting probability of the leg affected by regional precipitation. The rerouting probability of the leg may be denoted by , , or .
A leg can pass through three types of restricted thunderstorm areas: (a) the thunderstorm area is only the green area, (b) the thunderstorm area contains the yellow area, and (c) the thunderstorm area contains the red area.
As shown in the radar echo image in Figure 1, the area with a higher precipitation level is contained in the area with a lower precipitation level. Therefore, when the yellow area appears, it must be included in the green area. When the red area appears, it must be included in the yellow and green areas. The possible types of events are as in Figure 5.
The following notations are used in defining the probability of each event in Figure 5.
Definition 1. – represents the probability of Events 1–6, respectively.
Definition 2. represents the number of yellow curves that are surrounded by a green closed curve; represents the number of the yellow closed curves that are bypassed by the aircraft under the premise that they cross the green closed curve.
Definition 3. represents the number of red closed curves contained in the th yellow closed curve; represents the number of red closed curves, which are crossed when the aircraft is passing through the red closed curve contained in the th yellow closed curve.
According to these definitions, the probabilities of Events 1–6 can be expressed as follows: As demonstrated in Events 1–6, when combinations are formed, different flight paths may be executed when approaching the segment passing through restricted thunderstorm areas. With path as an example and as the number of events occurring in the case is denoted by , the probability of path isThus, when the nominal track passes through the restricted thunderstorm areas, the corresponding function and the restricted area will have a crossover point. The coordinates of each node can also be derived.
The abscissa of the crossover point that crosses the red area is denoted by and . Furthermore, the intersections of yellow and green are, respectively, denoted by and . All the coordinates of the points are ordered from small to large and constitute a sequence . Subsequently, the color of each value in the sequence is shown, and the following judgments can be made:
(a) The corresponding colors of are all green, and only Events 1 and 2 occur; .
(b) When the colors of are green and yellow, is green if the color of is yellow; or if the color of is green, then is yellow. Subsequently, the corresponding number of green crossover points can be extracted, and the total number of Events 3 and 4 is expressed as . The yellow crossover points among the adjacent green points can likewise be extracted, and the total number of Events 1 and 2 is expressed as .
(c) When the colors of contain red, is yellow if is red; or if is yellow, then is red. Subsequently, the coordinates of the first green crossover point in the direction can be determined subsequently and the coordinates of the first green crossover point can be found in the direction. The corresponding number of all green intersections can be extracted, and . Similarly, all the yellow nodes in the adjacent green intersections can be extracted, and the red intersections can be derived from the adjacent yellow nodes. The previous step can be executed to obtain the number of Events 1–4.
3. Rerouting Path Planning Model Based on Rerouting Probability in a Terminal Area
The assessment of controller workload is the basis for determining the capacity of terminal in the airspace sector [36]. The International Civil Aviation [37] has recommended a tool, DORATASK (Directorate of Operational and Analysis Task) method, to calculate the capacity of the airspace sector by using controller measures, particularly, the primary indicator, “time of directives issued by controller,” and the ideal outcome “less than 80% of worktime.” Subsequently, some scholars [38–40] have suggested that the statistical measure “controller workload” (in time units) be integrated to the airspace route structure. On the basis of past recommendations of ICAO and academic scholars, the present study analyzed the traffic delays of terminal airspace networks by means of determining the workload weights of each segment. In an attempt to update the original theory, the present study focused on rerouting of dynamic route topology of the terminal area.
The topological structure of the approach path in the entire terminal area can be regarded as a dynamic spatial structure network: , where node is a set of new route points in the rerouting path during operation and all the initial route points in the terminal area and arc is the flight phase among nodes in the arrival and departure routes. The overall best rerouting path that satisfies the overall objective of the terminal area is selected from the local optimal routes by setting the relevant constraints. The dynamic route network in the terminal area and the 0-1 hybrid model with low fuel consumption and time delay are used as the goals.
After the topological structure of the approach path is constructed, a multiobjective rerouting planning model is established for terminal airspace under convective weather conditions in accordance with the safety and economic principles of the route network. The optimization aims of this model are the minimal additional controller workload, minimal time cost loss of all aircraft, and the minimal cost of all length of the aircraft segment. The constraint condition covers additional fuel consumption, delay loss, time cost, and route length at terminal airspace. Finally, the four stages of work are established for the calculation of the planning model.
3.1. Parameters and Variables
The relevant parameters and variables of this model are defined in Notations.
3.2. Objective Function
In accordance with the safety and economic principles of the route network in terminal airspace, this multiobject optimization is proposed. Controller workload is one of the most important influencing factors of airspace capacity [11, 13]. So the first object function (formula (20)) is the minimal additional controller workload in rerouting model, and it complies with safety principle for route network. The second object function (formula (22)) is the minimal time cost loss of all aircraft; the third object function (formula (25)) is the minimal cost of all lengths of the aircraft segment. Formulas (22) and (25) comply with economic principle of the route network.
Formula (20) denotes the additional controller workload required for executing the rerouting path. In this formula, represents the aircraft in terminal airspace; represents the additional workload required when the controller commands the aircraft to fly across arc. If arc is in the planned leg, then is 0. According to the definition of probability density, when the number of aircraft involved in the flight mission is , then the number of aircraft able to execute the flight path is , and the number of aircraft unable to execute this flight path because of obstruction is . The elements used to define an arbitrary flight path contain attributes of the flight path and speed. That is, , , where is the flight path corresponding to the specified scenario, and is the average flight speed corresponding to the specified scenario.
Therefore, based on the weighted mean, the average flight speed for executing path is solved as follows: where is the linear distance between points and ; is the total execution time, and . Given , formula (21) can be expressed as . According to research on the controller workload, the controllers comply with the intention of the pilots and monitor the aircraft in an alternative flight segment. That is, the flight time of the aircraft is the additional monitoring load of controller . Therefore, when the flight distance is constant, the controller workload is . , , , and the scenarios of all the elements in the set are combined; the flight path will be executed to work with the additional load as ; the probability of the probability (safety parameter) is considered as the weight value; and the variables are summed. When all flight paths are considered at prediction time and when all aircraft pass through the terminal area of the route network , the workload consumed by the controller who commands the aircraft to pass the rerouting path composed of new route points can be obtained as follows: where time element is expressed as ; the objective function is the time cost loss of all aircraft ; time node is the time of arrival of each arc ; is the time cost of profit; and is the time cost of loss. Formulas (23) and (24), respectively, represent the starting and end points of the arc segment in arc ; the set of the arc segment is expressed as , where nodes and represent the beginning and end of the arc segment, respectively; is the set the initial arc; and . The specific elements in the set are .
When all the conditions in and are combined, the time differences and between the actual departure and initially planned times for each node are compared. The time losses and are considered as the weight values to sum up and derive the cost function of the topological network of the entire terminal area.In formula (25), time element is expressed as , where the objective function is the cost of all lengths of the aircraft segment; and are the nonnegative mixed variables that represent the -value of the actual leg length compared with the flight plan at node ; and are the profit and loss rates of the aircraft when flying in the sky, respectively; and and are the starting and end points of the approach procedure (or departure procedure) in the sector of the terminal area for the aircraft , respectively. is the flight time required for aircraft to execute the corresponding flight program in the terminal sector. When all the conditions in and are combined, the approach and departure routes in the topological network of the terminal area are selected by comparing the different values of the lengths in the approach and departure routes after a new waypoint is constructed with the original leg lengths and ; the loss of route lengths and are considered as the weight to sum up and derive the cost of the leg length.
3.3. Constraints
Formulas (26) and (27) represent the limit capacity of the airport in the terminal area and the limit capacity of the terminal airspace, respectively. Formula (26) indicates that when , the number of aircraft on the ground of airport does not exceed the theoretical value of the capacity of the ground. In both formulas, is the starting point of the approach procedure in the airport ; is the end point of the approaching procedure; is the starting point of the departure procedure in the airport ; is the end point of the departure procedure; represents the number of planes on the ground of the airport when ; and is the limit capacity of the ground in airport . Variable is used to sum , which represents the number of planes approaching airport at time . Variable is used to sum , which represents the number of planes for departure in airport at time .
In formula (27), is the total number of aircraft entering and leaving the terminal sector from the initial approach fix (IAF) at all airports in the terminal area at time ; is the total number of aircraft entering the terminal sector from the corridor; and is the total number of aircraft leaving the terminal sector from corridor. Formula (27) also indicates that the number of aircraft in the sector of the terminal area is not higher than the theoretical capacity of the airspace sector at time . Formula (28) represents that, depending on the 0-1 variable, aircraft can only perform one flight path at any moment .
The analytical method in formula (29) is the same as formula (30) to ensure the continuity among the nodes in each route. In formula (30), represents ; represents ; denotes the flight time of aircraft flying along the arc in the terminal sector; and and are the lower and upper values of the time interval , respectively, such that . Therefore, is expressed as the maximum value of the lower value of the time interval . When the flight path is performed with as the starting node, is the flight path that should be executed on the node to ensure that the two continuous legs perform the same flight program.
In formula (29), when the aircraft reaches point at time , and ; otherwise, . When the aircraft may reach the follow-up node at time , the value of is undefined.where , , and when aircraft enters from the outside of the sector in the terminal area. The aircraft must be handed off to the control tower at IAF point. When aircraft enters the sector from the IAF point, it is bound to leave the sector of the terminal area from the corridor, and control transfer is subsequently achieved. Formula (32) ensures the continuity of aircraft when flying across the leg . In this formula, is the decimal number that is sufficiently small, and , .
3.4. Algorithm Statement
The four stages of work should first be established prior to the calculation of the planning model, as shown in Figure 6. The first stage is the preparation and initialization of data, including traffic flow and air route information (Phase 1: Step 1). The second stage is the identification of constraint sets of relevant variables in their static route network conditions, including those related to flight plan, continuity and unilateral setup, and unique flight selection strategies; airspace and ground capacities in the terminal area are also determined (Phase 2: Steps 2–3). The third stage is the inclusion of weather information in thunderstorm conditions by using radar echo image recognition. This phase particularly aims to determine dynamic route networks and the probable distribution of route and workload weights of controllers (Phase 3: Steps 4–5). In the fourth stage, the LINGO software is used to calculate the optimal solution (Phase 4: Step 6).
The specific steps can be described and summarized as follows.
Step 1 (initialization). Data preparation and initialization mainly include historical data of traffic flow. Step 1 is conducted to determine the terminal region of every arrival and departure routes, the network nodes of each route, the ratio of aircraft flight sorties for the approach and departure routes, and the consumption ratio of mileage and fuel for each plane.
Step 2 (set constraints). Time parameters, the attributes of approach and departure legs, and properties and parameters of aircraft, and related variables are determined. The continuity and unilateral setup of legs (i.e., formulas (26) and (27)), as well as static route network and the ground capacity, are regarded as constraints, whereas the uniqueness of flight strategy (i.e., formula (25)) is considered as the limiting condition.
Step 3 (information loading). The flight plan is utilized to obtain the following information: voice and radar data, route network, workload mapping to the terminal airspace, and workload weights of the route network.
Step 4 (radar echo image recognition). From Sections 2.1 and 2.2, data on weather radar echo are established. Consequently, through image recognition, the coverage of the precipitation area on the route is obtained.
Step 5 (workload weights). Flight paths affected by weather are used to determine the weights of the rerouting paths. As mentioned in Sections 2.3 and 2.4, the rerouting probability of each route is calculated. Other relevant information needed for Step 5 are weather-affected arrival flight paths, route weight of workload in rerouting legs, and weather-affected dynamic route networks in the terminal area.
Step 6 (multiobjective programming). The objective function and its optimization are established to determine probable alternative routes. The LINGO software is used to calculate the workload objective function (i.e., formula (20)) and the economic objective function (i.e., formulas (21) and (22)) successively and to obtain the optimal solution of the multiobjective programming model. If the result is empty sets, considering safety requirements of the control operation, workloads were set as the objective function of the optimal solution. Pseudocode 1 describes the rerouting optimization process.
4. Results and Discussions
In this study, we first collected the radar echo data of the terminal area in Shanghai at 9:15 a.m. on August 12, 2015, acquired the corresponding image, and then determined the area of the thunderstorm activity. Second, the route for all levels across different weather conditions is determined according to the route structure, the optimal probability value for each segment is obtained, and the pilot with an effective flight path is established. Third, the additional controller workload is calculated and compared, as well as the dynamic capacity and flow distribution. Finally, for the selection of the optimal flight path, the flow control decisions and delay analysis are determined. Moreover, the process of acquiring simulation data includes these steps.(i)We investigate actual traffic of terminal area arrival/departure routes and obtain the distribution by fitting the data to a random number of flights set simulation program flow. Meanwhile, when flight goes through some waypoints, simulation program will trigger the controller workload (generate communication time). And, then, we can calculate the workload at some time interval based on simulation data.(ii)In order to ensure the accuracy of the simulation data, we take the terminal airspace capacity constraints into account, to obtain the weighted average of the 30 sets of simulation data in the peak hours of traffic as accurate simulation data.(iii)We analyze the correspondence between the rerouting scheme and controller workload by matching the simulation data and the weather data.
4.1. Processing of Radar Echo Image
The color attributes of the weather radar map acquired at the radar observation station are obtained. The map is converted into a binary image and the contour lines of the precipitation areas of , , and levels are acquired, as shown in Figures 7–9.
(a)
(b)
(c)
(a)
(b)
(c)
(a)
(b)
(c)
4.2. Contour Lines for Determining Whether the Route Passes through All the Thunderstorm Areas with Different Precipitation Levels
In the route network of the entire region affected by thunderstorms, the region that is most seriously influenced is the entry segment. The departure routes can simply wait on the ground during thunderstorms; thus, these routes are all off. The approach legs BELOP–XSY, TEMP4–XSY, JNZSS–ATRIP, and ATRIP–XSY in the Shanghai terminal area need to execute the rerouting path because of the thunderstorms.
The coordinate origin in pixel of PUD (−29908, 25181) is moved to the upper left corner of the radar echo image. Table 5 shows the function expressions of the legs after transforming the coordinates.
In the succeeding subsection, all the possible flight paths of the legs and their corresponding probabilities are presented. The rerouting probabilities of the different color indices are obtained using the survey statistics and shown in Table 6.
Based on the probability statistics and analysis of the rerouting on the legs of Shanghai terminal (i.e., BELOP–XSY, TEMP4–XSY, JNZSS–ATRIP, and ATRIP–XSY), the optimal probability value of each leg can be obtained, and pilots can be provided with effective flight paths.
4.3. Comparative Analysis of the Additional Controller Workload in Different Scenarios
4.3.1. Non-Weather-Related Theoretical Static Capacity
Terminal airspace capacity refers to the maximum service sorties per unit time of the terminal area. In this study, dynamic airspace capacity is defined as the capacity under the influence of dynamic factors, such as thunderstorms; otherwise, the capacity is called static capacity. Traffic surveys were utilized to obtain the following: traffic flow of the terminal area, departure route, and flight plan. The Shanghai terminal area was analyzed by implementing the DORATASK method (in a single-sector workload, 2880 seconds (3600 seconds × 0.8) for the sector capacity) and by applying the terminal area planning theory for the static capacity of workloads.
Samples from DORATASK were used for multiple linear regression analysis. The linear regression equation for the terminal area is , where is the total controller workload of the terminal area; is the total number of aircraft flying from WUXI to PUD approach; is the total number of aircraft flying from WUXI to SHA; is the number of aircraft arriving and departing by IBEG; is the total number of aircraft flying from PUD to PIKAS; is the number of aircraft departing from SHA to PIKAS; is the number of remaining aircraft arriving and departing. From these expressions of the regression equation and by referring to the traffic survey for the flow ratio on route, the static capacity of the Shanghai terminal area is 181 sorties.
4.3.2. Weather-Related Theoretical Dynamic Capacity
A traffic simulation is conducted to determine the controller workload and the time consumed for each flight path. The contour lines of the radar echo images of the precipitation areas of different levels show that the main approach procedures affected by convective weather conditions in the terminal areas of Shanghai are BELOP–XSY, TEMP4–XSY, JNZSS–ATRIP, and ATRIP–XSY. For BELOP–XSY, JNZSS–ATRIP, and ATRIP–XSY, the route points of their approach are covered. Thus, priority is given to finding the alternative to the approach procedure TEMP4–XSY. Seven alternative approaches based on the originally programmed segment can be derived based on the flight path of procedure, as shown in Figure 10.
(a) First flight path
(b) Second flight path
(c) Third flight path
(d) Fourth flight path
(e) Fifth flight path
(f) Sixth flight path
(g) Seventh flight path
The corresponding cost of each flight path can be obtained using the 0-1 hybrid model mentioned in the current paper and the cost calculation in the literature [41]. The alternative approach procedures are compared, and the results are shown in Table 7.
The target models for the additional controller workload for the departure and approach routes in the terminal area and the operating economic coefficient of the aircraft arewhere = 0.298, = 397.6, and represents the aircraft that flies from TEMP–XSY. When flight path 1 is executed, the variable value is 1; when the approach leg is implemented in the alternative segment, the variable value is 1.
Based on the results of the LINGO software, if the economic value of flight time is selected as the first objective function, then the first flight scheme is the optimal solution; if the additional controller workload is selected as the first objective function, then the fifth flight plan is the overall optimal solution (Table 7). Setting the additional controller workload as the first objective function results in the fifth flight scheme being affected by the thunderstorm because the controller workload is an important indicator for air traffic safety. Consequently, the aircraft would be diverted to avoid thunderstorm execution measures. Combined with the programming model in Section 3.2 on the effect of dynamic weather (i.e., dynamic theory), the capacity of the terminal area is 151 vehicles. Accordingly, 30 aircraft would be affected by service flow controls and have to do air-waiting outside the sector or wait on ground at the airport for the pretactical flow control strategy.
4.4. Comparative Analysis of the Dynamic Capacities and Flow Distributions in Different Scenarios
4.4.1. Comparative Analysis of Dynamic Capacity
Static delay was simulated by a computer. In the paper, static delay is characterized by the static topology in the terminal area of the total delay time in each air corridor. Meanwhile, dynamic delay is the total delay time in each air corridor; it was also determined by computer simulation. Finally, theoretical delay is the total delay time in each air corridor based on the dynamic terminal area topology (see Section 3.2).
Dynamic capacity was obtained from the simulated calculation of the dynamic route of the terminal airspace at peak hours. Consequently, dynamic simulation delay was determined by comparing the airspace static capacity from DORATASK but without the effects of weather. Conversely, static simulation delay was obtained by using static airspace simulation. The dynamic route planning model was compared with the static capacity of the airspace to determine the dynamic theory delay. Figure 11 shows the three delay analyses.
The DORATASK method is used to determine the load limit of the terminal sector according to the simulation data. A multivariate linear regression analysis was performed using the population samples. The dynamic capacities of the terminal areas in Shanghai in the different alternative flight paths are shown in Figure 12.
A computer simulation is conducted to compare a variety of flight paths with respect to the static capacity distribution and the theoretical capacity distribution.
4.4.2. Comparative Analysis of Flows
Based on the distributions of the aircraft flow along the departure and approach corridors, the number of aircraft waiting in the air and the number of aircraft waiting on the ground along the approach corridor can be calculated according to the simulation data. The approach and departure program codes are presented in Table 8.
In Figure 13, the blue area represents the waiting aircraft corresponding to the departure procedure; the red area represents the number of waiting aircraft in the approach procedure. The additional workload caused by the fifth flight path is less, its dynamic capacity is larger, and the planes waiting along the corridor and in the airport are less. This optimal condition further analyzes the terminal area route network.
(a) First flight path
(b) Second flight path
(c) Third flight path
(d) Fourth flight path
(e) Fifth flight path
(f) Sixth flight path
(g) Seventh flight path
4.5. Analysis of Flow Control Decisions and Delays in the Optimal Flight Path
Theoretical calculation and computer simulation data analysis are conducted to obtain approach and departure route capacity for the fifth flight path (Tables 9 and 10), time interval at approach corridors (Table 11), and departure time interval and number of waiting flights at airports and corridors (Tables 12 and 13). Subsequently, the flight delays along various departure and approach corridors can be determined (Table 14). Furthermore, the dynamic and the static simulation delay curves of air traffic can be generated (Figure 14).
In Figure 14, the static simulation delay curve of the aircraft is represented by the red line, whereas the dynamic simulation delay curve is represented by the blue line. The variation trends of the curves show that, under different conditions, the topological structure of the route in the terminal area has different degrees of delays. Compared with the static topology, when the sector in terminal area is affected by a thunderstorm, the delays in the aircraft are slightly larger, and the arrival rate of flights are low.
The following conclusions can be drawn from the theoretical calculation and simulation analysis of the traffic data in the terminal areas of Shanghai.
(1) Along the corridor WUXI (WUX), numerous aircraft are in service in a unit time. Therefore, the aircraft has to be considerably delayed along the corridor. When the dynamic capacity is predicted using theoretical methods, air traffic flow control measures are not implemented for the aircraft in WUX. Therefore, the total value of delays is basically consistent with the value of the static simulation delays determined by the theoretical capacity. In the dynamic delay simulation, the aircraft approaching WUX needs to reduce the longitudinal time interval to improve its arrival rate, and the total value of dynamic simulation delays is lower than that of the static simulation delays.
(2) Along the corridor Andong (AND), numerous aircraft are in service in a unit time. Therefore, the aircraft has to be considerably delayed along the corridor. When the dynamic capacity is predicted using theoretical methods, radar vectoring is applied to the aircraft in AND. The longitudinal time interval is longer, and the total value of the delays is higher. In the dynamic delay simulation, the aircraft approaching AND has to follow the provisions of the controllers. Part of the flight flow is distributed to the other corridors; as a result, the simulated value of the total delays is lower than the theoretical value of the delays but remains higher than the delays in the static conditions.
(3) Along the corridor Ning Po (NGB), numerous aircraft are in service in a unit time. Therefore, under the condition of static weather, the aircraft will not be delayed. When the dynamic capacity is predicted based on theoretical methods, shutdown measures should be implemented on the approach procedure of the aircraft in NGB. The aircraft will be hovering outside the sector, thereby causing considerable delays. In the dynamic delay simulation, the aircraft approaching NGB needs to follow the provisions of the controllers. Part of the flight flow is distributed to the other corridors; as a result, the simulated value of the total delays is lower than the theoretical value of the delays.
(4) Along the corridor DUMET, numerous planes are in service in a unit time. Therefore, under the condition of static weather, the aircraft will not be delayed. When the dynamic capacity is predicted based on theoretical methods, shutdown measures should be implemented on the approach procedure of the aircraft in DUMET. The aircraft will be hovering outside the sector, thereby causing considerable delays. In the dynamic delay simulation, the aircraft approaching DUMET needs to follow the provisions of the controllers. Part of the flight flow will be distributed to the other corridors; as a result, the simulated value of the total delays is lower than the theoretical value of the delays.
(5) In the Shanghai Pudong International Airport (PUD), numerous planes are in service in a unit time. Therefore, the aircraft has to be considerably delayed along the corridor. When the dynamic capacity is predicted based on theoretical methods, shutdown measures should be implemented on the departure procedure of the aircraft in PUD. The aircraft will wait on the ground in PUD, thereby causing considerable delays. In the dynamic delay simulation, the aircraft leaving PUD needs to follow the provisions of the controllers. The flight distance of the aircraft can be reduced; thus, the simulated value of the total delays is lower than the theoretical value of the delays, but it remains higher than the delays under static conditions.
(6) When the dynamic capacity is predicted using theoretical methods, shutdown measures should be implemented on the departure procedure of the aircraft in Shanghai Hongqiao International Airport (SHA). The aircraft will wait on the ground in SHA, thereby causing considerable delays. In the dynamic delay simulation, the aircraft which leaves SHA needs to follow the provisions of the controllers. The flight distance of the aircraft can be reduced; as a result, the simulated value of the total delays is lower than the theoretical value of the delays.
5. Conclusions
As the main factor of the flight delays, rerouting problems in terminal areas under convective weather conditions and airspace capacity prediction have become a topic of interest in the area of air traffic management. In this study, a rerouting planning model for terminal areas under thunderstorm conditions is established to provide a reasonable scientific basis for the flexible use of terminal airspace.
Instead of the rerouting planning model for solving route network design problems in the flight phase, the rerouting planning model for 0-1 hybrid dynamic route network in a terminal area is adopted in this study. A multiobjective programming model for additional fuel consumption, delay loss, time cost, and route length is established and the optimal planning scheme is identified by considering the changes in the speed and altitude of the aircraft in the terminal airspace when the route is under convective weather conditions. In the flight phase, the properties of the starting and end points of the route need not be considered in the rerouting problems. The proposed model handles the departure leg arc separately to solve the one-way problem of the route network in terminal areas because of the particularity of the departure leg.
A promising research direction could focus on path selection based on three dimensions. Furthermore, the procedure for optimal path selection takes the shortest set of paths as the feasible region on thunderstorm weather condition. Subsequently, the optimal solution is selected according to the safety and economic analysis. Another future research direction could be further analyzed by considering the departure procedures for ground holding and using the data from three-dimensional thunderstorm weather.
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Appendix
Extract codes of green contour in the radar echo images based on MATLAB are as follows.
“Green image recognition codes” Image=imread(C:∖old.jpg); Gray=rgb2gray(Image); R=Image(:,:,1); G=Image(:,:,2); B=Image(:,:,3); diff_R=0; diff_G=0; diff_B=0; Image_R=Image; RP_R=Image(:,:,1); RP_G=Image(:,:,2); RP_B=Image(:,:,3); XYR=~((G-R)>diff_G&(G-B)>diff_G); Mask=Gray(XYR); RP_R(XYR)=0; RP_G(XYR)=0; RP_B(XYR)=0; Image_R(:,:,1)=RP_R; Image_R(:,:,2)=RP_G; Image_R(:,:,3)=RP_B; imwrite(Image_R, C:∖Green image recognition.jpg,jpg);
“Code extraction of green binary image” Image=imread(C:∖old.jpg); Gray=rgb2gray(Image); R=Image(:,:,1); G=Image(:,:,2); B=Image(:,:,3); R=Image(:,:,1); G=Image(:,:,2); B=Image(:,:,3); diff_R=0; diff_G=0; diff_B=0; Image_R=Image; RP_R=Image(:,:,1); RP_G=Image(:,:,2); RP_B=Image(:,:,3); XYR=~((G-R)>diff_G&(G-B)>diff_G); Mask=Gray(XYR); RP_R(XYR)=0; RP_G(XYR)=0; RP_B(XYR)=0; Image_R(:,:,1)=RP_R; Image_R(:,:,2)=RP_G; Image_R(:,:,3)=RP_B; imwrite(XYR,C:∖green binary image.jpg,jpg);
“Code extraction of green contour” Image=imread(C:∖old.jpg); Gray=rgb2gray(Image); R=Image(:,:,1); G=Image(:,:,2); B=Image(:,:,3); R=Image(:,:,1); G=Image(:,:,2); B=Image(:,:,3); diff_R=0; diff_G=0; diff_B=0; Image_R=Image; RP_R=Image(:,:,1); RP_G=Image(:,:,2); RP_B=Image(:,:,3); XYR=~((G-R)>diff_G&(G-B)>diff_G); Mask=Gray(XYR); RP_R(XYR)=0; RP_G(XYR)=0; RP_B(XYR)=0; Image_R(:,:,1)=RP_R; Image_R(:,:,2)=RP_G; Image_R(:,:,3)=RP_B; A=im2bw(Image_R); B=bwfill(A,holes); BW1 = edge(A,sobel); imwrite(BW1,C:∖green contour.jpg,jpg);
Notations
Parameters in the 0-1 Hybrid Model for Rerouting Path| : | A set of corresponding flight paths for the flight phase () |
| : | Set of airports in the terminal area |
| : | Set of starting points and end points |
| : | Set of aircraft |
| : | Set of aircraft flows |
| : | Set of waypoints |
| : | Set of flight phase arcs |
| : | Set of the initial nodes in the flight phase |
| : | Set of initial flight phase arcs |
| : | Starting point of the flight program in the sector of the terminal area |
| : | Additional workload of the controllers who fly the execution route |
| : | Loss rate of aircraft when flying to node with respect to the original time point |
| : | Loss rate of the route in flying |
| : | Airport capacity in the terminal area at time node |
| : | Starting point of the flight phase arc |
| : | End point of the flight phase arc |
| : | Time of flight for aircraft flying from point to point |
| : | Planned departure time |
| : | Planned approach time |
| : | Flight time interval between adjacent aircraft |
| : | Flight time interval for aircraft flying from node to node |
| : | Planned time node of to reach node |
| : | Probability of aircraft performing flight instructions during the execution of flight path ; |
| : | End point of the flight program in the sector of the terminal area |
| : | Time delay for flying to node |
| : | Earning rate of aircraft when flying to node with respect to the original time point |
| : | Earning rate of route in flying |
| : | Capacity of sector airspace in the terminal area at time node . |
| : | When aircraft is at time node executing flight path and flying across arc, then the value of is either 1 or 0 |
| : | When aircraft is at time node , the time saved from the actual departure time at node with respect to the initially planned time |
| : | When aircraft is at time node , the overhead time from the actual departure time at node with respect to the initially planned time. |
Conflicts of Interest
The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Acknowledgments
This study was supported by the National Natural Science Foundation of China (no. U1233101, no. U1633119, and no. 71271113) and the Fundamental Research Funds for the Central Universities of China (no. NS2016062).