The demand of air traffic has grown rapidly in recent decades and will continue to increase. International Civil Aviation Organization (ICAO) recommended the planning, development, and implementation of CNS/ATM as the next generation air traffic control system at the 10th ICAO Air Navigation Conference in order to accommodate the increasing demand of air traffic and enhance safety and efficiency of air traffic management. The CNS/ATM is an abbreviation of Communication, Navigation, Surveillance and Air traffic management, which is based on the Global Navigation Satellite System (GNSS) and data link technology. In 1998, ICAO announced the Global Air Navigation Plan for CNS/ATM to encourage development and implementation in the global area [1]. In the Republic of Korea, Ministry of Land, Transport and Maritime Affairs have established and initiated the Air Transportation Advancement Program to answer the ICAO’s recommendation for the development and implementation of CNS/ATM.

Automatic Dependent Surveillance-Broadcast (ADS-B) is a key component of CNS/ATM. ADS-B is new surveillance technology for precise aircraft tracking based on GNSS and data exchange technology. It allows the aircraft or ground vehicle to broadcast its identification, position and velocity states, and operational information to other aircraft, ground vehicles and ground stations in the nearby region. ADS-B will enables cockpits display of traffic information, airborne conflict avoidance, conflict management and airspace deconfliction, air traffic service conformance monitoring, and other applications. In this paper, we propose ADS-B based trajectory estimation, the prediction algorithm, and the probabilistic conflict detection algorithm [2].

Aircraft tracking is a basic component of a surveillance system, which determines the position of aircraft from radar or ADS-B measurements. The α - β and α - β - r algorithms for aircraft tracking have been addressed in [3]. The multiple model estimation based algorithms considering the maneuvering of an aircraft have been proposed in [4, 5]. Trajectory prediction can be classified into nominal, worst-case, probabilistic approaches [6]. The nominal approach predicts the future position of an aircraft by propagating estimated states along a single trajectory [7]. The probabilistic approach deals with the potential changes in the future trajectory of an aircraft by the modeling and propagation of uncertainties [8]. The worst-case approach predicts the future trajectory under the assumption that the aircraft will move under worst-case maneuver defined by applications [6]. The probabilistic approach can reduce undetected conflict detection of the nominal approach and excessive false alarm of the worst-case approach. For such reasons, it is taken in this paper.

Based on the predicted trajectory of the aircraft by the probabilistic approach, one can compute the conflict probability between two aircraft encountering. The direct numerical integration algorithm can provide nearly an exact solution of conflict probability. However, it is intractable to implement a real-time application since it is computationally intensive. Paielli et al. have proposed an approximated analytic solution by extending the conflict zone under the assumption that the aircraft maintains constant velocity [9]. Hwang et al. have refined Paielli’s algorithm using the one block approximated conflict zone by coordinate transformation [10]. The Monte Carlo simulation approach to compute the conflict probability has been proposed in [11]. Daalen et al. have proposed a fast algorithm using the flow rate of conflict probability at the boundary of the conflict zone [12].

The contribution of this paper is two-fold. First, we address the multiple-model based trajectory prediction algorithm that leads to an accurate predicted conflict probability at a future forecast time. Second, we propose an efficient and accurate algorithm to calculate conflict probability based on the approximation of a conflict zone by a set of blocks. The performance of the proposed algorithms is demonstrated by a numerical simulation of two aircraft encounter scenarios.

The remainder of this paper is organized as follows. In Section 2, we briefly address the ADS-B system and information provided by ADS-B. Section 3 presents the trajectory estimation and the prediction algorithm using the multiple model estimation approach. In Section 4, we propose an efficient numerical algorithm to compute conflict probability. In Section 5, numerical simulations are presented to show the performance of the proposed algorithms.

ADS-B is an important component of CNS/ATM proposed by ICAO. It is a precise surveillance technique developed as alternatives of current surveillance techniques such as voice position reporting, primary surveillance radar (PSR), and secondary surveillance radar (SSR). The surveillance system by voice position reporting of the pilot, called procedural air traffic control, is slow and cumbersome. It requires highly skilled human resources and is prone to be corrupted by human error. The procedural air traffic control is a dependent surveillance system that requires pilot’s reporting of position information using voice channels. The PSR is non-cooperative and independent surveillance technique that requires no action from the aircraft. The PSR measures the aircraft’s range and azimuth angle by the returned signal radiated at the ground-based radar station. It eliminates the need for aircraft’s onboard equipment. The PSR is commonly used for aircraft tracking in terminal areas. The SSR is a partially cooperative and partially independent surveillance technique that uses ground-based radar but requires the aircraft’s onboard equipment to provide its own altitude and identity information. The SSR is commonly used in en route airspace [2].

