Probabilistic Conflict Detection Using Heteroscedastic Gaussian Process and Bayesian Optimization

Conflict detection plays a crucial role in ensuring flight safety and efficiency and is a critical component of an air traffic control system. Despite the availability of tools to support air traffic controllers in identifying potential conflicts, their quality, and accuracy remain limited due to the challenge of accurately accounting for uncertainty when predicting flight trajectories. To tackle this issue, researchers have explored various studies focused on using probabilistic techniques to model aircraft dynamics and trajectory uncertainty. However, these approaches share several common shortcomings, including their assumptions about uncertainty distributions and the high computational costs of detecting and calculating the risk of conflicts. In response to these challenges, we propose a data-driven approach combining a multi-output generative model with a Bayesian Optimization algorithm to effectively model the uncertainty of aircraft trajectories and rapidly identify the probability of a conflict. Our approach employs the Heteroscedastic Gaussian Process to capture complex trajectory patterns and uncertainty from historical data directly. The proposed predictive model can effectively capture heteroscedastic noise from real data, leading to improved predictions. It achieves Kullback-Leibler divergence around 1 to 1.3 for all dimensions which reduces by >45% for latitude, >24% for longitude, and 4% for altitude compared to the classical homoscedastic GP approach. The method also boasts high-performance predictions for 4D trajectories including descending, climbing, and en-route phases. To pinpoint when two aircraft are most likely to experience a conflict, the Bayesian Optimization algorithm is adopted, which shows good performance in terms of computational efficiency and flexibility for probabilistic conflict detection. The proposed model achieves percentage error <0.25% in estimating the conflict probability with computational cost <14s. By addressing the challenges of uncertainty and computational complexity, our method demonstrates great potential to enhance flight safety and efficiency.


