Radar tracking with an interacting multiple model and probabilistic data association filter for civil aviation applications.

The current trend of the civil aviation technology is to modernize the legacy air traffic control (ATC) system that is mainly supported by many ground based navigation aids to be the new air traffic management (ATM) system that is enabled by global positioning system (GPS) technology. Due to the low receiving power of GPS signal, it is a major concern to aviation authorities that the operation of the ATM system might experience service interruption when the GPS signal is jammed by either intentional or unintentional radio-frequency interference. To maintain the normal operation of the ATM system during the period of GPS outage, the use of the current radar system is proposed in this paper. However, the tracking performance of the current radar system could not meet the required performance of the ATM system, and an enhanced tracking algorithm, the interacting multiple model and probabilistic data association filter (IMMPDAF), is therefore developed to support the navigation and surveillance services of the ATM system. The conventional radar tracking algorithm, the nearest neighbor Kalman filter (NNKF), is used as the baseline to evaluate the proposed radar tracking algorithm, and the real flight data is used to validate the IMMPDAF algorithm. As shown in the results, the proposed IMMPDAF algorithm could enhance the tracking performance of the current aviation radar system and meets the required performance of the new ATM system. Thus, the current radar system with the IMMPDAF algorithm could be used as an alternative system to continue aviation navigation and surveillance services of the ATM system during GPS outage periods.

Introduction

According to the International Civil Aviation Organization's plan, the current technology trend is to develop and implement Communications, Navigation, Surveillance, and Air Traffic Management (CNS/ATM) systems based on the Global Positioning System (GPS) to replace legacy Air Traffic Control (ATC) systems based on ground-based radar [1]. The Next Generation Air Transportation System (NextGen) will modernize the air traffic control system in the United States, and many of the foundational elements required to meet the predicted capacity and efficiency improvements rely on widespread use of precision Positioning, Navigation, and Timing (PNT) services provided by the GPS [2]. The Federal Aviation Administration (FAA) needs to ensure a sufficient backup PNT capability to reduce the risks of aviation users if GPS becomes unavailable. As described in the FAA Flight Plan [2], the air traffic system of the future will be much more dependent on GPS. This strategic goal will not be achievable without a NextGen Alternate PNT (APNT), especially in the event of a GPS loss of service. The APNT program ensures alternate PNT services provided by the ATC system, minimizes the economic impact of GPS outages, and supports air transportation's timing needs. The existing legacy navigation and surveillance infrastructure, i.e., VHF Omnidirectional radio Range (VOR), Distance Measuring Equipment (DME), Non-Directional Beacon (NDB), and secondary radar, may not achieve the desired level of performance for NextGen operation [3].

To improve ATC tracking performance, one solution is to upgrade the ATC radar beacon system. However, this solution is expensive. A more cost-effective solution is to improve the ATC tracking algorithms, which does not require the implementation of new radar facilities. The present study thus proposes a radar tracking algorithm that meets the requirements of APNT systems. The major challenge and difficulty for ATC tracking is tracking target maneuvering in a cluttered environment [4-7]. Bar-Shalom and Blom [4] introduce a tracking algorithm called the Interacting Multiple Model (IMM) estimator, which provides tracking estimates with significant noise reduction and fast response to sequences of aircraft maneuver modes [4,5]. However, the tracking of an aircraft in a cluttered environment might be a challenge due to the several observations for a single aircraft under such environment [6,7]. That is, some tracking measurements do not originate from the target aircraft. Therefore, the present study utilizes the Probabilistic Data Association (PDA) filter [6,7] to assign weights to the validated measurements. The PDA filter can extend the tracking capability to a highly cluttered environment. In order to gain possible improvement on the tracking performance, this study combines the IMM estimator and PDA filters to create an IMMPDA filter (IMMPDAF). Most studies on radar tracking algorithms use simulated data for analysis. To further evaluate the IMMPDAF, this study uses real flight radar data collected by the Civil Aeronautics Administration (CAA) of Taiwan. A filter that uses the nearest neighbor method and the standard Kalman filter (NNKF) is commonly used for radar tracking systems. Therefore, this paper compares the tracking performance of the IMMPDAF and the NNKF in terms of positioning accuracy. Finally, the validation of the IMMPDAF under APNT requirements is discussed. Due to the complexity of the IMMPDAF, this paper also conducts a computational load study on the IMMPDAF and the NNKF.