ADS-B is developed to enhance situational awareness, conflict detection, conflict avoidance, runway incursion avoidance and surveillance in a non-radar environment. ADS-B is automatic since it operates without intervention of the pilot or air traffic controller. It is a dependent surveillance technique that requires a Global Navigation Satellite System (GNSS) signal or a Flight Management System (FMS) for positional data. ADS-B provides continual broadcast of the aircraft’s identity, position, velocity, and other information to other aircrafts, ground vehicles and ground stations. In contrast with PSR/SSR, ADS-B ground-based radio stations are simple and inexpensive to install and maintain. The ADS-B can provide more frequent updates and more accurate position information compared to other surveillance systems. Thus, ADS-B can increase safety, capacity and efficiency in air traffic control [2].

Two functional operations of ADS-B are ADS-B OUT and ADS-B IN. ADS-B OUT generates ADS-B messages containing positional data of the aircraft and transmits ADS-B messages to other aircrafts, ground vehicles, and ground stations. ADS-B IN receives and displays the ADS-B OUT information from another aircraft. ADS-IN also receives air traffic information and weather information services provided by the ground stations. The ADS-B system consists of three functional components: transmitter aircraft, receiver aircraft and ground station as shown Fig 1.

The transmitter aircraft has ADS-B OUT functionality and the receiver aircraft has ADS-B IN functionality. The ground station provides traffic information obtained by collecting ADS-B messages from aircrafts in the surveillance area. Data exchange in three components of ADS-B system entails standard reports, message format and transmission protocols.

ADS-B has three reports types: Surveillance State Vector (SV) Report, Mode-Status (MS) Report, Various On- Condition Reports [2].

The SV report contains information about current kinematic states of the aircraft, which includes geometric position, horizontal and vertical velocity, and heading. Measures of the SV element quality are contained in the MS report.

The MS report includes current operational information of transmitting participant, which contains Capability Class (CC) codes, Operational Mode (OM) codes, and SV report quality. CC codes indicates the capability of a transmitting ADS-B participant, which are TCAS/ACAS installation, CDTI display capability, Air Referenced Velocity(ARV) report capability and Target Status (TS) and Trajectory Change (TC) report capability. OM codes signify the current operational mode of a transmitting participant. The SV quality codes contain the Navigation Accuracy Category for position (NACp), Navigation Accuracy Category for velocity (NACv) and Surveillance Integrity Level (SIL). The NACp provides 95% horizontal and vertical accuracy bounds for reported geometric position in the SV report. The NACp helps surveillance applications to decide whether the reported position information is acceptable for the intended use. The NACv provides 95% horizontal velocity error and 95% vertical geometric velocity errors. The SIL is defined as the probability of the integrity containment the radius of NIC being exceeded without detection.

The Various On-Condition reports include the ARV report, TS report and TC report. The details of these reports are omitted here since we did not deal with them in this paper.

All ADS-B participants must provide the SV and MS report. The minimum requirements for ADS-B data exchange for equipage classes are provided in [2]. In this paper, we propose the ADS-B based conflict detection and avoidance algorithms between two encountering aircraft using information contained in the SV and MS reports.

In this section, we describe the estimation and prediction algorithm of aircraft trajectory based on the multiple model estimation algorithm.

The dynamic models described below are used by multiple model estimation to compute the current state and flight mode estimates of the aircraft. The motion of the aircraft can be decoupled into horizontal and vertical motion. The dynamic models correspond to the different flight modes that commonly occur in real flight. The modes in the horizontal and vertical plane chosen in this paper are as follows.

1) Horizontal modes

a) CV : constant velocity mode; straight cruise flight

b) CA : constant acceleration mode; maneuvers for heading changes with constant turn rate or speed changes with constant acceleration

2) Vertical modes