I. INTRODUCTION
Conflict detection and resolution (CDR) is critical in maintaining flight safety and efficiency in an air traffic control (ATC) system.The CDR problem is typically decomposed into two sub-problems: conflict detection (CD) and conflict resolution (CR).A conflict between any two aircraft is defined as the violation of vertical and lateral The associate editor coordinating the review of this manuscript and approving it for publication was Wentao Fan .minimum separations between them in en route airspace, e.g., s v ≤ 1000 feet or 10 flight level (FL) and s h ≤ 5 nautical mile (NM) [1].In the CD problem, the primary task is to foresee potential conflicts between aircraft.As such, this task heavily relies on the prediction of the aircraft's future locations or trajectories.The stochastic nature of the air transportation environment poses challenges in predicting these future locations.This stochasticity originates from different factors such as weather, wind, unexpected airspace closure or other events, pilots' intention in free flight mode, etc.Therefore, incorporating these uncertainties during trajectory prediction and conflict detection is important to achieve accurate conflict detection.
Several approaches for conflict detection in structured as well as free flight operations have been proposed in the literature.These approaches can be classified into three main categories: deterministic, worst-case, and probabilistic.Probabilistic approaches are considered a promising direction since they can incorporate different types of uncertainties in conflict prediction.In a probabilistic setting, some studies considered empirical distribution model of future aircraft positions [2], [3] or modeling the dynamics of aircraft [4], [5], [6] with certainty.Other works include a probabilistic approach to identify conflicts at the merging waypoints [7] and the use of random walks with unknown intents to calculate the conflict probability at each moment of discretized time [8].Two main approaches in probabilistic conflict detection are the analytical approach [3] and the Monte Carlo simulation [9], [10].The analytical approach requires an actual closed-form solution which is often impractical due to the complexity of the real-world scenarios.It also requires imposing certain constraints on the aircraft dynamics and the incorporated uncertainty for the computation of conflict probability.Monte Carlo simulation is a more practical approach because it can incorporate a more complicated state-space model with flexible uncertainty modeling (e.g., non-Gaussian, multi-modal).Especially when additional information, such as weather or pilots' intent, is available, new uncertainties can be included.However, the most noticeable disadvantage of the Monte Carlo simulation is the high computational resource required to run a huge number of iterations.Yang et al. [10] propose an algorithm to speed up Monte Carlo simulation by simplifying the aircraft dynamics model.In one of the recent works, researchers have used a combination of analytical and monte-carlo simulation to address the issues of computational cost, albeit in 2-D (latitude and longitude), without significant improvement in computational speed [11].The authors in [12], [13], [14], [15], and [16] also try to detect probabilistic conflict at a discretized time by computing the conflict probability at each moment with assumptions about uncertainty and aircraft dynamics.These approaches share common limitations, such as their assumption about uncertainty distribution and high computational cost for estimating conflict probability.Moreover, the computational cost of conflict detection is also a bottleneck in the study of conflict resolution and poses a challenge for evaluating a large number of candidate resolution maneuvers.To overcome those limitations, data-driven and machine-learning approaches have been investigated for trajectory prediction [17], [18], [19] without the need for an explicit model of weather and aircraft dynamics, etc.Since there are no explicit formulae to estimate conflict probability in these learning-based approaches, an efficient method based on the Monte Carlo simulation method is a potential candidate for estimating conflict probability between aircraft.
In this study, we combine the Gaussian Process (GP) to represent the uncertainty in flight trajectories with Bayesian Optimization (BayesOpt) to locate the point with the maximum conflict probability efficiently.The use of the Gaussian process for trajectory prediction is inspired by [20].First, the Gaussian Process is used as a generative model to capture the complex trajectory pattern and its uncertainty directly from trajectories.This generative model must be able to address two primary challenges: the ability to produce multi-output values and the heteroscedasticity of data.Our work incorporates the 3-D position of the aircraft (longitude, latitude, and altitude) and, thus, is a multi-output problem.To handle this attribute, we train three GPs, one for each dimension.Thus, our generative model is, in fact, a set of three GP models.The work in [21] develops a Homoscedastic GP model to capture the consistent variance in data.In our work, the variations in longitude, latitude, or altitude of aircraft positions vary with time.Therefore, the heteroscedastic GP [22] is adopted to capture the variance of the considered data, while the classical homoscedastic GP is used as the baseline for comparison.After the development of the generative model, the current flight information (conditional points) such as current position and reference route is provided to the model for training the trajectory prediction model with uncertainties.These uncertainties are obtained by considering data variation on conditional points using a second GP.We also apply the Bayesian optimization algorithm for probabilistic conflict detection.We aim to obtain similar results as the Closest-Point-of-Approach (CPA) but in a probabilistic manner by identifying the time at which two aircraft have the highest probability of conflict.This is called Probabilistic Closest-Point-of-Approach (P-CPA).Since conflict is a rare event, the cost of estimating conflict probability is high.Our proposal method is faster and more flexible in finding P-CPA compared to the grid-search approaches and the existing methods.The proposed approach and the main results are also reported as a part of the first author's doctoral dissertation [23].

II. ADS-B DATA
In this study, Automatic Dependent Surveillance-Broadcast (ADS-B) data is used for training and evaluating the proposed approach.An en-route sector (Sector 2E) within Kuala Lumpur FIR, managed by the Singapore Air Traffic Center (ACC) for providing air traffic service from FL120 to FL360 inclusive, is selected for further analysis.Figure 1 depicts the spatial characteristics of the selected sector.It takes about 5 minutes on average for a typical flight to cross the sector.The sector contains 8 waypoints and is crossed by 8 Air Traffic Service (ATS) Route.There is one crossing in the sector and one converging point in the south of the sector.The spatial characteristics and the airspace structure of the sector are simple, therefore the scenario of the sector can be easily identified.We used the ADS-B data for December 2016 for analysis.The original ADS-B data set is a large data set with noise and missing data points.Moreover, with the given spatial information of Sector 2E, only a subset of trajectories should be considered and investigated.Thus, some data pre-processing steps are required: 1) First, we apply a 2D spatial filtering to filter out all trajectories which do not pass through the sector.A total of 12,141 flights that passed through Sector 2E in Dec 2016 were extracted.2) A second 2D spatial filtering is applied on trajectories passing by Sector 2E to filter out data points outside of the sector.It is separated from Step 1 since the criteria for filtering can be changed in the future for different information extraction.3) To deal with missing data points, we remove all the flight trajectories that have less than three data points in the sector.Thus, the final dataset size for training and evaluation of our approach is 9,082 flights.Flight trajectories can be grouped using various metrics such as entry/exit points, maneuvers inside the sector, flight levels, etc. Different flight trajectory patterns between entry and exit points can be due to the maneuver instructions provided by the controllers or requests of pilots in response to traffic or to avoid regions of turbulent weather.From our observations and previous studies, given similar entry information, we can obtain similar 4D exit points for the flights.Thus, for en-route sectors, we propose grouping flights based on their entry and exit information to form flight routes, as shown in Figure 2. The variation, observed within each group of flights, is the uncertainty of a flight route.The final dataset contains 42 groups of flights which are utilized for training and evaluating the proposed approach.

