Improved interacting multiple model algorithm airport surface target tracking based on geomagnetic sensors

To avoid the inherent defects of current airport surface surveillance systems, a distributed non-cooperative surface surveillance scheme based on geomagnetic sensor technology is proposed in this article. Furthermore, a surface target tracking algorithm based on improved interacting multiple model (WIMM) is presented for use when the target is perceptible. In this algorithm, the weighted sum of the mean values of the residual errors, which is used to reconstruct the model probabilistic likelihood function, and the Markov model transition probability are updated using posterior information. When a target is imperceptible, its trajectory can be predicted by the target identified motion model and the adaptive model transition probability. Simulation results show that the WIMM algorithm can be used efficiently together with an observed small sample of velocity information for target tracking and trajectory prediction. Compared with the interacting multiple model and residual-mean interacting multiple model algorithms, the frequency of model switching and the rate of model identification were increased during the imperceptible period, and target prediction error was greatly reduced.


Introduction
With the rapid development of China's civil aviation transportation industry and the rapid growth in passenger and cargo throughput, the ground of large airports is becoming more and more crowded, and conflicts between aircraft and ground vehicles are more likely to occur. In response to increasing safety and efficiency problems on the airdrome surface, EUROCONTROL proposed an Advanced Surface Movement Guidance and Control System (A-SMGCS). 1 The system performs surveillance, routing, guidance, and control based on various sensor technologies. At present, airdrome surface surveillance relies mainly on Surface Monitoring Radar (SMR), Automatic Dependent Surveillance-Broadcast (ADS-B), and Multilateration (MLAT). These three surveillance approaches, however, have the following inherent defects: (1) SMR is susceptible to factors like buildings, ground clutter, and weather, and (2) MLAT and ADS-B can only monitor a target equipped with a transponder on the surface. 2 However, surveillance based on distributed noncooperative perception can fundamentally avoid these 1 defects. Honeywell has developed a dual infrared/ magnetic sensor and arranged a set of these sensors to detect aircraft. 3 Chartier and Hashemi 4 proposed that the position of an aircraft could be determined through information from loop sensors located at segment boundaries. Gao et al. 5 proposed to detect a magnetic target using a magnetic sensor network and to transfer the resulting data to A-SMGCS. Schonefeld and Moller 6 conducted comprehensive analysis on the performance of a runway intrusion prevention system, XL-RIAS, based on distributed sensors, and testified that the response rate was faster than that of ASDE-X (Airport Surface Detection Equipment, Model X). Cheng and Wu 7 proposed to install microwave detection sensors on the side of the runway threshold, combined with an enhanced electronic process strip system to prevent runway incursions. Tang et al. proposed a runway incursion prevention method that the discrete airdrome surface traffic situation can be reconstructed based on the fusion of object sensing events detected by a sensor network, 8 and an asynchronous target-perceiving-event driven surface target surveillance scheme based on the geomagnetic sensor technology was proposed. 9 However, these studies have focused only on discrete positioning of surface targets with noncooperative sensors; continuous target tracking and prediction were not further considered.
Maneuvering target tracking with changing dynamics has been widely investigated in various applications such as aircraft traffic control (ATC). To resolve this mode-changing problem, a variety of methods have been investigated, such as the interacting multiple model (IMM) algorithm. 10,11 Two main trajectory tracking and prediction algorithms are investigated in this study. The first is an algorithm based on parameter identification of aircraft dynamics and kinematics models. Gong 12 studied the taxi speed and acceleration characteristics of the aircraft and obtained the fixedstructure interacting multiple model (FS-IMM) by regression analysis. The FS-IMM algorithm may not match the actual situation in the target tracking, which easily leads to target tracking failure. Capozzi et al. 13 analyzed historical surface surveillance data and identified the parameters of a kinematics equation model. Yao and Wang 14 proposed to estimate real-time target acceleration with a state observer, benefiting from fast convergence and high estimation precision without target motion model identification. The second algorithm is based on optimal estimation theory. Conventional Kalman filters like aÀb and aÀbÀg are single-model tracking and filtering algorithms, which are not suitable for the variety and uncertainty of a target motion model on the surface. 15 Farina et al. 16 applied the motion constraint on the airport surface to a variablestructure interacting multiple model (VS-IMM), which improved tracking precision. To address the deficient performance of the FS-IMM algorithm in target tracking, Gong et al. 17 applied the VS-IMM algorithm to surface target tracking in combination with an airport map. Hwang et al. 18 redefined the means of the residual errors of a Kalman filter and combined them with a probabilistic model likelihood function to improve flight path prediction accuracy. However, continuous observation of the target motion state is required for these algorithms, which runs by a large sample of radar measurements data set for tracking and prediction. 19 To avoid the inherent defects of current airport surface surveillance systems, a distributed non-cooperative surface surveillance scheme based on geomagnetic sensor technology is proposed in this article. However, due to the low density of the deployed nodes, the continuous motion state of a target at adjacent nodes could not be perceived. To address this problem, an weighted interacting multiple model (WIMM) algorithm is presented, in which a small sample of target velocity information is used to accelerate model switching and increase the identification rate of the motion model. When the motion state of the target is imperceptible, the target trajectory can be predicted using the motion model identified during the perceptible period as well as an adaptive model transition probability.
This article is organized as follows: the ''Geomagnetic sensing-based surface target surveillance scheme'' section proposes a geomagnetic sensing-based surface target surveillance scheme and target instantaneous velocity acquisition method; in the ''Target tracking and prediction based on the WIMM algorithm'' section, based on a general target motion model, target tracking and prediction based on the WIMM algorithm is proposed; in the ''Simulation and analysis'' section, Monte Carlo simulation is employed to compare the IMM with WIMM algorithms based on trajectory tracking and aircraft prediction performance in both perceptible and imperceptible states.
Geomagnetic sensing-based surface target surveillance scheme

