Arrival Traffic Separations Modelling By Using Continuous Probability Distributions

This article presents a method that allows to analyze selected aspects of past arrival traffic by modelling distributions of time separations of arriving aircraft in a chosen navigation point of Terminal Manoeuvring Area with the use of continuous probability distributions. Modelling arriving aircraft time separations distribution with continuous probability density functions allows to apply various mathematical tools to analyze separations distributions. Moreover, by comparing distributions parameters, quantitative analysis of separations for days with various arrival traffic intensity can be performed. Assumptions, mathematical model, application in the exemplary experimental scenario with an airport and days with low and high traffic intensity, and results are presented in this article. Real air traffic data was used for the experimental scenario. Outcomes show that the method can be used for air traffic post-analysis, e.g assessment of maintaining separation.


INTRODUCTION
Constantly growing air traffic has a negative impact primarily on the environment, air operations safety, workload of air traffic controllers (ATC) and plane crew. Therefore, there is a need to seek for new methods which will allow for better and more robust analysis and organization of air traffic. A solution to these problems can be application of new air traffic sequencing methods. This will organize air traffic in an automatic and, at the same time, optimal way. Moreover, it will lower human workload, especially in departure and arrival phases that will result in a reduction of mental stress and its negative effects [1], [2]. Initial analysis showed that a solution for descent phase may be to combine Continuous Descend Operations (CDO) and Required Time of Arrival (RTA) assigned to each flight [3], [4]. Use of CDO allows to eliminate segments of horizontal flight. This reduces fuel consumption in aircraft descent phase, and thus lowers exhaust emission and flight cost. On the other hand RTAs allow to control trajectory in fourth dimension -time, what results in a better assessment of separations between aircraft. Therefore, it is less expected that radar-vectoring or go-around will be required due to too low separation. This concept could also contribute to the Aicraft Supervision System Concept [5] that assesses overall external and internal aircraft status. Implementation of this approach requires detailed analysis of past air traffic and study of wide range of possible scenarios, with a focus on separations between aircraft. The approach based on time separations, so-called Time-Based Separation (TBS) [6] is also required for 4D trajectories application.
Modelling of arriving aircraft time separations distribution at a given waypoint of Terminal Manoeuvring Area (TMA) for past air traffic with the use of one of known probability distribution will allow to use wide range of mathematical tools for separations analysis. When distribution parameters are obtained from archival air traffic, new test cases can be automatically generated by sampling data from obtained probability distribution.
Current studies on air traffic analysis focus mainly on assessing the complexity of a current air situation and possibility of overloading of air sectors. They also classify ways to assess the complexity of air situation and methods used for these purposes, such as graph theory, dynamic systems theory, geometric approaches, and others. A detailed overview can be found in [7] and [8]. These works also present models that allow to estimate the workload of ATC services. In [9] authors present the complex metric allowing to assess a current air situation, taking into account measurement uncertainty. The assessment of the values of aircraft separations was made in [10]. The authors used a simple classification in which separation below the required value is considered as a conflict, and otherwise a lack of a conflict. The test scenario covered the French airspace. It is worth noting that in none of these works the concept of four-dimensional trajectories [6] was used, and all of these works focus mainly on an overall assessment of a state of the airspace, often taking into account hypothetical cases, solved using sophisticated mathematical methods. This paper presents a method that allows to easily assess the air situation in the final stage of descent, in TMA, where the most intensive air traffic is.

SEPARATIONS IN AVIATION
Separations in aviation are standardized by ICAO Doc 4444 [11] which contains procedures for air traffic management. Vertical, lateral and horizontal separations can be distinguished. In descent phase, after entering TMA, every aircraft follows strictly defined trajectory, therefore lateral and vertical separations are preserved. However, there is a significant possibility of horizontal separation being too small, what is called a conflict. More precisely, conflict appears when separation between two approaching aircraft, one after another, is smaller than required minimum. Separation is determined basing on comparison of positions reported by both aircraft or on comparison of reported time of flying over given waypoint. If surveillance systems are used, such as radars, ADS-B and others, minimum separation shall not be smaller than 5 NM. This value can be slightly reduced after fulfillment of conditions described in ICAO Doc 4444. Additional requirement is separation resulting from presence of wake vortex turbulence [12] caused by preceding aircraft, which can be very dangerous for the following aircraft. Turbulence value is often assessed with the use of numerical models [13], which usually express only a part of reality, so there is still some level of inaccuracy. Therefore, in this case aircraft are classified by their maximum take-off mass.
Increasingly, TBS concept is used instead of distance to assess horizontal separations. It is caused by the fact that TBS allows to achieve better synchronization of air traffic, due to fact that aircraft have different speed values during descent. This is not addressed while analyzing distance separation, while time separations reflect this. Moreover TBS allows for more effective air traffic control in windy conditions. When wind speed is significant, wake vortex turbulence disperses faster and this physical phenomena is time-based. Finally, TBS is part of 4D trajectory concept, thus all range of features offered by this concept may be used. This results in optimized trajectories and shorter flight times.