III. OVERVIEW OF THE PROPOSED METHODOLOGY
The proposed approach (Figure 3) consists of two components: a predictive model and a Bayesian optimization model.In the first step, historical trajectories are clustered in terms of the locations of their entry to and exit from the airspace sector.Trajectories from the same cluster are considered to share a common entry-exit pair; however, due to uncertainty, their actual flight paths from the entry to the exit are varied.As discussed, the trained generative model is a set of three Gaussian processes.To perform predictions for the incoming flights, recent positions and future reference positions of the flights are provided as conditional points.Given these conditional points, the dataset is filtered to identify a subset of trajectories that pass through all the conditional points.We then train our GP predictive models using filtered trajectories and learned flight position uncertainty.
The current work is focused on two aircraft conflict detection.Thus, predictive models are trained to predict the future locations of each aircraft and to be used for estimating their conflict probability.This conflict detection problem is formulated as an optimization problem, in which we identify a timestamp t when the conflict probability is the highest.The proposed Bayesian optimization algorithm is discussed in more detail in section V.For a given a timestamp t, we estimate conflict probability P(t) using the Monte Carlo simulation method.This result is stored in the memory and also used to update a probabilistic model, called the surrogate.The surrogate function is an approximation function of the original function but needs less time to evaluate.Then the acquisition function, a mathematical technique for guiding the exploration of the parameter space during Bayesian optimization, uses the updated surrogate function to make the prediction of the next value of t to be used for evaluation.The searching process stops when the number of iterations i reaches the maximum number of evaluations Max_Eval.
This proposed approach is designed to perform probabilistic conflict detection in the short and medium term with a time horizon varying between 5 and 20 minutes.The predicted conflict probability can be provided to air traffic controllers (ATCO) to improve their situation awareness.A binary conflict detection can also extracted from this probability by providing a threshold.Since air traffic control and air traffic are dynamic, the appropriate threshold should be adjustable and depending on each ATCO, which will be investigated in our future work for conflict detection with human input.Moreover, in this approach, the predictive models learn the heteroscedastic noises from the historical realization of the planned trajectory (e.g., conditional points), sufficient data is required to provide accurate estimation.In situations where the planned trajectory is non-nominal, the default uncertainty distributions can be used to generate the data for the input of this approach.In that case, the performance of this approach approximates the accuracy of the default uncertainty distributions which can be adopted from state-of-the-art studies for trajectory uncertainty modeling.

IV. PREDICTIVE MODEL FOR FLIGHT TRAJECTORY WITH LEARNED UNCERTAINTY A. PROBABILISTIC TRAJECTORY PREDICTION USING GAUSSIAN PROCESS
Recently, machine learning methods such as hidden Markov models and random forests have been used for trajectory  prediction.However, when considering a predictive model with uncertainty, the Gaussian process is a promising approach because it can effectively model the dispersion of a trajectory.The Gaussian process is a collection of random variables, any finite number of which have (consistent) joint Gaussian distributions [24].A Gaussian process f ∼ G(m(t), k(t, t ′ )) is fully specified by its mean function m(t) and covariance function k(t, t ′ ) which need to be symmetric and satisfy Mercer's condition.Since our work involves 3-dimensional positions of aircraft over time, we train 3 Gaussian processes with time as input and latitude, longitude, and altitude as outputs of Gaussian processes.

1) HOMOSCEDASTIC NOISE AS THE INPUT INDEPENDENT NOISE
We consider a single cluster with H trajectories.For each trajectory i we will have N i input vectors t paired with outputs y.Thus, the training data is of the form: D N = {(t j , y j )} j=1..N with N = H i=1 N i .Gaussian process regression is a machine learning technique for inferring likely values of y for an input t as y(t) = f (t) + ϵ where f is a Gaussian process, ϵ ∼ N (0, σ 2 ) is an input noise, and σ is a positive constant.
For the generative trajectory model, we train three one-dimensional models with t as the time and y as either the longitude, latitude, or altitude of the aircraft at time t.In this manner, we can specify the noise variance for each dimension independently.However, one drawback of this approach is that it does not model the correlation among the noises of the three dimensions.