a) CH : constant altitude mode; steady level flight

b) CAD : constant altitude change mode; maneuvers for ascending or descending with constant altitude change rates (e.g. takeoff and landing)

The state of dynamic models in horizontal and vertical planes are defined as

The aircraft dynamics for each mode are represented by a discrete-time linear dynamic system as:

where j ∈ {CV,CA] is the mode in the horizontal plane or j ∈ {CH,CAD] in the vertical plane. The process noise is assumed to be a zero-mean Gaussian random noise with a mode-dependent covariance matrix Q_j. The measurement noise is assumed to be a zero-mean Gaussian random noise with mode-independent covariance matrix R.

Let us define the following matrices useful for notational simplicity.

The mode-dependent system matrices and process noise covariance matrices for each mode are defined as follows.

CV mode:

CA mode:

CH mode:

CAD mode:

The ADS-B message includes position and speed information in the horizontal and vertical planes. The modeindependent output matrix is given by

The measurement covariance matrix is obtained from the ADS-B system performance parameters such as NACp and NACv in the MS report.

The flight mode dynamics are modeled as a first order Markov chain. The probabilities of the switching from mode j at time k to mode at time k+1 are assumed to be known and denoted as

The system described by Eqs. (2), (3) and (10) is called a Markov-jump linear system.

The interacting multiple model (IMM) algorithm is used to estimate the position and velocity of an aircraft in this paper. The IMM is a computationally efficient suboptimal estimation algorithm for multiple model estimation [4]. It consists of a bank of Kalman filters matched to each mode and yields final estimates by the fusing of estimates for all Kalman filters. The routine of IMM has four steps: interaction, mode-matched filtering, mode probability update, and combination as follows

Step 1: Interaction

The mixing probability μ_ij(k|k) is the probability that mode i was in effect at time k given that mode j is in effect at time conditioned all measurements up to time k.(Ed. Note: This is confusing. Please clarify) It is obtained as

where μ_i(k) is a mode probability for mode i at time k and

is a normalizing constant.

The initial condition for the filter matched to mode j is calculated by the mixing of estimates for all filters at time k with the mixing probabilities obtained in Eq. (11)

Step 2: Mode-matched filtering

The mode-matched Kalman filter j yields its own estimate

and its covariance P_j(k+1) using Eq.(2) and (3).

Step 3: Mode probability update

The likelihood functions corresponding to mode is given by

where r_j and S_j are the residual and the innovation covariance matrix of the mode-matched Kalman filter j. The function

represents an n-dimensional Gaussian probability density function (PDF) with mean

and covariance Σ.

The mode probability can be updated as

where

is a normalizing constant.

Step 4: Combination

The final estimate of states and the covariance matrix of IMM are obtained as a weighted sum of the mode-matched filter

Fig 2. presents the structure of trajectory estimation in the horizontal and vertical planes using the IMM algorithm.

The future trajectory of an aircraft at forecast time T_F can be predicted from the current state estimates and current mode probabilities from the IMM algorithm. The routine of trajectory prediction is analogous with the IMM approach with the exception of measurement update in modematched filtering and mode probability update using the likelihood functions.

Step 1: Interaction (equal to step 1 of IMM)

Step 2: Mode-matched prediction

The aircraft states and covariance are propagated as

Step 3: Mode probability update

The mode probability vector

is updated as

where ？ is state transition matrix with elements π_ij.

Step 4: Combination (equal to step 4 of IMM)

Repeat steps 1-4 to propagate the aircraft state up to the forecast time T_F

In this section, we propose an efficient numerical method to calculate the conflict probability for two aircraft in a nearby region.

Let the estimated or predicted position of aircraft A and B be

Let x_AB = x_A - x_B be the relative position of two aircraft. Under assumption of independency of the trajectories of the two aircraft, the PDF of relative position is given by

where

The horizontal conflict zone of an aircraft in en route airspace is given by a circle of radius d_min=5nm by the current standards of FAA. A conflict C between two aircraft occurs if the relative position is within the conflict zone. Since the position of two aircraft has uncertainty, the conflict is expressed by probability. The conflict probability at time k is defined as

If the conflict probability is higher than the predefined threshold, the system d|e|c|l|a|res a conflict.