FITTING PROBABILITY DISTRIBUTION
Values of separation between pairs of aircraft flying behind each other are obtained as a result of random process, which can be described by one of known probability distributions. Usually, gamma distribution is used to determine probability of occurrence of the same event after given amount of time. However, often other similar distributions are also used, i.e. weibull and lognormal. All of mentioned probability distributions have two parameters referring to scale and shape of the probability density function. As sampled data is limited to some number of observations, and results are often inconclusive, sometimes it is hard to assess which of these distributions should be chosen, and decision about selection of the distribution is made basing on the knowledge of the process [14].
To check if an obtained set of observations (x 1 , x 2 , ..., x n ) belongs to the population of given probability distribution f (x, φ) where φ is parameters vector to be estimated, assessment is made by distribution fitting and evaluation of those outcomes.
Fitting probability distribution consists in finding appropriate function, which describes behavior of a random variable in as good way as possible. For a random variable with known probability distribution f (x|φ) and samples (x 1 , x 2 , ..., x n ) unknown parameters vector φ has to be estimated. To estimate unknown values of distributions parameters Maximum Likelihood Estimation (MLE) method [15][16][17] was used. This method requires defining likelihood function L for data samples (1): The value of likelihood function provides information about probability of obtaining the same data sample, with an assumption that sampling will be from chosen probability distribution. Fitting probability distribution with the use of MLE can be presented on the example of the gamma probability distribution [18], defined as follows (2): where x is a random variable, α, λ are the gamma distribution parameters, respectively shape parameter and inverse scale parameter, and Γ stands for gamma function [18]. By combining (1) and (2) and denoting L = L(x 1 , x 2 , ..., x n |α, λ), likelihood function for the gamma probability distribution takes the following form (3): Estimation of unknown parameters φ can be done by finding such aφ which maximizes the likelihood function L. In practice, minimum of negative log-likelihood function is seeked (4): To solve equation (4), analytic methods, like comparing partial derivatives to zero, may be used, but for more complicated forms of the likelihood function iterative methods are selected. Initial estimate and, at the same time, a starting point for maximum likelihood method may be obtained by using method of moments. By using this method, initial estimates of shape and scale parameter can be calculated using the following equations (5): where x is the mean and σ 2 is the variance.