2) HETEROSCEDASTIC NOISE AS THE INPUT INDEPENDENT NOISE
Instead of inferring y for an input t as y(t) = f (t) + ϵ where f is a Gaussian process, ϵ ∼ N (0, σ 2 ) is the input noise with σ is a positive constant, we consider the model: This is because it can capture the different noise levels along the trajectories.Furthermore, we use kernel regression to fit the noise-level function σ (t).

3) AIRCRAFT TRAJECTORY PREDICTION GIVEN CONDITIONS POINTS
From the generative model, the mean function m * and kernel k * are obtained and used to sample K next locations of the aircraft DK = (t i , ŷi ).This could be done by using the posterior distribution: Equation ( 1) can also be used to sample trajectories with given conditional points by setting D M = (T M , y M ) as the set of conditional points.A conditional point is either the current position of the aircraft or reference points such as waypoints in the flight route.Instead of providing the model with the current positions only, we consider scenarios where flights have information on the reference routes from flight plans or from air traffic controllers.This setting can be used in conflict resolution where each maneuver can be considered as a reference flight route or set of conditional points and input into our model.

B. EXPERIMENT 1
This experiment is designed to evaluate the ability of the proposed probabilistic trajectory prediction model to capture the uncertainty from data using the heteroscedastic Gaussian process.The historical data will be utilized for training and evaluating the proposed approach.Its performance will be compared with a homoscedastic GP approach, adapted from [21], for accessing the advantages of our approach in modeling the heteroscedastic noise in real trajectory data.
The training process of the probabilistic trajectory prediction model is designed and implemented following Algorithm 1.The first part of the algorithm involves training the Gaussian processes to model the trajectories in each flight group with corresponding uncertainties.This step trains one Gaussian process for each data dimension (longitude, latitude, and altitude).In the second part, the prediction Gaussian process models are trained using the trained kernel of previous models as initial kernels and conditional points as input.The algorithm requires the following inputs for predictions: processed trajectories for each flight route (grouped trajectories), minimum number of trajectories (N T ), maximum number of data points for a subset of trajectories (NP MAX ), number of conditional points (N CP ) and initial radius value of each dimension for filtering trajectories with conditional points ( ).Here, N T is used to guarantee sufficient input data for training.However, due to the limitation of computational power, NP MAX can be adjusted to generate a smaller subset of data for training the three Gaussian processes (high complexity O(n 3 )).Furthermore, is crucial for building training data for the predictive model.It is multiplied with normalized learned variance (Normalized_V ) to compute the filtering region (within radius Radius Filtered ) for each conditional point.Only trajectories that pass through all filtering regions of conditional points remain in the new data set, which is then used to train predictive models.These input parameters can be modified depending on the particular application.In this experiment, several trials have been conducted for the appropriate parameters, and as the results, the selected parameters are [N T = 10,NP MAX = 500, N CP = 4, = [1NM , 1NM , 0.5FL]].We use GaussianProcessRegression from SKlearn for learning and in the case of heteroscedastic noise, a heteroscedastic kernel, adopted from [25], is used.
The historical trajectory data, described in Section II, are used with the train/test ratio of 70/30 for each flight group.The training data is used to train the generative models which capture the uncertainty of trajectories within each group.Noting that the selected training data of each group is sufficient for training purposes with more than N T = 10 trajectories and NP MAX = 500 data points.For the testing data, the conditional points are obtained by extracting the turning points from each trajectory using the Douglas-Peucker algorithm [26].A set of conditional points for each trajectory will include the start and end locations and the set of extracted turning points.The size of those sets may vary between 3 to 5 in our dataset, depending on For comparison, a baseline probabilistic trajectory prediction model using classical homoscedastic GP adapted from [21] is implemented, trained, and tested with the same set of data.The performances of those models are accessed both qualitatively (e.g., visualizing and observing the variation of the predicted trajectories from both approaches) and quantitatively (e.g., comparing the Kullback-Leibler (KL) divergence obtained from both approaches).The Kullback-Leibler (KL) divergence, also known as relative entropy, is a measure of how one probability distribution diverges from another.It is used in this experiment to quantify the difference between the probability distributions of the prediction model and the historical data.The formula for KL divergence between two probability distributions P and Q is defined as: where D KL (P||Q) represents the KL divergence from distribution Q to distribution P; P(i) and Q(i) are the probabilities of an aircraft position i occurring in distributions P and Q, respectively.