In most cases, there exists no analytic solution to the integral of Eq. (24). Paielli et al. have approximated an analytic algorithm to calculate this integral by modifying the conflict zone [9]. Hwang et al. have proposed a refinement of this algorithm using coordinate transformation [10]. We propose a modification of Hwang’s algorithm to provide a more accurate integral solution based on the simple numerical integration approach.

Fig 3(a). demonstrates the relative encounter geometry of two aircraft.

The relative position uncertainty expressed by error covariance P_AB is allocated to one of the aircraft, called the stochastic aircraft. The ellipse centered on the stochastic aircraft represents the error covariance P_AB in Eq. (23). The conflict zone D is centered on the other aircraft, called the reference aircraft. The conflict probability is the sum of probability regarding the relative position in the shaded circular conflict zone by the definition in Eq. (24). The elliptic covariance centered on the stochastic aircraft can be transformed to a circular covariance by coordinate transformation. The coordinate transformation matrix T_C is obtained by Cholesky decomposition as

where L is a lower triangular matrix. Note that the circular conflict zone is also transformed to an elliptical conflict zone as shown in Fig 3(b).

In the mathematical sense, this transformation converts the two dimensional Gaussian PDF

to the product of two independent one-dimensional Gaussian PDFs N₁(x_η;)N₁(x_ξ;). Hwang and et al. have proposed the additional rotational transformation to align the major axis or minor axis of the elliptic conflict zone with the axis of reference coordinates. They considered a block conflict zone D_B circumscribing the elliptic conflict zone as shown in Fig 4(a). The conflict probability by one block approximation proposed in [10] is then given by

where δ_ξ and δ_η are the relative position of the stochastic aircraft in the new ξ - η - coordinates, and Δ_ξ and Δ_η are the half-width and half-height of the block conflict zone D_B after rotational transformation. Ψ(δ;ㆍ) denotes a one-dimensional Gaussian cumulative distribution function (CDF). This algorithm provides a simple analytic solution of the integral in Eq. (24) using a well-tabulated one-dimensional Gaussian CDF. They have further extended this idea to better approximate the elliptical zone as shown in Fig 4(b).

However, the approximation of the conflict zone using circumscribing rectangular can induce excessive conflict probability compared with the true probability from the elliptic conflict zone. The error of conflict deteriorates as the error covariance of the relative position is smaller than the conflict zone.

To improve the over-approximation problem in one block approximation, we propose a set of blocks approximation as shown Fig 4(c).

The ellipse conflict zone can be approximated by a set of blocks that covers the ellipse. This set of block reduces the superfluous area at the corners of one block approximation. The conflict probability by a set of blocks approximation is then given by

where P[C(k)D_Bi] is the probability in block B_i and is calculated similar to Eq. (27).

In this section, we demonstrate the performance of the proposed ADS-B based trajectory estimation and prediction algorithm and probabilistic conflict detection algorithm through two examples regarding a scenario of two aircraft encountering.

The first encounter scenario is shown in Fig 5.

In this scenario, aircraft A flies heading toward the +x direction at (0 nm, 0 nm) and then turns left at 50s to avoid aircraft B. After avoiding aircraft B, aircraft A flies heading toward +y direction. Aircraft A maintains forward constant velocity during the whole flight and a turn rate as 1.5deg/s during the conflict avoidance turn maneuver. We assume that aircraft B is stationary at (7 nm, -2 nm) to show the performance of the proposed algorithm clearly. The position uncertainty of aircraft B is assumed as P_B=diab(1nm², 1nm²) and the position uncertainty of aircraft A is calculated by the IMM algorithm. The radius of the conflict zone is set to 2 nm. The total simulation time is 130s. The time steps for trajectory estimation and prediction, and conflict detection are set to be 1s since ADS-B messages are provided every second. The forecast time is set to 20s.

The IMM algorithm is used to estimate and predict the trajectory of aircraft A with CV and CA linear models in the horizontal plane. The mode transition probability matrix between the CV and CA mode is assumed as

The process noise covariance matrix for the CA mode is modified to apply the heading of the aircraft as

where Ψ(k) denotes the heading of the aircraft A. Fig 6. shows the mode probabilities of CV and CA modes computed by IMM algorithms.