Surface target surveillance scheme
In general, Moving targets in airport scenes include vehicles and aircraft, which are ferromagnetic objects that detect the position of non-cooperative targets based on magnetically induced anisotropic magnetoresistance geomagnetic sensors. 20 By using a wireless geomagnetic sensor network, non-cooperative motion targets can be effectively detected with high precision, small volume, low cost, no wiring, and deployment flexibility. Due to the large surveillance area, users of the geomagnetic technology-based surveillance scheme must consider the method of deployment and the number of geomagnetic detection nodes to reduce tracking cost and communication redundancy. It is known that with a given set of restrictions on a moving target, the target motion characteristics can be predicted in various airport areas. A new deployment mode for nodes with more space between them has been designed in combination with the characteristics of the moving target and prior information on the airport surface, as shown in Figure 1 (taking a taxiway section as an example).
The taxiing route of a moving target on the surface can be divided into n sections, L = fL 1 , L 2 , . . . , L i , . . . , L nÀ1 , L n g. Geomagnetic detection nodes, N = fN 1 , N 2 , . . . , N i , . . . , N nÀ2 , N nÀ1 g, are deployed on the runway or taxiway centerline of the segmentation boundary between adjacent sections. Each section is described by four parameters. For instance, section i can be defined as (L i , N i , N i + 1 , long i ), where L i is the section number, N i is the start node, N i + 1 is the terminal node, and long i denotes length. The section information is preserved in geomagnetic detection nodes for distributed computation after the nodes perceive a target. In this method of deployment, nodes can accurately perceive the target position as a target passes among them and can modify their previous position information. Velocity information can also be modified instantaneously by means of the target velocity obtained from the nodes. Of course, the use of geomagnetic sensors also has limitations. First, the moving target can only be a substance that has an influence on the geomagnetic sensor. Second, the target should not deviate too far from the runway and the taxiway.