C. EXPERIMENTAL RESULTS AND DISCUSSION
Using three models to handle multi-outputs may lead to issues with correlations among different dimensions.However, Figure 4 shows our prediction result for trajectory in 2 dimensions (longitude and latitude).In which, the mean predicted trajectories (green lines) of both models are smooth and fit to the trajectory data.The results show that combining predictive values from our trained models can still capture the distribution of multi-dimensional trajectory data.Figure 5 illustrates the results of the generative models for three dimensions (longitude, latitude, and altitude) for one flight route.The figures in the first column present the normalized values of the trajectories in each dimension over normalized time.Based on algorithm 1, the sampling step is applied, i.e. keeping 10% of data ≈ 500 data points.This will speed up the training of the Gaussian process and is also similar to the case where the number of data points is limited for fewer traffic routes.Two columns with black scatter plots on the right-hand side show the learning results.As we can see in columns 3 and 4 of Figure 5, the Gaussian process with homoscedastic kernel results is consistent in its variance over time for the three dimensions while the one with a heteroscedastic kernel replicates the distribution of raw data very well.Additionally, we can observe a bigger dispersion in latitude dimension than the others which is also reflected in the differences between the shapes of simulated data for both types of models.The simulated latitude values of the heteroscedastic model have a very small variance at the beginning, which significantly increases over time.
The results for predictive models using the heteroscedastic model for ADS-B data are illustrated in Figure 6.The second column of figures shows the results of the filtering process.Because the filtered radius will be proportional to the variance of the group at given time stamps, the filtered trajectories can also have big dispersion by keeping the trajectory quite far from conditional points.The results for predictive models can be observed in the third column of Figure 6 are interesting.The mean predicted trajectory can be used to anticipate the flight trajectory that flows the reference route (through conditional points) since it fits well with the shape of the trajectories in the data.Another interesting observation is the learned variance or uncertainty of the model.The sampling data from the model has shown its ability to capture not only the dispersion of filtered trajectories but also the full data.Moreover, the model with a heteroscedastic kernel 109346 VOLUME 11, 2023 Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.provides good results in capturing the uncertainty from trajectory data as compared to the homoscedastic model.The comparison of two types of models using Kullback-Leibler (KL) divergence is represented in Figure 7. From the left to the right, we can observe the results for different dimensions (longitude, latitude, and altitude).The heteroscedastic model achieves robust performance with KL divergence around 1 to 1.3 for all dimensions.The performance of the two altitude predictive models is similar with the difference ≈ 4%, even though the median value of the heteroscedastic model is slightly smaller.A possible explanation is that the variance of altitude dimension in raw data is consistent over time which eliminates the benefits of our proposed approach.Our proposed approach outperforms the classical Gaussian process in modeling both longitude (> 24% reduction) and latitude (> 45% reduction), especially since there is a significant improvement in the performance of the model for latitude.We can conclude that heteroscedastic modeling can significantly better capture model trajectory data than homoscedastic modeling.