The rest of this paper is organized as follows: Section 2 introduces the algorithms of the IMM estimator and the PDA filter. Then, the IMM estimator and the PDA filter are integrated. The dynamic model setting of the IMMPDAF for the ATC scenario is also discussed in Section 2. The implementation of the IMMPDAF to a radar system is discussed in Section 3. The tracking performance of the IMMPDAF is compared with that of the NNKF. In addition, the computation loads of these filters are evaluated. Section 4 concludes this work.

IMMPDAF

IMM Estimator

The IMM estimator is a good candidate for a radar tracking system [4]. The IMM estimator provides significant noise reduction and a fast response. Previous investigations have found that the IMM estimator is a cost-effective technique for tracking a maneuvering target [4]. In the IMM estimator, several Kalman filters are used in parallel. Each Kalman filter uses a different dynamic model. The IMM estimator can be separated into four steps. In the first step, the state estimate and covariance are mixed. These mixed estimates and covariances are the inputs of the Kalman filter. The equations are:

In Equations (1–3), the subscripts i, j, and l indicate different Kalman filters. There are up to r different dynamic models of Kalman filter, the total number for u, , , and is r, and the total number for U and p is r2. u (k − 1) represents the model probabilities in the previous time (k − 1). Since the initial model probabilities have little impacts on the results over time, they can be conveniently chosen. In this study, the initial model probabilities are set to be rectangular functions. p is an element of the Markov transition matrix; there is probability p that the target will transit from model i to model j at each time step. U represents the conditional probability of the target in the j model state, which transited from the i model state. According to [5], the IMM algorithm performance is robust to the choice of transition matrix. and are the updated state and covariance, respectively, from the Kalman filter in the previous time step (k − 1). and are the mixed state and the mixed covariance, respectively, for each Kalman filter as inputs at time k. Then, the second step implements the multiple models of the Kalman filter:

The details of Kalman filter can be found in [8]. In Equations (4,5) the subscript j indicates it is the j th Kalman filter, and Equations (4,5) are applied for each Kalman filter. In Equation (4)F is the dynamic model and Q is the covariance of the process noise. and are the uncorrected predicted state and covariance, respectively. In Equation (5)H is the observation matrix, R is the covariance of the measurement noise, and K is the Kalman gain. Z(k) is the measurement at time k and I is the identity matrix. and are the corrected predicted state and covariance, respectively. In the third step, the updated model probabilities are calculated as: In Equation (6), the subscript j indicates it is the jth Kalman filter; exp[] means the exponential operation and || represents calculating the determinant of the matrix. The likelihood function Λ is the probability density function of the estimated measurement H · , which is a normal distribution with mean Z(k) and covariance S(k). In Equation (7), the subscripts i and j indicate different Kalman filters. u(k) represents the latest model probabilities at time k. In the last step, the corrected predicted state and covariance from the Kalman filters are combined. The weighting factor is the latest model probabilities u(k):

X(k) and P(k) are the combinations of the system state and covariance from 1 to rth Kalman filters, respectively. Using Equations (1–9), the IMM estimator calculates the system state and covariance from time k − 1 to k.

PDA Filter

The PDA filter is designed for tracking a target in a cluttered environment based on the Kalman filter [6,7]. The PDA filter obtains an estimator which incorporates all the returning measurements that might originate from the target of interest rather than select only one of them. The resulting estimator sets a gate to determine whether the measurements are valid. The association probabilities of the target being tracked for each validated measurement are then calculated. These probabilities are assigned to all of the validated measurements with different weights. Next, the validated measurements are combined with different weights based on location, and then the combined measurement is used to update the state estimate of the target. The PDA filter can extend the tracking capability to a highly cluttered environment.