It shows that the mode probability of the CV mode is larger than the CA mode in non-maneuver flight and the mode probability of the CA mode becomes larger during a turning maneuver. Using current estimates and mode probabilities of aircraft A, one can predict the trajectory of aircraft A after the forecast time T_F as mentioned in section 3.

Fig 7. presents the instantaneous and predicted conflict probabilities computed by the proposed set of blocks algorithm.

The instantaneous conflict probability described solid line is computed from current estimates of aircraft A. The predicted conflict probabilities are computed from the predicted estimates by three models: CV model, CA model and combined model by the IMM algorithm. The predicted conflict probabilities are illustrated by shifting +20 second in Fig 6. to clearly compare with the instantaneous probability. The predicted conflict probability by the CV model is overestimated in comparison to the instantaneous probability. It may induce an excessive false alarm of the conflict detection, which annoys pilots and air traffic managers. The predicted conflict probability by the CA model is underestimated compared to the instantaneous probability from 52s to 72s. It may induce unreported conflicts endangering aircrafts. The predicted conflict probability by the combined model through IMM is also overestimated, but its error is smaller than the CV model. The prediction by IMM can derive more correct conflict probability prediction and reduce the excessive false alarm of conflict detection.

The second encounter scenario is shown in Fig 8.

In this scenario, aircraft A flies heading toward the +x direction at (-10 nm, 0 nm) with a constant velocity of and aircraft B flies heading toward the +x direction at (-5 nm, -2 nm) with a constant velocity of . Aircraft A overtakes aircraft B since it is faster. The radius of the conflict zone is set to be 2 nm. The conflict event does not happen during overtake since the clearance between the two aircraft in the y-axis is 2 nm. To demonstrate the performance of the proposed set regarding the block approximation algorithm, we consider two cases with different navigation accuracy of aircraft : a) P=diag(1nm², 1nm²) and b) P=diag(0.1nm², 0.1nm²).

Fig 9. shows the comparison of conflict probabilities computed by the one block approximation algorithm in [10] and the proposed set of blocks algorithm with a varied number of blocks.

The true conflict probability as a reference is calculated by the direct numerical integration of Eq. (24) with 10⁶ uniformly spaced points in a circular conflict zone. To obtain an accurate integral solution, the direct numerical integration requires a functional evaluation on numerous integral grids (or points). It is computationally expensive and intractable to implement real-time applications. The one block approximation algorithm yields relatively poor conflict probability compared to the set of blocks approximation algorithm. The over-approximation at corners of the circumscribed block leads larger probability than true probability. The proposed algorithm yields accurate conflict probability by overcoming this drawback of one block approximation by a set of blocks approximation. As the number of blocks increases, the set of blocks approximation algorithm provides more accurate conflict probability as shown in Fig 9. The improvement of conflict probability calculation is more significant in the good navigation accuracy condition b) P=diag(0.1nm², 0.1nm²) than bad navigation accuracy aa) P=diag(1nm², 1nm²). If the position uncertainties of aircraft are large, the conflict probability in the superfluous area at the corners of the circumscribed block is relatively small with conflict probability in the true conflict zone. However, if the position uncertainties are small, the conflict probability in the superfluous area is relatively large. Thus, a well-approximated conflict zone is required to calculate accurate conflict probability if the aircraft has good navigation accuracy. The number of blocks remains a tuning factor for the designer. The proposed algorithm requires more computation load than the oneblock approximation algorithm in order to calculate the probability of each small block. However, it is not serious since the proposed algorithm requires simple arithmetic and loading of the value regarding one-dimensional Gaussian CDF, which is well-tabulated. The proposed algorithm is accurate and much more efficient than the direct numerical integral algorithm.

The ADS-B based trajectory estimation and prediction algorithm and conflict detection algorithm have been developed. The IMM algorithm is used to estimate the trajectory of an aircraft from information contained in ADS-B messages. The future trajectory of an aircraft is accurately propagated from current estimates and flight modes probability based on the IMM algorithm. The set of blocks approximation algorithm is proposed for the computation of conflict probability. This algorithm is more accurate than the one block approximation algorithm and much more efficient than the direct numerical integration approach. The proposed conflict detection algorithm can be suitable for conflict detection in not only air traffic control but also other application areas such as vessel traffic management.