Bayesian Optimization is a class of machine-learning-based optimization methods that focuses on solving the problem max t∈A F(t),
where A is the feasible set and F is a ''black-box'' objective function or is costly to be evaluated.
In contrast to GS, the BayesOpt algorithm keeps track of the past evaluation results to update a probabilistic model which maps the given parameters to the probability of a score s from objective function F(t).The probabilistic model P(s|t) is called a ''surrogate'' function with an input parameter, t, and a predicted score s.Bayesian optimization consists of five main steps: 1) Building surrogate probability model of objective function: The surrogate function is the probability representation of the objective function.There are several different forms of the surrogate function including Gaussian processes, Random Forest regression, and Tree-structured Parzen Estimator [27] which is used in this study.Instead of directly representing P(s|t), the Tree-structured Parzen Estimator (TPE) applies the Bayes rule to get surrogate probability P(s|t): P(s|t) = P(t|s)P(s) P(t) where P(t|s) is the probability of the parameter t given the score s of the objective function and is defined as: Here, s * is the threshold value of the objective function.l(t) and g(t) are two Gaussian Mixture Models (GMM) in which l(t) is fitted using the data points in which their corresponding objective values are less than s * and g(t) is fitted using the remaining data points.2) From the surrogate function, finding the next parameter t to be evaluated using Acquisition Function: Given the surrogate function, the acquisition function A EI (t) is the criteria to find the next parameter.This study uses the most common choice of this criteria, which is the expected improvement: In this equation, t is the proposed parameter, and s is the actual value of the objective function using parameter t.The aim of this method is to maximize the Expected Improvement with respect to parameter t or find the best parameters t under the surrogate function P(s|t).
If EI s * (t) > 0, the parameter t is expected to yield a better result than the threshold value.3) Apply the new parameter to the original objective function.4) Update the surrogate model incorporating the new result.5) Repeat steps 2-4 until maximum iterations or time is reached.

B. PROBABILISTIC CLOSEST POINT OF APPROACH -PCPA
We model the uncertainty in two trajectories by two sets of Gaussian processes f 1 and f 2 which are fitted to data, with the format f (t) = (GP longitude (t), GP latitude (t), GP altitude (t)).
Our purpose is to find the time in which two aircraft have the highest probability of conflict: where , GP latitude (t)), and f 3 (t) = GP altitude (t).In this work, we set s h = 5NM and s v = 10FL.
Taking the advantage of Gaussian process, we can approximately compute P(t) by Monte Carlo simulation: where N s is the number of MC runs, A i (t) is the event As discussed air traffic conflicts are rare events.As a result, estimating the conflict probability at each time t requires a large number of simulation runs N s .We also need to evaluate the maximum value of P(t) with respect to t (refer equation ( 3)).Therefore, the computational load is very high.To overcome this situation, we apply the BayesOpt technique by choosing Tree Parzen Estimators (TPE) for the surrogate probability model and using Expected Improvement for the selection function.

C. EXPERIMENT 2
This experiment's main purpose is to evaluate BO's advantages for quickly estimating accurate conflict probability between a pair of aircraft.Algorithm 2 describes the conflict probability estimations of two flights (with uncertainty) using Bayesian optimization (BayesOpt).At a given timestamp t, we perform the sampling with size N s for both flights: Then, for each pair, the distance between two aircraft Dist(P 1 k , P 2 k ) is compared with the minimum separation allowed (constraints).The conflict probability is computed directly from the number of pairs which violate the constraints over N s .The output of these algorithms is the maximum conflict probability (CProb), the time to probabilistic closest point of approach (Time − to − PCPA), and the computational cost (Cost).The hyperopt Python package is used for the BayesOpt algorithm.Since we want to avoid any assumption about the search space for t, a uniform distribution is used to sample t.It is worth noting that the used package tries to minimize our defined objective function while our objective is to maximize conflict probability.Thus, the returned probability is multiplied by a negative one in the objective function.In BayesOpt, there are two terminating conditions: the convergence of maximum objective values or maximum number of evaluations (max_evals) is reached.In this study, we use max_evals as the terminating criterion.