Target instantaneous velocity acquisition method
In the system proposed here, a geomagnetic detection node is equipped with dual sensors to obtain the instantaneous speed of the target. Let the distance between the two geomagnetic sensors be D. When a target passes the geomagnetic detection node, the two geomagnetic sensors have asynchronous access to the geomagnetic induction signal, the offset of which depends on the target's real-time velocity. The geomagnetic induction signal of a certain surface target obtained from a geomagnetic detection node is assumed to be as shown in Figure 2.
As can be seen in Figure 2, induction signals 1 and 2 of the two geomagnetic sensors on the same detection node share a few common feature points, including the detection starting point onset, peak, and trough, the extreme point, and the detection ending point. Assuming that the two sensors have the same parameter index, the corresponding feature points of the two induction signals can be regarded as the geomagnetic induction signal acquired when one target cross section passes the two geomagnetic sensors. The instantaneous velocity sequence of several sections can be obtained from the sampling offset of corresponding feature points, and the velocity sequence can be used as an instantaneous velocity sequence when a target passes a geomagnetic detection node.
Assuming that ds i is the offset of the ith corresponding feature point of signals 1 and 2 in the sampling sequence and that the sampling frequency of the two  sensors is f s , the estimation formula for the instantaneous velocity v(t i ) can be given as The , v(t n )g can be obtained through a single geomagnetic detection node when a target passes it.
Target tracking and prediction must be conducted on the time sequence using the same time interval, but because the interval between key points on a given geomagnetic induction signal curve is variable, it is necessary to transform the instantaneous velocity sequence into an equal-interval sequence. This research used the quadratic Lagrange interpolation method to perform interpolation on the acquired instantaneous velocity sequences so that their time intervals are equal. Given that q 1 = t 1 , q n = t n , the equal interval is defined as Dq, and the formula can be written as The interpolation formula can be expressed as where t jÀ1 ł q i ł t j + 1 , i = 1, 2, . . . , n À 1, n ; j = 2, 3, . . . , n À 2, n À 1.
After interpolation, the equal-interval velocity sequence can be written as

Dynamic model
The target's motion model is generally expressed as a dynamic model and a measurement model. 21 The dynamic and measurement models can be expressed as where k À 1, k represent sampling times, X is the state vector, F is the state transition matrix, Z is the observed values, H is the measurement matrix, w is motion process noise that obeys a Gaussian distribution with mean 0 and variance Q, and v is observation noise that obeys a Gaussian distribution with mean 0 and variance R.
Moving targets in airport scenes include vehicles and aircraft, which are ferromagnetic objects that detect the position of non-cooperative targets based on magnetically induced anisotropic magnetoresistance geomagnetic sensors. State transition matrices for these three models can respectively be expressed as follows In the constant turn (CT) model, the state variable X = ½x, _ x, y, _ y 0 is chosen, assuming that € x and € y are treated as random noise, that is, € where t is the sampling interval. The motion state vector of the aircraft is X = ½x, _ x, € x, x