The probabilities are used for weighting the correction and covariance in the PDA filter. The basic assumptions and theories for the PDA filter can be found in [6,7]. Here, a brief introduction of its functions and equations is given. Because the PDA filter is based on the Kalman filter, the first equations are the same as those for the Kalman filter: Similarly, the first step is to predict the system state and covariance. In Equation (10), the subscript j indicates it is the jth PDA filter, is the calculated state at time k − 1, F is the dynamic model, and is the predicted state at time k. H is the observation matrix and zˆ(k) is the predicted measurement at time k. is the calculated covariance at time k − 1 and Q is the process noise covariance. The second step is to determine which measurements can be used to update the state and covariance:

The main purpose of Equation (11) and (12) is to validate the measurements by Chi-Square test, and it includes all the measurements returned by all the radars. In Equation (11), the total number of measurements is unknown; the subscript i indicates it is the ith measurement. is the innovation of ith measurement for jth PDA filter, z(k) is the ith measurement at time k and R is the covariance of the measurement noise. Equation (12) is the validation equation. γ is the threshold corresponding to the gate probability. γ can be obtained from Chi-Square tables for a chosen gate probability. Once the ith measurement passes the Chi-Square test in Equation (12), it can be utilized in the rest of the PDA filter. Then, the association probabilities of each validated measurement are calculated as:

In Equation (13), the subscript j indicates it is the jth PDA filter, V(k) is the volume of the validation region, m is number of total validation measurements, and N[z(k);zˆ(k);S(k)] represents the normal distribution, centered at zˆ(k)with variance S(k). is the likelihood of the ith measurement for jth PDA filter, is the associated probability of the ith measurement for jth PDA filter. In Equations (13) and (14), P and P are the target detection probability and the gate probability, respectively. β(k) is the association probability of the validated measurements. Note that β0(k) is the probability that all measurements are not valid [7]. The suggested values of P and P can be found in [5]. P is required to be larger than 90% for primary surveillance radar sensor and 98% for secondary surveillance radar sensor. P is the probability that the true measurement from the target is detected and validated; and P is usually chosen to be at least 95%. If P is set too low, the gate γ will be too small that no measurement can pass the validation and it could be end up losing track. If P is set too high on the other hand, the gate γ will be too large that causes too much false alarms (i.e., not true measurement from the target) and reduction on the accuracy [5]. The last step is to update the system state and covariance. In this step, only the validated measurements and their association probabilities are used:

In Equations (15,16), the subscript j indicates it is the jth PDA filter, v(k) is combined innovation; W(k) is Kalman gain; P″(k) and P‴ (k) are the covariances of state updated. In Equation (16)X+(k) and P+(k) are the calculated system state and covariance, respectively at time k. The propagation of the PDA filter is shown in Equations (10–16).

IMMPDA Filter

The IMM estimator and PDA filter can be combined as shown in Figure 1. The Kalman filters in the IMM estimator are replaced by the PDA filters. At the top, Equations (1–3) are implemented to deal with the data from the previous epoch. Then, the PDA filters [Equations (10–16)] use the data from Equations (1–3) as the initial values and generate the calculated states and covariances. The next step is to use Equations (6,7) to obtain the latest model probabilities. The last step is to utilize Equations (8,9) to combine the states and covariances with different weights, which are related to the latest model probabilities. The combinations of states and covariances, X(k) and P(k) are the outputs of the IMMPDAF. In Figure 1, the IMMPDAF has three different models of PDA filter. The number of models depends on the number of dynamic motions of the target. In this paper, the tracking target is a commercial airline aircraft. The selection of suitable motion models for an ATC system is discussed in the next section.

Figure 1.

Flow chart of the IMMPDAF with three models.

Tuning of IMMPDAF

As described above, the IMMPDAF can be considered as a modification of the Kalman filter. The Kalman filter is tuned by changing the process noise covariance Q matrices and the measurement noise covariance R matrices. The matrices Q and R adjust the Kalman gain. The R matrices represent the accuracy of the measurement device; larger R matrices mean that the measurements can be trusted less. In most practical situations, matrices Q and R are unknowns. In this paper, a run-time estimate method [9] is used to obtain the measurement noise covariance R. The adaptive measurement noise covariance is estimated in Equation (17). v(n)is the innovation in time epoch n, obtained from Equation (15); and v(k) is average innovation value from time 1 to time k:

From the conclusions of [9], the run-time estimate method can handle changes of measurement noise covariance, and converge faster than that of a standard Kalman filter.

Dynamic Models of Commercial Airline Aircraft

The key to successful target tracking is selecting suitable dynamic models that fit the target motion [10,11]. Ideally, each motion should have a correct dynamic model. In an ATC scenario, the motions of commercial airline aircraft can be summarized as ascent, descent, en-route, turning, and change of speed. These motions can be described as three basic dynamics: constant velocity motion, accelerated motion, and turning motion [11,12]. Following [11-15] and our previous work [16], this study uses the constant velocity model, the acceleration model, and the turning rate estimator model, respectively given in Equations (18–20), respectively, for the IMMPDAF. For ATC radar tracking, the measurement is the target 3-dimensional position (x, y, z) of the tracking target in WGS-84 coordinates. In Equations (18–20), X is the system state, x, y, z are the positions, ẋ, ẏ, ż are the velocities, and ω is the turning rate. F represents the dynamic model and t is the time interval. Generally speaking, the time epoch is close to 5 seconds. Q is the covariance of process noise and σ is zero-mean Gaussian white noise. The σ value is tuned by experience and testing, and the value is related to the measurement noise covariance:

Experimental Results and Discussion

The experiment data are from the CAA of Taiwan. The data formats include radar data (ASTERIX AC048) and ADS-B (ASTERIX AC021). The radar data were recorded from local time 12:00:00 midnight on 28 March 2012. Hundreds of commercial airline aircraft were tracked and recorded. Two cases are used here to investigate whether the performance of the IMMPDAF significantly varies in the area surrounding Taiwan. The first case target was from east head to Taiwan and the second case target flight away from Taiwan to north. The radar data was processed using the IMMPDAF and the NNKF. The ADS-B data have the target ‘Callsign’ and GPS calculated positions; these were set as the initial values of the two methods and regarded as the true positions. As described in the previous section IMMPDAF includes three filters of different dynamic models shown in Equations (18–20). The parameters of initial model probabilities and transition matrix are shown in Equation (21). The initial values of model probabilities have minor effect on the performance, and the IMMPDAF is very robust to the choice of transition matrix. Since the values set for P is 0.9 and P is 0.999, and the degree of freedom is 3, the gate threshold obtained from Chi-square table is approximately 16.

The performance requirements for APNT can be found in [17]. The minimum performance requirements for APNT are shown in Table 1. The FAA categorized the airspace into several zones [17]. En-route is the airspace at flight level (FL) 180 to FL 600 (18,000 to 60,000 feet). Terminal is the airspace starting at 500 feet above the ground and extending out to 5 statute miles from the airport; it then goes up at a 2-degree angle to FL 180. Under 500 feet, the APNT system needs to support the LNAV/non-precision approach. Accuracy (95%) is defined as the errors in the horizontal 95% under the accuracy requirement. That is, according to a normal distribution, 95% of error (error mean + 2 standard deviations) has to be under the accuracy requirement.

Table 1.

Required performance for APNT.

	Navigation	Surveillance
	Accuracy (95%)	Separation	NACp
En-route	10 nmi	5 nmi	0.1 nmi
	4 nmi
	2 nmi
Terminal	1 nmi	3 nmi	0.05 nmi
LNAV	0.3 nmi

Case 1

In Figure 2, the target heads to Taiwan Taoyuan International Airport (IATA: TPE, ICAO: RCTP) from east to west. The interval between epochs is 5 s. Note that a lot of surveillance radar systems have a scan time of close to 5 s. For example, ASR-12 surveillance radar has a scan time of 4 to 6 seconds. The tracking performance is shown in Figure 3. In case 1, the target starting at en-route airspace, so the accuracy requirement begins from 2 nmi. As shown in Figure 3, the NNKF experiences the instability at the beginning and the IMMPDAF performs more stable. The reason of the difference is the strategy to deal with the low quality measurements. While the NNKF just chooses to trust one of the measurements and updates, the IMMPDAF gives all the validated measurements weightings and updates with all of them. From Figure 3, the tracking performance of both the IMMPDAF and the NNKF meets the accuracy requirement. The IMMPDAF outperforms the NNKF, as shown in Table 2. The legends for Figures 3, 4 and 5 are illustrated below.