Algorithm 2 Bayesian Optimization Method
Input: N, GP1_2, GP2_2, Duration, Threshold Result: CProb, Time-to-PCPA, Cost 1 def objective(t): Algorithm 3, which is based on Grid-Search (GS) is developed as the baseline for comparison.The algorithm has steps similar to Algorithm 2, but the difference is how they find the value of t for evaluation.GS algorithm evaluates all values of t in a given range.Since the time interval (R) is given, only a predefined set of values t is evaluated.This approach is mentioned and commonly used for probabilistic conflict detection, especially when finding the maximum conflict probability [12], [13], [14], [15], [16].On the other hand, BayesOpt determines a sequence of values t depending on historical evaluations.This sequence of values of t can be continuous.By employing BayesOpt, the sample size is reduced, and consequently, lower computational cost is required.For accessing the performance of those approaches, a scenario of two aircraft, illustrated in Figure 8, is generated by randomly sampling a pair of trajectories from two different groups of trajectories from historical data, mentioned in Section II.The first group, called group A, is the group of north-south flights that start near waypoint UKASAUKASA, passing by waypoint PEKLA and heading to waypoint PIBAP to exit the sector.The second group, called group B, is the group of east-west flights that start from waypoint LENDA, passing by waypoint PEKLA and heading to waypoint VMR to exit the sector.The travel time of flights in each group varies between 12 to 16 minutes.The starting time of the east-west flight is modified by adding an appropriate delay to ensure a certain level of conflict probability.In this case, the time differences at the crossing points of those two trajectories are limited to 1 minute.This is a crossing scenario where the conflict happens because of the uncertainty in the position of the aircraft at the given timestamp t.The learning process is applied and at the end, we can obtain two predictive models for two routes (called A, and B).These models are then used for evaluating BayesOpt as a probabilistic conflict detection algorithm.
Detecting aircraft conflict is a safety-critical task where the accuracy of the algorithm may have a strong impact on the decision of ATCO.In this work, a large number of MC simulations (e.g., sampling size N s from [0.5M , 1M , 2M ]) will be used to accurately estimate the potential conflict between any pairs of aircraft under high trajectory uncertainty.Besides, 5 different values of sampling interval R ([1s, 5s, 10s, 30s, 60s]) of the Grid-Search algorithm are These parameter values are increased by a factor of 2 to perform parameter tuning and the following experimental results have shown their appropriateness for investigation of our approach's performance.For the convenience of discussion, all experimental models are named with a similar format {Method}{Parameter}, e.g., GS01, GS05, or BO100.To distinguish models from the same methods for different sampling sizes, a corresponding prefix will be added to the name, e.g., 0.5M -GS01 vs. 2M -GS01.In this study, each model is evaluated 100 times to obtain the average model's performance.
The computational cost and final results for P-CPA (e.g., Time to P-CPA and maximum P-CPA) will be reported and compared between different settings.Because the exhaustive search has been performed, the results of the GS with sampling interval 1s and sampling size 2M are chosen as the baseline for comparison in terms of the Time to P-CPA and Maximum P-CPA.The difference (e.g., percentage error PE) in the conflict probability (P-PCA) between the other models' estimation V Predicted and the baseline V Baseline will be extracted and discussed to highlight the advantages and disadvantages of the proposed BO approach.

PE
The experimental is conducted as described and the average T-PCPA, P-CPA and computational cost are captured for discussion.First of all, all models can achieve the average T-PCPA values which approximate the T-PCPA at 2M-GS01.Moreover, the GS01 shows the sign of convergence at N s = 1M, in which the PE of 1M-GS01 and 2M-GS01 is just 0.05%.All BO models approximate the baseline with (N s ≥ 1M ) with the PE < 0.25%.On the contrary, all GS methods with R ≥ 5s have high PE and only achieve PE > 0.65% at N s = 2M.Since with R ≥ 5s, the GS models can only evaluate the conflict probability at those discrete points (as their predicted T-PCPAs) with the expected deviation (approximately R 2 ) from the actual T-PCPA.Thus, the fine level of the grid or the value of R has a strong impact on the accuracy of the GS models' estimation.The PEs for GS30 and GS60 are very high, e.g., 15% and 30% respectively, even at N s = 2M.Those convergence behaviors can be observed in Figure 9. Since the PE values of GS30 and GS60 are too large, they are excluded from the figure for better visibility.Furthermore, with the increase in sampling size, the uncertainties of the estimation also decrease.For comparison, the standard deviation (std) of 2M -GS01 is used as the baseline (e.g., 100%) for computing the relative percentage of the std of the other models.Figure 10 depicts the reduction in relative std or uncertainties in the estimation of the evaluated models.Expectedly, the observed reductions are proportional to the increase in the sampling size of the MC simulation.F Since the computational cost strongly depends on the sample size of MC simulation, the results, in Figure 11,  show a linear increase in terms of computational costs for all algorithms when the number of samples (N s ) is increased.When R increases, there are more values of t that should be evaluated, and hence, the cost of GS models also increases.Therefore, even though GS01 can achieve high accuracy in estimating P-CPA, its computational cost is much higher than the rest of the models.It needs more than 330 seconds to estimate the conflict probability with N s = 1M while the others just need less than 120 seconds.For the given conflict scenarios, with 20-minute trajectories, the cost for BayesOpt with 100 evaluations (BO100) is close to that of GS10 (with 120 evaluations).Moreover, the BO50 can achieve high performance with significantly less computational costs (e.g., <14s seconds at N s = 1M).The only models that execute faster than BO50 are GS30 and GS60 but with a significant trade-off in model performance.For the GS05 model, it still needs 4.5x the running time of BO50 and 2.3x the running time of BO100 with a 2.6x increase in PE.Those trade-off between PE and computational costs, at N s = 1M, is depicted in Figure 12.Since GS01 has a too-large computational cost (>320s), GS30 and GS60 have too-large PE (>13.5% and >29.5%), they are omitted from the figure for better visibility.As observed from the figure, BO50 and BO100 outperform the other models in both speed and accuracy for estimating the probability of conflict between a pair of aircraft.In short, the 1M -BO50 model can be considered as the ''best'' model with PE < 0.25% and computational cost < 14s.