WIMM algorithm
The WIMM algorithm improves the IMM algorithm in two ways. First, a weighted sum is performed on the means of the residual errors, and the model probabilistic likelihood function is reconstructed, thus increasing the probability of identifying the true motion model. Second, the model transition probability is updated for self-adaptation using the model posterior probability, thus accelerating model switching as well as increasing the model identification rate. Figure 3 shows a schematic diagram of the WIMM algorithm. This algorithm consists of the following five steps: input interaction, Kalman filtering, model probability update, model transition probability self-adaptation, and output fusion.
Input interaction. Assuming that a model set consists of r motion models, the state estimation value and the covariance matrix of each model at time k À 1 are respectively as follows:X j (k À 1jk À 1) andP j (k À 1jk À 1), j = 1, 2, . . . , r. u ij is the probability of converting the i model into the j model. After interaction, the input in model j at time k can be expressed as follows Among them, the filtering of the IMM algorithm is realized by Kalman filtering and will not be described here.
Model probability update. In the IMM algorithm, the maximum likelihood function at time k for model j is as given in equation (8) The residual sequence obeys a Gaussian distribution with mean 0 and variance S j (k), and the motion model set can contain all motion models during the operation. However, due to factors like target movement uncertainty, surface restrictions, and maneuverings, the target motion model may exceed the model set in the algorithm. In addition, motion models cannot be added endlessly to the conventional IMM algorithm for convenience in computation. Multiple models will also tend to compete with each other, thus reducing trajectory tracking accuracy. Therefore, if a motion model is not added, the model probabilistic likelihood function is reconstructed and is no longer assumed to obey a Gaussian distribution with mean 0 and variance S j (k).
Let us assume the true motion model of a moving target on the surface to be as follows The model state transition matrix error is as follows: DF j = F T À F j . Then the error of target state estimation 15 can be defined as follows An expression for the state estimation error after input interaction 16 can be obtained as From the true motion model, the state estimation error after input interaction, the expression in equation (13) can be obtained e j kjk ð Þ= I À K j k ð ÞH À Á F T e 0j k À 1jk À 1 ð ÞÀk j k ð Þv T k ð Þ The residual error of model j is obtained as The residual error r j (k) can be expressed using the state estimation error after input interaction as The mean value obtained from equations (13) and (15) can be expressed as follows Because of the uncertainty of the true motion model of a surface target, quantization of F T cannot be performed, and therefore the mean of the residual errors cannot be obtained. In this article, each model in the motion model set is regarded as the real motion model of the moving object, and a weighted sum is calculated using the posterior probability of each model. The expression in equation (18) can be obtained by summing the means of the residual errors where e 0j (k À 1jk À 1)= P r i= 1 u ij (k À 1jk À 1) e i (k À 1). Then the maximum likelihood function for model j can be obtained as The posterior probability is Markov model transition probability self-adaptation. In the IMM algorithm, the input interaction process is generally a Markov process, and the initial transition probability of the model is determined by humans and does not change with actual motion conditions. When the deviation between the actual motion of the target and the initial transition probability is large, the artificially determined Markov matrix cannot accurately reflect the actual motion mode transition of the target due to the distortion of the situation. In view of this, the article uses the runway and taxiway in the airport. A priori information such as the probability of running is corrected in real time to make it fit the actual motion state of the target. Therefore, the Markov model transition probability p ij is updated using the posterior information in the WIMM algorithm to solve this problem.
Assuming that the probability in model j at time k À 1 is u j (k À 1) and that at time k is u j (k), the probability differential value of the same model at adjacent times reflects the change in the degree of matching between model j and the true motion model. The rate of change of the posterior probability in model j can be defined as Let the probability that the target motion model i transitions from model j at time k À 1 be p ij (k À 1), and update the Markov model transition probability with Du j (k). Then the expression in equation (21) is obtained In equation (21), a is an adaptive parameter, l is a control parameter for adaptive velocity control, and both values relate to the practical situation of the moving target on the surface. Whether the Markov transition probability model is fixed or self-adaptive, it must satisfy the two basic properties given in equation (22) 0\p ij \1, i, j = 1, 2, . . . , r P r j = 1 Then normalization must be performed on p 0 ij (k), and then the updated transition probability p ij (k) can be obtained as As can be seen from equation (23), when one motion model transitions to model j, the updated p ij (k), i = 1, 2, . . . , r increases as the posterior information u j (k) increases, meaning that model j plays a critical role in input interaction at the next time period.
Output fusion. The interactive output results at time k can be expressed as followŝ

Trajectory prediction in an imperceptible state
When a target is moving in the section between two adjacent nodes, since the distance between the two nodes exceeds the length of the aircraft, meaning that velocity information cannot be obtained; therefore, extrapolated prediction is needed. At the last moment when the target is perceptible, the WIMM algorithm provides identification of each model in the model set, namely the model posterior probability u j (k), j = 1, 2, . . . , r. Then, given the self-adapting Markov model transition probability p ij (k), the expression for predicting the probability of each model in the model set when a target is not perceptible at time k + 1 can be obtained After substituting the predicted state valuê X j (k + 1jk) and the predicted model probability u j (k + 1jk) of each model into equations (25) and (26), the predicted state value and the prediction filtering covariance matrix when the target is imperceptible can be defined aŝ Simulation and analysis