Figure 2.

Flight path for case 1 (target heading toward Taiwan).

Figure 3.

Performance of IMMPDAF and NNKF for Case 1.

Table 2.

Overall mean accuracies of the horizontal positioning for IMMPDAF and NNKF.

Accuracy (95%)	IMMPDAF	NNKF
Cases 1	0.311 nmi	0.921 nmi
cases 2	0.397 nmi	0.582 nmi
AVERAGE	0.354 nmi	0.751 nmi

Figure 4.

Flight path for Case 2 (target flying away from Taiwan).

Figure 5.

Performance of IMMPDAF and NNKF for Case 2.

Case 2

In Case 2, the target takes off from Taiwan Taoyuan International Airport and flies northeast as depicted in Figure 4. The tracking performance is shown in the upper plot of Figure 5, and the bottom plot of Figure 5 is the zoom-in plot of the upper one which shows the tracking performance for the LNAV/non-precision approach region. Note that the NNKF results initially exceed the accuracy requirement of APNT whereas the IMMPDAF results remain under the accuracy requirement. That is, the NNKF method does not meet the accuracy requirement (95%) in this case. As a result, the IMMPDAF outperforms the NNKF.

The mean accuracies of the horizontal positioning are summarized in Table 2. The IMMPDAF outperforms the NNKF for all tested cases. In comparison to the use of the NNKF, the average improvement on the positioning (tracking) performance gained from the proposed IMMPDAF is about 50%. The overall mean accuracies of the horizontal positioning are also listed in Table 2.

Computation Load

For all the experiments, the computations were performed in MATLAB (R2008a). A PC with an Intel Core 2 Duo E7500 2.93-GHz CPU and 2 GB of RAM was used. The computation performance was measured using a MATLAB function. The computation performance results are listed in Table 3. The values are total seconds that running a case needed. Each test runs 10 times and to obtain the average result. The average result is shown that both algorithms in the single target ATC radar tracking scenario are below the radar scan time an epoch (5 seconds).Without coding optimization, the IMMPDAF requires about nine times than that of the NNKF.

Table 3.

Computation performance of IMMPDAF and NNKF for Cases 1 to 2. Lower values are better.

Cost Time (s)	IMMPDAF	NNKF
CASE 1	19.754 s	2.249 s
CASE 2	25.110 s	2.734 s
AVERAGE	22.432 s	2.491 s

Conclusions

This paper presented a filter based on the Interacting Multiple Model (IMM) estimator and the Probabilistic Data Association (PDA) filter, and this filter is thus called IMMPDAF, which was applied to the radar tracking system. The real flight radar data was used to evaluate the tracking performance of the IMMPDAF, and the conventional Nearest Neighbor Kalman Filter (NNKF) was used as the baseline to show the performance gained from the proposed filter. Experimental results show that the tracking performance of the IMMPDAF is better than that of the NNKF. In one test case (i.e., Case 2), the IMMPDAF meets the accuracy requirement of the Alternate Positioning, Navigation, and Timing (APNT), but the NNKF results exceeded the APNT accuracy requirement for some periods of time. Although the computation load of the IMMPDAF is nine times that of the NNKF, the IMMPDAF is still suitable for ATC radar since processing time of each epoch is below the radar scan time (5 s). The proposed IMMPDAF gained about 50% improvement on tacking performance than that of the NNKF. Importantly, the radar tracking with the proposed IMMPDAF met the performance requirements of the APNT. That is, the current radar with the proposed IMMPDAF is capable to provide the navigation and surveillance services for the Air Traffic Management (ATM) system during periods of GPS outages.