VI. CONCLUSION
The results demonstrate the potential of our approach for trajectory prediction and conflict detection.Heteroscedastic Gaussian process regression can learn the uncertainty from data and perform better aircraft position prediction.Although some recent approaches try to increase the confidence of prediction by using more data and advanced machine learning techniques, the uncertainty from data still exists because of various reasons including human intervention and incomplete information.Our approach for predictive modeling can successfully capture heteroscedastic noises in trajectory data with robust performance.It achieves KL divergence around 1 to 1.3 for all dimensions which reduces by > 45% for latitude, > 24% for longitude, and 4% for altitude compared to the classical homoscedastic GP approach.Moreover, the results in experiment 2 have highlighted the benefits of Bayesian Optimization over the classical Grid-Search method in terms of speed and accuracy.The algorithm can work with continuous time, which means it can locate the conflict position without the impact of the selection of the sampling interval.The proposed approach can achieve less than 0.25% percentage error compared to the baseline with only 1/23 of the computational cost.Especially, when computing the conflict probability is expensive, the contribution of this method is more significant.
For future work, we believe that improving the multi-output problem could enhance the performance of the predictive model and speed up the training.Besides, the BayesOpt approach for probabilistic conflict detection can be generalized to work for the whole sector with multiple flights and conflicts.

FIGURE 1 .
FIGURE 1.The layout of Sector 2E, with its waypoint and ATS routes.

FIGURE 2 .
FIGURE 2. Examples of grouping trajectories based on entry and exit points.

FIGURE 3 .
FIGURE 3. Illustrating approach for probabilistic conflict detection using Bayesian optimization.

FIGURE 6 .
FIGURE 6. Prediction results from predictive model with ADS-B data.

FIGURE 8 .
FIGURE 8. Illustration of the two nominal routes (referred to as route A in black and route B in yellow) for generating the scenarios for experiment 2. Route A represents a group of north-south flights, called group A, that pass by waypoints UKASAUKASA, PEKLA and PIBAP.Route B represents a group of east-west flights, called group B, that pass by waypoints LENDA, PEKLA, and VMR.Groups A and B are groups of flights from historical data.

FIGURE 9 .
FIGURE 9. Illustration of the convergence of the P-CPA estimation with the increases in sampling size of Monte Carlo simulation (from 05.M to 2M).The PE values of GS30 and GS60 are excluded from the figure (e.g., PEs ≥ 13.5% and 29.5%, respectively).

FIGURE 10 .
FIGURE 10.Illustration of the reduction in standard deviation (std) of the model's estimation with the increases in sampling size of Monte Carlo simulation.The baseline (e.g., 100%) is the std of the Grid-Search 1s (GS01) with 2M MC simulation.

FIGURE 11 .
FIGURE 11.The computational cost (s) of all evaluated methods with the increases in sampling size.

FIGURE 12 .
FIGURE 12. Illustration of the trade-off between the computational cost to find the maximum probabilistic and the of the P-CPA estimation with the increases in sampling size of the Monte Carlo simulation (from 05.M to 2M).The PE values of GS30 and GS60 are off the chart and excluded from the figure (e.g., respectively 13.5% and 29.5%).