Simulation condition and parameters
An aircraft passing a certain geomagnetic detection node on the taxiway is taken as an example in this research and compared with the IMM, residual-mean interacting multiple model (RMIMM), and WIMM algorithms with regard to trajectory tracking and aircraft prediction performance in both perceptible and imperceptible states using Monte Carlo simulation. Assume that the process of aircraft operation is as follows. First, when the aircraft is in the perceptible period: (1) 0-4.5 s for the CT model; (2) constant acceleration (CA) is performed at 0.45 m/s 2 from 4.8 to 12 s;(3) constant velocity (CV) is performed at the velocity obtained from the first two steps from 12.3 to 15 s. Second, when the aircraft is in the imperceptible period: (4) the aircraft maintains CV for about 30 s at the velocity obtained from step 2 and then operates to the next detection node. Figure 4 shows the actual position of the aircraft according to the operating procedure.
The simulation parameters are selected as follows: the noise covariance in each model during the estimation process is Q =  Simulation results and analysis during the perceptible period Figures 5 and 6 show simulation results for aircraft position and velocity during the period when the aircraft is perceptible. Figures 5 and 6 illustrate that the IMM, RMIMM, and WIMM algorithms are all suitable for tracking of aircraft motion state during the period when the aircraft is perceptible. However, the WIMM algorithm is superior to the other two algorithms for tracking. To show the advantage of WIMM more clearly, the position error curve and the velocity root mean square error (RMSE) curve for the WIMM and IMM algorithms are respectively plotted in Figures 7 and 8. Figure 7 shows that during the period when the aircraft is perceptible, the maximum position error using the IMM algorithm is 2.13 m and that using the RMIMM algorithm is 1.30 m, but when using the proposed algorithm, the maximum position error is only 0.42 m. Figure 8 shows that the maximum RMSE error using the IMM algorithm is 0.059 m/s and that using the RMIMM algorithm is 0.356 m/s, but when using the WIMM algorithm, the maximum RMSE error is  only 0.023 m/s. As can be seen from these results, tracking precision is greatly improved when using the WIMM algorithm. Figure 9(a) shows the IMM algorithm, Figure 9(b) shows the RMIMM algorithm, and Figure 9(c) shows the WIMM algorithm and gives the selection probability curves of the CV, CA, and CT models. Figure 9 demonstrates that the IMM algorithm cannot identify each motion model clearly and that the intersection points of the three selection probability curves are close together. Moreover, switching among motion models is slowed down, and flexibility is also reduced. The RMIMM algorithm is more advantageous than the IMM algorithm, yet still needs to be improved. The WIMM algorithm, however, can solve the problems encountered by the first two algorithms, not only facilitating estimation of aircraft motion state using observed data but also effectively increasing the degree of model identification. Table 1 shows a comparison of the highest identification rate of each motion model and the model switching time among the three algorithms.

Simulation results and analysis during the imperceptible period
Simulation results of trajectory prediction during the period when the aircraft is not perceptible are displayed in Figures 10 and 11.   Figure 11. Position error curve. Figure 10 illustrates that as aircraft operation time increases, the error between the positions predicted using the IMM, RMIMM, and WIMM algorithms and the actual position increases. However, the predicted position is closest to the actual position when using the WIMM algorithm. Figure 11 illustrates that at the last moment of position prediction, the prediction error has accumulated to 10.32 m when the IMM algorithm is used and 5.72 m when RMIMM is used, but only 2.16 m with the WIMM algorithm. It is apparent that the WIMM algorithm outperforms the IMM and RMIMM algorithms in terms of extrapolated trajectory prediction, particularly during the period when the aircraft is not perceptible.

Conclusion
This article proposes a distributed non-cooperative scene target monitoring scheme based on geomagnetic sensing technology. In addition, based on the WIMM algorithm, the trajectory tracking and prediction algorithm for moving objects in airport scenes is derived. The algorithm improves the IMM algorithm by calculating the weighted sum based on the mean of the residuals and reconstructing the model probability likelihood function. The model transition probability is updated for self-adaptation with posterior information, thus accelerating model switching as well as increasing the model identification rate. Finally, this article presents simulation results for target tracking and prediction both when a target is perceptible and when it is not perceptible using two algorithms. From simulation results, both from periods when the aircraft is perceptible and when it is imperceptible, it is apparent that the WIMM algorithm is superior to the conventional IMM and RMIMM algorithms in terms of trajectory tracking and prediction, especially when the aircraft is imperceptible.

Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.