EXPERIMENTAL SCENARIO
Archival air traffic data was obtained from the Eurocontrol Data Demand Repository 2 (DDR2) database [19]. DDR2 contains, among others, planned and archival air traffic for European airspace. The data is stored in Aeronautical Information Regulation And Control [20] files that contain air traffic for 28 consecutive days. To obtain data from the DDR2, the Eurocontrol's software, Network Strategic Tool was used. Then, data on incoming traffic were exported. In this way, a set of two trajectories was obtained for each day: m1, the trajectory reflecting a last filled flight plan and m3, the trajectory obtained on a basis of radar measurements, determined after the flight. By comparing descent trajectories for m1 and m3, a decision was made that m3 data will be used for further analysis, as these are real trajectories that may differ significantly from planned ones, especially in the descent phase. In addition, they allow calculating the real distribution of the flight separation values. The test scenario was prepared for Barcelona El-Prat airport (Spain), ICAO code: LEBL. The airport is known as one of the most busiest airport in Europe. The airport has three runways and five approach directions. The airport chart [21] is presented in Figure 1. The direction from which the plane will approach depends on a current configuration of the airport, which depends on a direction of wind and a time of a day. For Barcelona El-Prat airport, as defined in the Standard Arrival Chart (STAR), three airport operational configurations are possible. The most commonly used is Conventional / West Configuration, which corresponds to the STAR3 descent procedure [22]. The descent procedure can be started from eleven directions, which then are grouped in four Initial Approach Fixes (IAFs), where instrument approach begins. Then, from four IAFs, the air traffic is entirely directed to the Intermediate Fix (IF) (this is the navigation point at which the aircraft finished descent and starts the ILS procedure for landing), that is shown in the Instrument Approach Chart (IAC) [23] for this configuration (Figure 2). For this configuration, the direction is RWY 25R, for which IF is TEBLA. Due to the fact that in this point all aircraft arriving at the Barcelona El-Prat airport are grouped, this point was selected for the analysis of the separation values.  Figure 2. Instrument Approach Chart for RWY 25R [23] To represent distribution fitting, two days were selected from the analyzed data set. During those days all planes followed the STAR3 procedure and the RWY 25R direction. Selected days: January 27 th , 2016, a day with low arrival traffic -229 landings took place during that day and April 3 rd , 2016, a day with high arrival traffic -380 landings took place during that day. For each pair of aircraft flying one after another, the time separation values were determined at TEBLA, by calculating absolute difference between times at which each of aircraft passed TEBLA fix. In that way time separations for whole day were obtained and were used as an input data for distribution fitting. Eighty-nine continuous probability distributions [24] were fitted to the data and distribution with the lowest value of Sum of Square of Errors (SSE) was chosen as the best fit.  Figure 3 show that only small amount of distributions among tested is close to reflect characteristics of time separations distribution. Figure 4 presents the best fitted distribution for January 27 th , 2016. The horizontal axis shows the value of time separation in seconds and the vertical axis shows probability value. The histogram bars present actual separations data, while the continuous line presents best fit probability distribution. Results presented on Figure 4 show that the distribution which fits best to the data is Johnson's S U with shape parameters a=−0.93, b=0.69, location parameter loc=102.46 and scale parameter scale=29.63. The SSE for presented fit is 1.49 · 10 −6 .   Figure 5 show that only some distributions fits were close to the data distribution. Figure 6 presents the best fitted distribution for April 3 rd , 2016. Markings are analogous to Figure 4. Results presented in Figure 6 show that the distribution which fits best to the data is Johnson's S U with shape parameters a=−1.65, b=0.83, location parameter loc=70.47 and scale parameter scale=22.41. The SSE for presented fit is 6.17 · 10 −6 .
Further analysis that consisted of ten days in which all planes followed the STAR3 procedure and the RWY 25R direction were performed. In eight of ten analyzed days the best fit was for the Johnson's S U distribution, whilst for remaining two days Johnson's S U was 2 nd and 5 th best fit (Table 1).  Linear regression was applied to obtained estimates of every of the Johnson's S U distribution parameters versus number of arrivals in a given day. Results are presented in Figure 7. In all of the inset figures obtained estimates of parameters are marked with dots and linear regression is presented with the continuous lines. Horizontal axes present number of arriving aircraft, while vertical axes present parameter values.  From the plots in Figure 7, it can be observed that shape parameters a and b, as well as a scale parameter are increasing with air traffic density, whilst a location parameter is independent of air traffic density. From the regression analysis it is also possible to determine the separation distribution parameters for any traffic intensity. However, this requires more data points (e.g in order to select polynomial order) and is a subject of a separate study.
Results showed that regardless of the volume of air traffic only few of dozens distributions can be used to fit distribution of time separations, while the best fit can be achieved with use of the Johnson's S U distribution. By comparing the results obtained for different air traffic intensities, it can be observed that shape and scale parameters change, while location parameter is approximately constant. It is due to the fact that for a smaller intensity of air traffic distributions are slightly more flattened, while the maximums of the functions do not move significantly.

CONCLUSIONS
In this article method for modelling of arrival traffic separations distribution for a chosen waypoint with the use of continuous probability distribution was presented. Realistic experimental scenario was prepared and presented method was applied, taking into account assumptions. Results showed that presented method can be used successfully to model separations distribution.
In the future more sophisticated methods may be considered to achieve more accurate models, especially for parameters regression. Moreover, creation of generalized model can be considered, which will allow to estimate distribution parameters as a function of a number of flights, airport configuration, wind speed and other arguments.