Hybrid 4-Dimensional Trajectory Prediction Model, Based on the Reconstruction of Prediction Time Span for Aircraft en Route

Abstract

This paper presents the results from a test of the performance of several general trajectory prediction methods and proposes a hybrid trajectory prediction model that aims to increase the safety of flights en route and improve airspace management capabilities by predicting the aircraft’s four-dimensional trajectory (4DT) more accurately. The automatic dependent surveillance-broadcast (ADS-B) data from 589 trajectories of cruising aircraft from the Guangzhou area were extracted for experiments. Numerous trajectory prediction methods, including velocity trend extrapolation, long short-term memory (LSTM), stateful-LSTM, back propagation (BP) neural network, a one-dimensional convolutional neural network (1D-ConvNet), Kalman filter, and flight plan interpolation were used for prediction experiments, and their performance at different time spans of prediction is obtained. By extracting the best methods using different time spans of prediction, a hybrid prediction model is proposed based on the reconstruction of these methods. For the data in this paper, the mean squared error (MSE) of the hybrid prediction model is significantly reduced compared to other methods in different time spans of prediction, which has great significance for future trajectory prediction in a structured airspace.

1. Introduction

A series of air traffic management problems, such as flight delays, airspace congestion, controller overload, and airspace flow imbalance, have arisen as a result of the rapid development of civil aviation around the world. In order to promote the balance of airspace flow capacity, reduce the workload of controllers, and improve the safety and efficiency of air traffic management (ATM), Eurocontrol proposed a project called the “Single European Sky ATM Research” (SESAR) in 2004. One of the most important targets of SESAR is to realize the transformation from an airspace-based operation to embrace the concept of trajectory-based operation (TBO), further enhancing the refinement, coordination, and predictability of ATM. In 2005, the United States proposed their next-generation transportation system (NextGen), which treats TBO as one of the basic features necessary for an advanced air traffic management system. In 2016, the civil aviation administration of China (CAAC) adopted TBO as one of the critical concepts in civil aviation. There is an increasing number of aviation organizations in the world that regard TBO as a new operating concept and attach great importance to it. In 2021, the German air traffic control and navigation service company, Deutsche Flugsicherung (DFS), and the China Electronics Technology Group Corporation (CETC) started joint research on the TBO and the Free Route Airspace (FRA) control system, which aims to pave the way for the implementation of TBO in China. Trajectory prediction is one of the most important technologies of TBO because the accuracy of trajectory prediction will directly affect the performance of trajectory planning, cross-sector air traffic control, conflict detection, and many other aspects.

The methods used for trajectory prediction in the current research can be mainly divided into three categories [1,2]:

• State estimation model: This type of trajectory prediction method usually makes predictions based on the aircraft dynamics model, motion characteristics, and the Markov property of the aircraft’s state. This type of prediction method has strong solvability and offers excellent short-term prediction performance due to large aircraft usually having a small rate of change in terms of motion parameters. The state of the aircraft at the next moment is only approximately related to the aircraft’s state at the current moment, based on the Markov theory [3]; therefore, the aircraft’s trajectory can be predicted according to the current state of the aircraft, when combined with its cybernetics model. At the same time, such methods are often mixed with theories on aircraft performance, dynamics, and interacting multiple models (IMM) [4]. Extrapolation of the speed trend is used in the short-term trajectory prediction of air traffic control (usually 6s on the radar screen) at present, which technique can be classified as a simple state estimation model [5]. However, the prediction error of these state estimation models usually increases rapidly as the time span of the prediction gets longer.
• Machine learning model: The machine learning model used for trajectory prediction generally falls into the category of supervised learning. By learning from a large amount of historical trajectory data, the universal characteristics of machine learning models could be extracted by the model. These methods can also use aircraft performance, air-route structure, and environmental parameters as extra inputs to predict the trajectory (data to data). Machine learning models always achieve good prediction performance and benefit from extensive trainable parameters [6]. These kinds of methods generally extract the characteristics of data changes via data analysis to predict the trajectory, using filter-based methods [7,8] or deep learning-based methods [9]. In recent years, the long short-term memory (LSTM) model, which is improved by a recurrent neural network, has been shown to offer good performance in processing time-series data [10,11]. Some trajectory prediction models have combined various neural networks, such as convolutional neural networks(CNN) and LSTM, to predict trajectory data and have achieved a good performance [12]. However, the machine learning model usually only learns a specific type of trajectory; if there is a need to predict the trajectory of different aircraft types on different air routes, researchers must retrain the model, adding or changing the parameters of the machine learning model, and the cost of retraining in terms of time can sometimes be very high.
• Flight plan-based model: Flight planning is a pre-tactical trajectory planning method; the preparation of flight plans involves coordination between individual flights, air traffic control units, and airports [13]. Therefore, pilots and controllers will try to make sure the required time of arrival (RTA) of the aircraft is consistent with the flight plan [14]. Although the contribution of machine learning to trajectory prediction is great, we should not ignore the important influence of flight-plan information on aircraft trajectory. It is also a general medium- and long-term trajectory prediction method that is widely used by dispatchers and controllers at present [15]. This method is susceptible to external influences from exceptional circumstances (such as diversions, alternate, and delays caused by air traffic flow management), but the prediction performance based on the flight plan is more stable compared to other methods analyzing the whole duration of flights [16].

The above-mentioned various trajectory prediction methods have different prediction performances with different time spans of prediction. This paper aims to use the same data set and the same prediction performance evaluation standard to explore the performance of different trajectory prediction methods, thus constructing a hybrid prediction model that works under different time spans of predictions. Firstly, the methods that are used include velocity trend extrapolation, LSTM, stateful-LSTM, BP neural network, one-dimensional convolutional neural network (1D-ConvNet), Kalman filter, and flight plan interpolation to predict a trajectory and then evaluate their performance. Secondly, based on the prediction performance of those methods, a hybrid four-dimensional trajectory prediction model is proposed that switches from one model to another strategically, thereby taking advantage of the benefits offered by different models. Finally, by comparing the MSE of different prediction methods with a flight plan interpolation, the longest time span over which those methods can predict trajectory accurately can be determined.

This paper describes the construction idea, methods, and performance evaluation indicators of the hybrid 4D trajectory prediction model shown in Section 2. We use the ADS-B 4D trajectory data to conduct experiments on different trajectory prediction models and explore their prediction performance in Section 3. Then, a hybrid trajectory prediction model is constructed, and its effectiveness is verified in Section 4. Finally, the conclusions of this research are in Section 5.

2. Hybrid 4-D Trajectory Prediction Model

The process of constructing a hybrid 4-D trajectory prediction model can be described as follows: (1) firstly, it is necessary to process the trajectory data to make it adaptable to general trajectory prediction methods (Section 2.1). (2) To ensure the comparability of the prediction performance of these different methods, it is necessary to determine a general trajectory prediction model and a performance indicator of different trajectory prediction models (Section 2.2 and Section 2.3). (3) Then, a candidate pool that includes multiple types of trajectory prediction methods should be built for the subsequent process of reconstruction (Section 2.4). (4) Finally, the hybrid trajectory prediction model should be constructed, as described in Section 2.5.

2.1. Structure and Preprocess of 4-D Trajectory Data

2.1.1. Basic Structure of the Trajectory Data

In order to make the trajectory prediction more accurate, we extracted environmental data from the European Center for Medium-range Weather Forecasts (ECMWF), which includes the direction and speed of the wind in the area of the trajectory data at each longitude, latitude, and height (GRIB file).

The typical trajectory intention of an aircraft is generally designed along the center of the route; the waypoints are entered into the flight management computer (FMC) before the aircraft takes off. When an aircraft is cruising, it often deviates from the centerline of the air route due to the reduced positioning performance and high-altitude wind, which is a common phenomenon of trajectory data. We define the trajectory intention deviation $\sigma$ as the extra input for trajectory prediction. Mathematically, the trajectory intention deviation $\sigma$ is:

$$\sigma =\min \left( \frac{\left|A_{i}x+B_{i}y+C_i\right|}{\sqrt{A_i{}^2+B_i{}^2}}\right)$$

(1)

In Equation (1), $A_i, B_i$ and $C_i$ are the linear parameters of the $i$-th trajectory’s nominal air route, $x$ and $y$ are the longitude and latitude of the trajectory point; $\sigma$ is used to measure the deviation of the aircraft’s position from the nominal air route.

Therefore, the structure of the input data for prediction is:

$${\mathit{X}}_{it}={\left[x,y,z,v_{\mathrm{aircraft}},v_{\mathrm{diraircraft}},v_{\mathrm{wind}},v_{\mathrm{dirwind}},\sigma \right]}_{it}$$

(2)

In Equation (2), for the t-second trajectory data of the i-th trajectory, $x$ is the longitude of the trajectory point, $y$ is the latitude of the trajectory point, $z$ is the query of normal elevation (QNE) of the aircraft, $v_{\mathrm{aircraft}}$ is the speed of the aircraft (m/s), $v_{\mathrm{diraircraft}}$ is the heading of the aircraft (degree), $v_{\mathrm{wind}}$ is the wind speed of the current trajectory point (m/s), $v_{\mathrm{dirwind}}$ is the wind direction (degree) of the current trajectory point, and $\sigma$ is the trajectory intention deviation (degree) between the trajectory point and the standard route.

2.1.2. Data Preprocess

(1) Interpolation: The current ADS-B trajectory data mainly have update frequencies of 20 s, 10 s, 5 s, and 1 s (the interval of the ADS-B data in this paper is 10 s). In order to improve the universality of the trajectory data, the trajectory data need to be processed into intervals of 1 s. In order to get a more accurate trajectory prediction performance from the different methods, a smooth interpolation based on the Bezier curve is performed on each attribute of the ADS-B data, which is shown in Equation (3):

$$p\left(t\right)=p_0{(1-t)}^3+3p_{1}t{(1-t)}^2+3p_2{t}^2\left(1-t\right)+p_3{t}^3,t\in \left[0,1,\mathrm{step}=0.1\right]$$

(3)

In Equation (3), t is the time step that needs to be interpolated, and $p_0,p_1,p_2,p_3$ are four adjacent trajectory points. $p\left(t\right)$ is the interpolated trajectory point after smooth interpolation.

(2) Normalization: After interpolation, the ADS-B data is updated every 1 s. We randomly select 70% of the total trajectory data as the training data (I = 412 in this paper, see Equation (17)), and the remaining 30% as the validation data (I’ = 117 in this paper, see Equation (11)). All training data are normalized according to their means and variances. Mathematically, the process is:

$$\widehat{\xi }=\frac{\xi -u_{\xi }}{{\delta }_{\xi }+\epsilon }\xi \in X$$

(4)

Similarly, the predicted trajectory data needs to be de-normalized. Mathematically, the process is:

$$\xi ={\delta }_{\xi }\cdot \widehat{\xi }+u_{\xi }$$

(5)

In Equations (4) and (5), $\xi$ is one of the trajectory attributes, $u_{\xi }$ is the mean value of the trajectory attributes, $\xi$, in the training data set, ${\delta }_{\xi }$ is the variance of the trajectory attributes $\xi$ in the training data set, and $\epsilon$ is a small positive number ($1\times {10}^{-14}$ in this paper).

(3) Timestamp: All of the trajectory samples regard the first trajectory point as the first trajectory point where the timestamp t = 1, and the interval after the interpolate process is set at 1 s. Thus, the second trajectory point has the timestamp of t = 2 and so on.

In order to meet the training and fitting needs of each prediction method, trajectories from the training data set are processed into the following two-dimensional form, as shown in Figure 1.

To prevent errors caused by data splicing, the time span of the prediction will end at the shortest trajectory running time ($\mathrm{T}$ = 2046 s in this paper, see Equation (11)). At the same time, the input data and output data should be avoided, to cross the splicing boundary while constructing training data for the machine learning model.

2.2. Trajectory Prediction Mode

Taking the data and methods of this paper as an example, trajectory prediction methods use the first four trajectory points as the starting data input to predict the next point of the trajectory data. All these predictions are based on single-step, bootstrap, and global modes: that is, each prediction predicts the trajectory data of the next second (single-step), the predicted result of the trajectory data is used as the subsequent input parameter (bootstrap), and the process is repeated until it is possible to predict the range of the prediction’s time span (global). The flight plan interpolation does not use the prediction data structure in Figure 1 but instead uses the flight plan to predict the subsequent trajectory points.

For all of the subsequent machine learning prediction methods in this paper, the next input parameter needs to be kept in the form of 1 × 8 but the output prediction result is in the form of 1 × 3, in the process of continuous bootstrap prediction. Therefore, the output result should append the aircraft speed, direction, wind direction, wind speed, and trajectory intention deviation parameters, which are set in advance according to the position of the aircraft. As a result, the trajectory prediction algorithm, based on the single-step, bootstrap, and global modes, can be described as in Algorithm 1.

Algorithm 1 Trajectory prediction based on single-step, bootstrap, and global modes
Input: The first four trajectory points near ODOPT in the validation data set:$\left\{X_{it}\right\} i\in \left[0,\mathrm{I}’\right],\mathrm{t}\in \left[1,4\right]$. The real trajectory points are in the validation data set {${\widehat{p}}_{it}$},$i\in \left[0,\mathrm{I}’\right]$,$t\in \left[1,\mathrm{T}\right]$. The trajectory prediction model:$\mathrm{model}$.
Output: The predicted trajectory points $\left\{p_{it}\right\} i\in \left[0,\mathrm{I}’\right],\mathrm{t}\in \left[5,\mathrm{T}\right]$.
Hyper-parameters: $\mathrm{T}=2046$,$I^{\prime }$ = 177 for the data in this paper, wind data:$\mathrm{data} \mathrm{frame}$
Initialize: $i=1,t=1$. While$i$< $I^{\prime }$: #($I^{\prime }$ = 177 in this paper, the trajectories in the validation data)
While $t$ < $\mathrm{T}$: #($\mathrm{T}$ = 2046 in this paper)
$p_{it+1}=\left(x_{t+1},y_{t+1},z_{t+1}\right)=\mathrm{model}\left(\left\{X_{it}\right\}\right)$	(6)
$v_{\mathrm{aircraft}}{}^{t+1}={\left|\left|p_{t+1}-p_t\right|\right|}_2=\sqrt{{\left(x_{t+1}-x_t\right)}^2+{\left(y_{t+1}-y_t\right)}^2+{\left(z_{t+1}-z_t\right)}^2}$	(7a)
$v_{\mathrm{diraircraft}}{}^{t+1}=\frac{\left|\left(x_{t+1}-x_t\right)+\left(y_{t+1}-y_t\right)+\left(z_{t+1}-z_t\right)\right|}{{\left(x_{t+1}-x_t,y_{t+1}-y_t,z_{t+1}-z_t\right)}_2}$	(7b)
$\left(v_{\mathrm{wind}}{}^{t+1},v_{\mathrm{dirwind}}{}^{t+1}\right)=\mathrm{data} \mathrm{frame}\left(x_{t+1},y_{t+1},z_{t+1}\right)$	(8)
$\sigma =\min \left( \frac{\left|A_{i}x+B_{i}y+C_i\right|}{\sqrt{A_i{}^2+B_i{}^2}}\right)$	(9)
$X_{t+1}=\left[x_{t+1},y_{t+1},z_{t+1},v_{\mathrm{aircraft}}{}^{t+1},v_{\mathrm{diraircraf}}{}^{t+1},v_{\mathrm{wind}}{}^{t+1},v_{\mathrm{dirwind}}{}^{t+1},{\sigma }^{t+1}\right]$	(10)
$t$ = $t$+ 1 $i=i+1$

In Equation (6),$p_{t+1}$ is the aircraft’s position data at the next time step predicted by the $\mathrm{model}\left(.\right)$, which is a trajectory predict function, based on the prediction model; $p_t$ is the aircraft position data at the current time step. $X_t$ is the input data for trajectory prediction; it has the same structure as Equation (2). Equation (7a,b) indicates that the speed of the aircraft is calculated by the predicted position of the aircraft. The wind information used in $X_t$ is read from the data frame of the ECMWF (Figure 2) according to the position and the time of the aircraft in Equation (8). $\sigma$ is the trajectory intention deviation that has been explained in Equation (1). As a result, the input data for the next step of trajectory prediction is shown as Equation (10).

2.3. Indicator of Trajectory Prediction Performance (Objective Function)

Each trajectory prediction method generates the optimal trajectory prediction model from the training data set (for the deep learning model, this is achieved by exploring the best structure and using early stopping to establish the optimal model). Then, we use the optimal prediction model to predict the trajectory data in the validation data set, to evaluate the prediction performance. All trajectory prediction methods use the MSE in the validation data set as the evaluation standard for prediction performance. Mathematically, MSE is:

$$\mathrm{MSE}\left(\mathrm{degree}\right)=\frac{1}{I^{\prime }·T}·{{\displaystyle \sum }}_{i=1}^{I^{\prime }}{{\displaystyle \sum }}_{t=1}^T\parallel {\widehat{p}}_{it}-p_{it}{}_2\parallel =\frac{1}{I^{\prime }·T}·{{\displaystyle \sum }}_{i=1}^{I^{\prime }}{{\displaystyle \sum }}_{t=1}^T{\left({\widehat{p}}_{it}-p_{it}\right)}^2$$

(11)

In Equation (11), $I^{\prime }$ is the number of the trajectories ($I^{\prime }=$177 in this paper). T is the range of the prediction time span ($\mathrm{T}=$ 2046 in this paper).$p_{it}$ is a vector that indicates the predicted position of the t-th point of the $i$-th trajectory.${\widehat{p}}_{it}$ is a vector that indicates the actual position of the $t$-th point of the $i$-th trajectory.

There are several reasons why the MSE is used to evaluate prediction performance: on the one hand, MSE is an indicator for measuring prediction accuracy. On the other hand, with the latitude and longitude distance conversion formula, it can reflect the Euclidean distance error between the predicted position of the trajectory point and the actual position of the trajectory point directly, as shown in Equations (12)–(14):

$$p_{it}=\left(x_{it},y_{it},z_{it}\right), {\widehat{p}}_{it}=\left({\widehat{x}}_{it},{\widehat{y}}_{it},{\widehat{z}}_{it}\right)$$

(12)

$$\mathsf{\omega }=\arccos \left(\sin \left(y_{it}\right)·\sin \left({\widehat{y}}_{it}\right)·\cos \left(x_{it}-{\widehat{x}}_{it}\right)+\cos \left(y_{it}\right)·\cos \left({\widehat{y}}_{it}\right)\right)$$

(13)

$$\mathrm{MSE}\left(\mathrm{km}\right)=\sqrt{{\left(\frac{\pi \mathrm{R}\omega }{180}\right)}^2+{\left(z_{it}-{\widehat{z}}_{it}\right)}^2}$$

(14)

In Equations (12)–(14), $p_{it}$ is a vector that indicates the predicted position of the t-th point of the $i$-th trajectory; $p_{it}$ consists of $x_{it},y_{it},z_{it}$, where $x_{it}$ is the predicted longitude, $y_{it}$ is the predicted latitude, and $z_{it}$ is the predicted QNE. ${\widehat{p}}_{it}$ is a vector that indicates the actual position of the $t$-th point of the $i$-th trajectory, ${\widehat{p}}_{it}$ consists of ${\widehat{x}}_{it},{\widehat{y}}_{it},{\widehat{z}}_{it}$, ${\widehat{x}}_{it}$ is the actual longitude, ${\widehat{y}}_{it}$ is the actual latitude, ${\widehat{z}}_{it}$ is the actual QNE. $\mathrm{R}$ is the radius of the earth (6371 km). $\mathrm{MSE} \left(\mathrm{km}\right)$ is the MSE expressed in kilometers, which is used as the indicator for the performance of trajectory prediction models.

2.4. Candidate Pool of Prediction Models

The core idea of constructing the hybrid trajectory prediction model is to select the combination of trajectory prediction methods with the best prediction performance for different time spans of prediction. Then, we reconstruct a hybrid trajectory prediction model, based on the best method for the different time spans of prediction.

After determining the trajectory prediction mode and the indicator of trajectory prediction performance, it is necessary to select a series of trajectory prediction methods for trajectory prediction experiments. These trajectory prediction models constitute the candidate pool; the subsequent hybrid trajectory prediction methods are composed of models in this candidate pool. The reason why multiple trajectory prediction models are placed in the candidate pool is that for an unknown trajectory dataset, we are not sure which model will be better and which will be worse. Generally, the candidate pool should contain different types of trajectory prediction models because the models of different types may show a significantly different performance. Therefore, we can select the sub-model with better performance in a certain time span of prediction from among them. The resulting hybrid trajectory prediction model is reconstructed based on the best-performing model in the pool; the better the prediction performance of the models in the candidate pool, the better the trajectory prediction performance of the hybrid trajectory prediction model will be.

Seven trajectory prediction methods from three categories mentioned in Section 1 are chosen to construct the candidate pool, then these methods are used to predict the trajectory and calculate its performance using prediction experiments. The methods in the candidate pool are shown in

Table 1. Trajectory prediction methods in the candidate pool.

Type of Method	Name of Method
State estimation model	Velocity trend extrapolation (Section 3.2.1)
Machine learning model	BP neural network (Section 3.2.2)
	LSTM (Section 3.2.3)
	Stateful-LSTM (Section 3.2.4)
	1D-ConvNet (Section 3.2.5)
	Kalman filter (Section 3.2.6)
Flight plan-based model	Flight plan interpolation (Section 3.2.7)

.

Several trajectory prediction methods need to be placed in the candidate pool because we do not know how well these methods perform in this dataset before the prediction experiment starts. However, too many methods in the candidate pool will also make the construction of a hybrid trajectory prediction model more complicated and increase the workload. Therefore, choosing an appropriate number of the different types of trajectory prediction methods will help construct a hybrid trajectory prediction model.

2.5. Hybrid Model Reconstruction

Before building a hybrid trajectory prediction model, all the trajectory prediction methods in the candidate pool must be tested to explore their prediction performance, then select one of the best trajectory prediction methods at each time span of prediction and recombine to form a hybrid trajectory prediction model, according to the performance of different prediction methods on different time spans of prediction. The process of reconstructing the hybrid trajectory prediction model is:

$$\mathrm{Hybrid}\_\mathrm{Model}=\left\{{\mathrm{Model}}_1{}^{t1},{ \mathrm{Model}}_2{}^{t2},\dots ,{ \mathrm{Model}}_k{}^{ti}\right\}$$

(15)

$${\mathrm{Model}}_k{}^{ti}=\underset{\mathrm{Model}}{\mathrm{argmin}}\left\{MS{E}_k{}^{ti}\right\}, Mode{l}_k\in \mathrm{candidate} \mathrm{pool},t_i\in \mathrm{T}$$

(16)

Equations (15) and (16) mean that the hybrid prediction model is constructed by the different models with the best performance over different prediction time spans. In Equations (15) and (16), ${\mathrm{Model}}_k{}^{ti}$ indicates the k-th model of the prediction time spans $t_i$, and $MS{E}_k{}^{ti}$ indicates the mean square error (MSE) of the k-th model of the prediction time spans $t_i$. All of the models are in the candidate pool and the time span of prediction is in the range of $\mathrm{T}$. The schematic diagram of the construction of the hybrid trajectory prediction model is shown in Figure 3.

Due to using the same data processing method, the same evaluation criteria of the trajectory performance, and the same prediction mode selected by the different methods, the prediction accuracy of each method becomes comparable. The hybrid trajectory prediction model is constructed by the sub-model with the best performance in the different time spans of prediction, so the hybrid model will at least perform no worse than any of the individual models that comprise it.

3. Trajectory Prediction Experiments

In Section 3.1, 589 ADS-B trajectory data points generated by aircraft cruising were extracted as the original trajectory data set, which were divided into a training data set (412 of them) and a validation set (177 of them) after the data processing method described in Section 2.1. In Section 3.2, 7 trajectory prediction methods are used for experiments and their prediction performances at different prediction time spans are obtained. Part of the data for experiments can be obtained in the Supplementary Materials of this paper.

3.1. Information on Trajectory Data and Weights Update Method for Deep Learning

3.1.1. Information of Trajectory Data

The trajectory data are as follows: Extract 589 ADS-B cruising trajectory data of aircraft A320 from ODOPI to P101 (the solid red line in Figure 4) in sectors AR01, AR05, and AR04 of Guangzhou area in May 2019. The ADS-B data mainly include the flight call sign, aircraft type, UTC time (YMD h:m:s), longitude (^o), latitude (^o), query-normal elevation (QNE, m), ground speed (km/h). The air route of these specific trajectory data is shown in Figure 4.

The selected trajectory data are mainly adjacent to the nominal air route, while a few of them deviate from it. To make the obtained trajectory prediction model more suitable for actual application needs, as well as to guarantee its generalization ability, these trajectories have not been deleted. Part of the trajectory data is shown in Figure 5.

3.1.2. Weights Update Method for Deep Learning

In order to guarantee the performance comparability of various deep learning methods, the same loss function L is used by all machine learning models (LSTM, stateful-LSTM, BP neural network, 1D-ConvNet), which used to measure the MSE of the training data. Mathematically, the loss function L is:

$$\mathrm{minL}\left(W\right)=\frac{1}{\mathrm{I}·\mathrm{T}}{{\displaystyle \sum }}_{i=1}^{\mathrm{I}=412}{{\displaystyle \sum }}_{t=1}^{\mathrm{T}=2046}\parallel {\widehat{p}}_{it}-p_{it}\parallel ,p_{it}={\psi }_W\left(X_{it}\right)$$

(17)

$$\parallel {\widehat{p}}_{it}-p_{it}\parallel ={\left(x_{it}-{\widehat{x}}_{it}\right)}^2+{\left(y_{it}-{\widehat{y}}_{it}\right)}^2+{\left(z_{it}-{\widehat{z}}_{it}\right)}^2$$

(18)

All of the deep learning methods use the Adam weights update method and its initial weight generation method. The Adam weight update method [17,18] is shown in Algorithm 2.

Algorithm 2 weights update method based on Adam for deep learning

Input: the training data batch of the trajectory:$p_{it}$,${\widehat{p}}_{it}$,$i$∈ [0,412],$t$∈ [0,2046]

Output: trajectory prediction neuron network model ${\psi }_W$.

Hyper-parameters: training epochs K, β1 = 0.9, β2 = 0.999, η = 0.001.

Initialize: Initialize weights, W, as random numbers; Initialize intermediate variables $v_k$, $m_k$, $v{’}_k$, $m{’}_k$ as random numbers; $k$ = 0.

While $k$ < K:

L(W) = ${\widehat{p}}_{it}-p_{it}$

$m_k={\mathsf{\beta }}_1\cdot m_{k-1}+\left(1-{\mathsf{\beta }}_1\right)\cdot \frac{\alpha L}{\alpha W}$

$v_k={ \mathsf{\beta }}_2\cdot v_{k-1}+\left(1-{\mathsf{\beta }}_2\right)·\frac{\alpha L}{\alpha W}⨀\frac{\alpha L}{\alpha W}$

${m^{\prime }}_k=m_k/\left(1-{\mathsf{\beta }}_1\right)$

${v^{\prime }}_k=v_k/\left(1-{\mathsf{\beta }}_2\right)$

$W=W-\eta {m^{\prime }}_k/\sqrt{v{’}_k}$

$k$ = $k$+1

In Equations (17) and (18),$\mathrm{I}$ is the number of the trajectories in the training data set, T is the range of the prediction time span, $p_{it}$ is a vector that indicates the $t$-th point of the $i$-th trajectory, predicted by the prediction model ${\psi }_W$, ${\widehat{p}}_{it}$ is a vector that indicates the $t$-th actual point of the $i$-th trajectory. Equation (17) is the objective function for deep learning, and Equation (18) is the Euclidean distance between the predicted and true trajectory points. In Algorithm 2, ${\beta }_1$ and ${\beta }_2$ are hyperparameters setting the weight update iteration, where ${\beta }_1$ = 0.9, ${\beta }_2$ = 0.999, $\mathrm{L}$ is the loss function, $\eta$ is the step size, $W$ is the weight matrix of the hidden layers for deep learning, $k$ is the number of iterations, $m_k$, $v_k$, ${m^{\prime }}_k$, ${v^{\prime }}_k$ are intermediate variables that play the role of passing parameters, and $\frac{\alpha L}{\alpha W}$ is the partial derivative of the loss function L to the weight matrix $W$. The $⨀$ is the dot product symbol between two matrices.

3.2. Trajectory Prediction Experiments for the Candidate Model

3.2.1. Velocity Trend Extrapolation

Method introduction and experiment: the method of speed trend extrapolation predicts the position of the aircraft in the subsequent t second, based on the position and the ground speed of the aircraft’s first trajectory point ($v_0$). This method has the advantages of small computational complexity and little data dependence; the formula is:

$$p_{it}=p_{i\left(t-1\right)}+v_0$$

(19)

In Equation (19), $v_0$ is the ground speed of the aircraft’s first trajectory point, while $p_{it}$ is a vector that indicates the $t$-th point of the $i$-th trajectory predicted by the prediction model.

The MSEs of velocity trend extrapolation of the trajectory in the validation data set, changing with the time span of prediction (per second), are shown in Figure 6.

Results analysis: As a basic trajectory prediction method, speed trend extrapolation has a preferable prediction accuracy over a short period, with low computational complexity. It mainly benefits from the small speed-change rate of high-altitude and fixed-wing aircraft because there are always structured routes in the control areas, and the aircraft’s flight intentions and attitude are stable over a short time span of prediction. Its shortcoming is also obvious: once the aircraft changes its movement intention, such as normal route-turning or speed adjustment, the prediction error by speed trend extrapolation will increase rapidly. Therefore, this method is not suitable for medium- or long-term trajectory prediction because of the change in the aircraft intention.

3.2.2. BP Neural Network

Method introduction: the BP neural network is one of the classic deep learning methods. Since the data processing method, input and output data structure, loss function and weights update method have been explained earlier in Section 2.3 and Section 3.1, this section focuses on a description of the structure of the BP neural network and the comparison between the prediction performance of the BP neural network using different parameters. The BP neural network of this paper is a fitting model composed of error back-propagation layers, sigmoid, and ReLU activation functions [18]; the basic structure of the BP neural network prediction model is shown in Figure 7.

Excluding the loss function and weight update method, the BP neural network structure parameters mainly include: (1) the depth of the BP neural network (layers)—the greater the number of layers, the stronger the fitting performance of the neural network, but the weaker the generalization performance. (2) The number of neurons in each layer (neuron number)—the greater the size of the neuron number, the stronger the fitting performance of the neural network, but the weaker the generalization performance. (3) The training epochs of early stopping (training epochs)—in order to avoid overfitting, the generalization performance of the neural network is usually guaranteed by controlling the number of training epochs. Since it is impossible to test all of the possible structures of the BP neural network as the time needed for training a deep learning model is so long, this paper uses a random search, which tests several structures of the BP neural network to determine the best hyperparameters, as shown in

Table 2. The hyper-parameters of the BP neural network structure.

Name	Layers	Neuron Number	Training Epochs
BPmodel1	1	200	150
BPmodel2	1	300	150
BPmodel3	1	400	200
BPmodel4	2	100	150
BPmodel5	2	200	200
BPmodel6	2	300	250
BPmodel7	3	100	150
BPmodel8	3	200	150
BPmodel9	3	300	200

.

The MSE of the prediction time span (per minute) is used for comparing data visually because there are too many data sample points and fluctuations in the data, which are converted by the mean value of MSE per second. The prediction performance of the BPmodel1–BPmodel9 neural network (per minute) is shown in Figure 8.

Figure 8 is used to observe the generalization performance of the different models. It can be seen that BPmodel1 and BPmodel4 are under-fitting due to the lack of neurons, while BPmodel5 demonstrates the best prediction performance. We chose a prediction model with a dual hidden layer and epochs of 200 as the optimal model (BPmodel5, red dotted line in Figure 8). The MSEs of the optimal BP neural network in the validation set change with the time span of prediction (per second), as shown in Figure 9.

Results analysis: The prediction performance of the BP model gets worse as the time span of the prediction increases. Firstly, this is because the BP neural network structure cannot memorize time-series information according to its structure. Secondly, the prediction styles in this paper (one-step, bootstrap, global trajectory prediction styles) lead the MSE of the BP model to accumulate along with the time span of prediction getting longer. However, compared to the method of velocity extrapolation, it can be seen that the BP neural network has the ability to learn the non-linear features of trajectories, which suppress the rapid growth of MSE according to the time spans of the predictions getting longer.

3.2.3. LSTM

Method introduction: The LSTM model is a recurrent neural network (RNN) composed of LSTM memory units, so it has the same basic structure as an RNN but it is more complex than that. The gate structures are composed of the sigmoid function and the tanh function in the LSTM unit. The gate structures (f_t, i_t, c’_t, and o_t in Figure 10) set in the LSTM unit can give it an advantage in the processing of sequence data. The good performance of LSTM in processing time-series data and memorizing key information has been demonstrated previously by many scholars [19,20,21,22,23]. LSTM neurons consist of a hidden gate, input gate, output gate, and forget gate. The basic structure of the LSTM regression model [24] is shown in Figure 10.

Experiment: The input and output data structure, loss function and weights update method, and other information have been described in Section 2 and Section 3.1. In order to explore the optimal LSTM model for trajectory prediction, our experiments mainly adjust the number of LSTM units and the training epochs of early stopping. According to the input data structure, the trainable parameters that each LSTM model contains comprise:

$$train\_num=4\times \left(size+hidden\right)\ast hidden+hidden$$

(20)

where $size$ is the input data dimension, while $hidden$ is the number of LSTM units. In order to avoid over-fitting, the number of LSTM units is reduced appropriately. Several hyper-parameters of the LSTM structure that were tested in experiments were the same as in the random search method in Section 3.2.2; these are shown in

Table 3. The hyper-parameters of the LSTM structure.

Name	Number of Neurons	Training Epochs
LSTMmodel1	128	250
LSTMmodel2	256	250
LSTMmodel3	384	250
LSTMmodel4	128	350
LSTMmodel5	256	350
LSTMmodel6	384	350

.

After learning from the training data, the prediction performance of the models in Table 3 is shown in Figure 11.

It can be seen from Figure 11 that the prediction performance of the LSTM is better than the BP, and the model of LSTMmodel1 is under-fitting (compared to model4) while the model of LSTMmodel5 is over-fitting compared to LSTMmodel2. Chose LSTMmodel2, with the best performance, as the optimal model. The MSE of the optimal LSTM model in the validation data changes with the time span of prediction is in Figure 12.

Results analysis: the MSE of the LSTM model remained at a low level when the time span of prediction is less than 400s but, as the time span of prediction becomes longer, its prediction performance tends to deteriorate significantly (MSE increases significantly). To explore the reason for this, we changed the LSTM regression model to the LSTM multi-class classification model, where the output is the probability of the aircraft position. The structure of the LSTM regression model and the LSTM multi-class classification model is similar and the only difference between them is the structure of the output data [24]. The structure of this multi-class classification model is shown in Figure 13.

The output of the LSTM multivariate classification model is the probability of the predicted position. We trained this multi-classification model, then predicted the trajectory position probability with a prediction time span of 500 s and 1000 s. The Monte Carlo method [25] was used to output 5000 trajectory positions, then the Gaussian kernel function was used to calculate the probability density of each position. The probability of the predicted position is shown in Figure 14.

It can be seen in Figure 13 that although the LSTM model has excellent time-series data memory capabilities, its predictive ability will still be restricted by the chaotic phenomenon [26,27,28,29] in the process of trajectory prediction. The chaotic performance is this: as the time span of prediction becomes longer, the output probability density of the predicted trajectory position gradually becomes uniform. This phenomenon stems from the uncertainty of the aircraft’s operation. The chaotic phenomenon makes it difficult for most prediction models to predict the position of the aircraft over a long time span accurately. Therefore, except for the systematic accumulation of prediction errors in the single-step, bootstrap, and global prediction style in this paper, the chaotic phenomenon of the aircraft is also an important reason for the decrease in prediction performance as the time span of the prediction increases.

3.2.4. Stateful-LSTM

Method introduction: the basic structure of the stateful-LSTM is as same as the LSTM model in Section 3.2.3. “Stateful” means that the last state for each sample of index t in a batch will be used as the initial state for the sample of index t in the following training batch, which preserves the dependencies between trajectory series data [30]. There is previous research [30,31] that shows that the stateful-LSTM has advantages in processing long time-series data, especially when there are strong dependencies between adjacent time-series data sets.

Experiment: Since the basic structure of LSTM and stateful-LSTM are similar, the structure parameters for stateful-LSTM experiments that are used are shown in Table 3 as well. After learning from the training data, the prediction performance of each stateful-LSTM model in Table 3 is shown in Figure 15.

We chose the optimal model of stateful-LSTMmodel2 as the one with the best prediction performance. The MSE of the optimal stateful-LSTM in the validation data changes according to the time span of prediction (per second), as shown in Figure 16.

Results analysis: when using the training data to train the stateful-LSTM model, the training epochs to achieve the minimum MSE in the validation data are more than in LSTM (Figure 17). The prediction performance of the stateful-LSTM is more sensitive to the training epochs, due to the transfer of state parameters in the stateful-LSTM. Users of the stateful-LSTM can learn from the training data in more detail, meaning that the trajectory prediction performance of the stateful-LSTM is better than the LSTM. However, the stateful-LSTM model is also more prone to overfitting than the LSTM model in experiments. The generalization ability declines when there are too many training times chosen or the model structure is too complex, which is manifested as an increase in the loss of the validation set at the later epochs of training. This is also the reason why models with different hyperparameters have been designed in this paper for comparison. Comparing Figure 9, Figure 12 and Figure 16, it can be seen that the stateful-LSTM has a good suppression effect on the trend of the prediction error, increasing with the time span of the prediction. Generally, the stateful-LSTM offers a preferable performance for medium- and long-term trajectory prediction.

3.2.5. The 1D-ConvNet

Method introduction and experiment: the 1D-ConvNet is a convolutional neural network, one of the deep learning methods. The basic neural network unit of the 1D-CNN is also the convolution-pooling mode. Although the 1D-CNN has a similar structure to traditional convolutional neural networks, it is more compact than a 2D-CNN. A 1D-CNN can be applied successfully to the time series analysis of trajectory data. The 1D-CNN in this paper consists of a convolution-max pooling layer and a fully connected output layer; the basic structure of the 1D-ConvNet is shown in Figure 18 [32]. The input and output data structure, loss function, weights update method, and the prediction mode have been described earlier in Section 2 and Section 3.1.

Experiment: excluding the loss function and weights update method, the 1D-ConvNet structure parameters mainly include the number of convolution max-pooling layers, the number of convolution and max-pooling units in each layer, and the training epochs of early stopping. Since it is impossible to obtain an optimal structure accurately, this paper only tests those structure hyper-parameters that are usually better, as shown in

Table 4. The hyper-parameters of the 1D-ConvNet structure.

Name	Layers	Conv-Max-Pooling Units	Training Epochs
1D-Convnet1	1	128	150
1D-Convnet2	1	128	200
1D-Convnet3	1	256	200
1D-Convnet4	2	128	150
1D-Convnet5	2	128	200
1D-Convnet6	2	256	300

.

After the model’s learning from the training data, the prediction performance of each 1D-ConvNet model in Table 4 is shown in Figure 19.

It can be seen in Figure 19 that the models of 1D-ConvNet1 and 1D-ConvNet2 show the phenomenon of under-fitting compared to the 1D-ConvNet6, according to the prediction performance of these models. The prediction model1D-Convnet4 is selected as the optimal model (the red dotted line in the figure above) because it offers a better prediction performance. The MSE of the optimal 1D-ConvNet in the validation data, changing with the time span of prediction (per second), is shown in Figure 20.

Results analysis: The performance of the 1D-ConvNet decreases as the time span of prediction gets longer, as with other trajectory prediction methods; the trajectory prediction performance of the 1D-ConvNet is better than the BP neuron network because of its model structure, which has the ability to learn features from time-series data. The 1D-ConvNet method has good prediction performance but it is inferior to the prediction performance of the optimal stateful-LSTM model, overall.

3.2.6. Kalman Filter (KF)

Method introduction and experiment: the prediction model based on KF is widely used for short-term trajectory prediction. There are two main hypotheses while using the KF method to predict the trajectory: firstly, the relationship between the independent variable and the dependent variable in the system is linear. Secondly, the systematic error and measurement error demonstrate the features of Gaussian distribution [24]. Generally, the prediction model of KF is:

$$\mathrm{minMSE}=\frac{1}{\mathrm{I}·\mathrm{T}}·{{\displaystyle \sum }}_{i=1}^{\mathrm{I}=412}{{\displaystyle \sum }}_{t=1}^{\mathrm{T}=2046}\parallel {\widehat{p}}_{it}-p_{it}\parallel$$

(21)

$$p_{i,t+1}=p_{i,t}⨀{\mathit{A}}_t+{\mathit{B}}_t⨀{\mathit{X}}_t+w_t$$

(22)

Equation (21) shows the objective function of the KF trajectory prediction model, where $p_{it}$ is a vector with the size of 1 × 3 that indicates the $t$-th position of the $i$-th trajectory, predicted by the prediction model based on KF, and ${\widehat{p}}_{it}$ is a vector that indicates the t-th actual position of the i-th trajectory. In Equation (22), ${\mathit{X}}_t$ is a vector with a size of 1 × 8 that includes the parameters described in Equation (2), while ${\mathit{A}}_t$,${\mathit{B}}_t$ are the vectors of parameters to be solved. $w_t$ is the error that fits with the Gaussian distribution. The $⨀$ is the dot product symbol between two matrices.

After fitting using the training data set, the MSE of the Kalman filter model in the validation data changes with the time span of prediction (per second), as shown in Figure 21.

Results analysis: The model based on the Kalman filter is widely used in the field of trajectory prediction. It can be seen from Figure 21 that, compared with the trajectory prediction performance of other methods, the prediction accuracy of the Kalman filter is satisfactory in terms of the short-term prediction process. In addition, the Kalman trajectory prediction model is more suitable for the mode of single-step prediction, while the performance of global trajectory prediction, especially for the mid- to long time span of prediction is not good under the mode of prediction used in this paper. Although the Kalman filter successfully extracts the discipline of trajectory changing, its prediction error accumulates very obviously with the increase in the prediction time span. This conclusion is close to some of the previous research results based on the Kalman filter [3,8].

3.2.7. Flight Plan Interpolation (FPI)

Method introduction and experiment: firstly, we interpolate the flight plan according to the estimated time and waypoint of the original flight plan to obtain the corresponding trajectory point of each second [33,34]. Then, we obtain the flight plan corresponding to each trajectory in the data set, extracting the waypoints ODOPI, XEBUL, LMN...P101 in the flight plan, and linearly interpolate the estimated time of each trajectory point. Mathematically, the flight plan interpolation method can be expressed as:

$$p_t={\widehat{p}}_{t_k}+\frac{{\widehat{p}}_{t_{k+1}}-{\widehat{p}}_{t_k}}{t_{k+1}-t_k}\left(t-t_k\right)$$

(23)

where $p_t$ is the vector of the predicted trajectory position at t s, $t_k$ is the estimated passing time of the $k$-th flight plan waypoint, and ${\widehat{p}}_{t_k}$ is the vector consisting of the longitude, latitude, and altitude of time $t_k$ in the flight plan.

The MSE of the flight plan interpolation in the validation data changes with the time span of prediction (per second), as shown in Figure 22.

Results analysis: using flight plan interpolation to predict the trajectory at different time spans does not require the input of aircraft speed, wind direction, and other parameters. Its prediction performance is not greatly affected by the time span of the prediction because the flight plan reflects the long-term movement characteristics of the aircraft during normal operation. Taking the data in this paper as an example, the average MSE predicted by the flight plan interpolation method is 7.884 km, which means that the flight plan interpolation method generally has an error of about 7.9 km when predicting the position of the aircraft, which is unacceptable in short-term trajectory prediction. However, it shows a better prediction performance than other methods mentioned above in long-term trajectory prediction because the pilot and controller jointly control the RTA time to keep the flight on schedule.

4. Hybrid Trajectory Prediction Model

Section 4.1 describes an example of constructing a hybrid trajectory prediction model based on the above seven trajectory methods. Section 4.2 details the trajectory prediction experiments and performance analysis of the hybrid trajectory prediction model, to demonstrate its rationality.

4.1. Hybrid Prediction Model Constructing

After the experiments using the above-mentioned trajectory prediction methods, it can be seen that each method has its advantages in terms of the different time spans of prediction. Therefore, the optimal model under the corresponding prediction span is used to predict the trajectory, and a hybrid trajectory prediction model is further proposed. The hybrid trajectory prediction model can perform short-term, medium- and long-term trajectory predictions while theoretically offering a good prediction performance.

The data of the short-term MSE, changing with the time span of prediction (per second) and the data of the mid- to long-term 34-minute MSE, changing with the time span of prediction (per minute) of different prediction methods, are shown in

Table 5. Comparison of the MSE of different prediction methods in the first 45 s.

Time Span (s)	VTE	BP	LSTM	KF	1D-CNN	S-LSTM	Time (s)	VTE	BP	LSTM	KF	1D-CNN	S-LSTM
1	0.052	--	--	--	--	--	23	0.113	0.260	0.135	0.199	0.127	0.093
2	0.023	--	--	--	--	--	24	0.131	0.242	0.154	0.185	0.156	0.137
3	0.028	--	--	--	--	--	25	0.097	0.262	0.138	0.200	0.142	0.129
4	0.034	--	--	--	--	--	26	0.103	0.272	0.151	0.207	0.161	0.106
5	0.036	0.175	0.139	0.177	0.109	0.081	27	0.107	0.274	0.167	0.224	0.154	0.127
6	0.061	0.172	0.136	0.157	0.102	0.076	28	0.142	0.275	0.155	0.182	0.161	0.125
7	0.045	0.177	0.097	0.164	0.113	0.066	29	0.102	0.277	0.156	0.183	0.172	0.129
8	0.046	0.163	0.136	0.133	0.103	0.081	30	0.110	0.284	0.148	0.210	0.168	0.125
9	0.063	0.192	0.100	0.165	0.133	0.088	31	0.120	0.278	0.168	0.200	0.147	0.127
10	0.062	0.208	0.130	0.162	0.118	0.071	32	0.140	0.293	0.188	0.212	0.170	0.139
11	0.048	0.202	0.124	0.149	0.117	0.085	33	0.145	0.302	0.170	0.203	0.176	0.128
12	0.083	0.222	0.113	0.199	0.154	00067	34	0.153	0.313	0.139	0.217	0.170	0.132
13	0.070	0.215	0.125	0.149	0.115	0.010	35	0.132	0.292	0.161	0.211	0.182	0.141
14	0.104	0.222	0.152	0.179	0.115	0.107	36	0.164	0.312	0.198	0.225	0.175	0.126
15	0.097	0.223	0.118	0.155	0.115	0.121	37	0.157	0.314	0.199	0.186	0.164	0.153
16	0.067	0.224	0.121	0.148	0.133	0.113	38	0.154	0.321	0.198	0.218	0.187	0.150
17	0.120	0.241	0.162	0.179	0.152	0.100	39	0.165	0.319	0.203	0.236	0.194	0.142
18	0.104	0.258	0.122	0.156	0.133	0.133	40	0.172	0.327	0.174	0.225	0.188	0.149
19	0.100	0.225	0.148	0.160	0.128	0.074	41	0.163	0.335	00168	0.232	0.189	0.172
20	0.107	0.263	0.159	0.171	0.134	0.104	42	0.151	0.343	0.201	0.207	0.186	0.168
21	0.116	0.257	0.139	0.198	0.144	0.131	43	0.169	0.371	0.198	0.259	0.167	0.145
22	0.073	0.272	0.131	0.185	0.136	0.135	44	0.184	0.331	0.182	0.264	0.211	0.170

and

Table 6. Comparison of the average MSE per minute of the different prediction methods.

Time Span (min)	VTE	FPI	LSTM	KF	1D-CNN	S-LSTM	Time (min)	VTE	FPI	LSTM	KF	1D-CNN	S-LSTM
1	0.05	8.06	0.13	0.07	0.08	0.04	18	9.04	7.72	12.72	17.29	7.69	3.28
2	0.09	8.12	0.32	0.10	0.12	0.08	19	18.84	7.84	13.80	19.93	8.07	3.65
3	0.29	8.11	0.66	0.39	0.33	0.13	20	28.86	7.99	14.27	21.24	8.59	4.79
4	0.75	7.59	0.98	0.52	0.76	0.31	21	40.41	8.24	15.39	22.58	9.56	5.80
5	1.08	7.25	1.59	0.66	1.04	0.49	22	53.90	8.23	16.02	24.09	10.88	6.99
6	1.54	7.42	2.04	0.85	1.38	0.67	23	63.17	8.00	16.94	25.39	12.19	8.05
7	2.58	7.61	2.50	1.08	1.73	0.93	24	78.45	7.81	17.80	26.62	13.09	9.07
8	3.65	7.78	3.02	1.39	2.31	1.06	25	89.36	7.60	18.51	27.98	14.20	10.03
9	3.85	7.97	4.03	1.69	2.85	1.31	26	94.75	7.65	19.50	29.49	15.46	11.29
10	4.41	7.91	5.22	2.06	3.54	1.43	27	--	7.81	20.17	30.84	16.69	12.41
11	4.27	7.89	6.12	2.34	4.12	1.55	28	--	7.90	21.41	32.40	17.96	13.34
12	4.49	7.86	6.77	3.25	4.67	1.83	29	--	8.03	22.27	34.16	19.64	13.98
13	5.91	7.95	7.53	4.64	5.21	1.94	30	--	8.17	23.13	35.59	20.69	14.55
14	7.38	8.04	8.24	6.21	5.66	2.21	31	--	8.19	24.11	37.46	23.73	15.25
15	8.24	8.10	9.42	8.86	6.13	2.47	32	--	7.25	25.26	39.13	24.84	16.86
16	6.61	8.14	10.50	11.71	6.61	2.75	33	--	7.94	26.47	40.72	26.29	17.67
17	4.71	8.01	11.58	14.50	7.18	2.95	34	--	7.72	27.43	42.36	27.64	18.49

and Figure 23 and Figure 24.

From the MSE comparison results of each prediction method in the first 45 s in Figure 23 and Figure 24 (the FPI model is not displayed in Figure 23 due to its large MSE), it can be seen that the speed trend extrapolation trajectory prediction method is better than other trajectory prediction methods in the first 25 s, while the stateful-LSTM has better prediction performance from 25 s to 1320 s (22 min), and the flight plan interpolation has relatively better prediction performance after 22 min. In addition, after 22 min, the prediction performance of all the above prediction models is close to or inferior to the flight plan interpolation prediction performance. It can be concluded that the prediction performance of all prediction methods after 22 min is not good.

The hybrid trajectory prediction model is composed of various advantageous prediction methods in different time spans of prediction (as in Figure 25). In terms of the different time spans of prediction, 0–25 s uses the speed-trend extrapolation prediction method, 25 s to 1320 s uses the stateful-LSTM prediction method, and after 22 min, the flight plan interpolation prediction method is used. The MSE of the hybrid trajectory prediction model, changing with the time span of prediction (per second), is shown in Figure 26 and Figure 27.

Results analysis: The hybrid model combines the advantages of these prediction methods, and it can be seen that its MSE over time is lower than any other single prediction method. It can be seen that even in a prediction time span of more than 1500 s, the MSE is only about 8 km, while this hybrid trajectory prediction model has the calculation complexity of flight plan interpolation. What is more, this hybrid trajectory prediction method is suitable for both short-term and long-term trajectory prediction. This hybrid trajectory model also has low algorithm complexity and good prediction performance when combining the flight planning model to predict long-term trajectories.

At the same time, the maximum predictable time span of a specific trajectory could be determined in the process of building the hybrid model. When the time span of prediction exceeds this, the prediction methods will lose their good prediction performance. Their prediction performance will be close to or even inferior to the flight plan interpolation prediction method, while the computational complexity is still higher than FPI.

4.2. Model Verification

For the model verification, Extract 10 of the new trajectory data (testing data set) from the sectors of the Guangzhou area that have not been included in the training data set and validation data set is used for further verification of the performance of the hybrid trajectory prediction model. Then, we used the hybrid prediction model, stateful-LSTM, LSTM, 1D-ConvNet, Kalman filter, and BP neural network to predict it and test the prediction performance. One of the predicted trajectories of the hybrid trajectory prediction model in the sectors of the Guangzhou area is shown in Figure 28.

It can be seen from Figure 28 that the prediction result of the hybrid trajectory prediction model does not produce an unacceptable trajectory prediction error as the time span of the prediction increases. The hybrid model makes full use of the advantage of velocity trend extrapolation that benefited from aircraft state stability in the short-term trajectory prediction, the advantage of memory ability regarding the time series data of the stateful-LSTM model in mid-to-long-term trajectory prediction, and the advantage of the stability of FPI in long-term trajectory prediction. These advantages bring the predicted trajectory obtained by the hybrid prediction model closer to the actual trajectory. The MSEs of various methods in the test trajectory data set are shown in

Table 7. Performance comparison of prediction methods.

Name	MSE/km	Deviation/%
Hybrid	4.26	--
BP	20.13	78.80%
LSTM	13.86	69.26%
stateful-LSTM	8.72	51.46%
1D-Convnet	12.15	64.94%
KF	16.01	73.39%
FPI	6.32	32.60%

.

It can be seen from Table 7 that the prediction accuracy of the hybrid trajectory prediction model (MSE) is significantly smaller than any other original trajectory prediction method, and this method can keep a smaller MSE in the different time spans of trajectory prediction. Finally, this trajectory prediction method and the hybrid trajectory prediction model are theoretically applicable to other different sectors, different air routes, and different aircraft types, with quite a good generalization performance.

It should be pointed out that the trajectory data used in this paper are the cruising aircraft data at a specific altitude (9200 m). Although the aircraft’s altitude is used as input and output, the aircraft’s altitude is a fixed value. Unlike the approaching aircraft, the cruising aircraft do not show a clear intention in their altitude. The aircraft’s altitude adjustment completely follows the controller’s order, and the controller’s order may be related to the resolution of a conflict [35,36]. Therefore, it is unreasonable to use machine learning methods to predict the trajectory of a single aircraft in a vertical direction. For this reason, this research is more inclined to trajectory prediction in the dimensions of longitude, latitude, and time. The research regarding trajectory prediction in the vertical dimension should consider the interactions between aircraft theoretically, which can offer a further research topic in trajectory prediction.

5. Conclusions

(1)

For the trajectory data set in this paper, within the time span of prediction of 0–25 s (short-term), the velocity trend extrapolation trajectory prediction method has a minimum MSE, and within the prediction time span of 25 s-22 min (mid- to long-term), while the flight procedure interpolation method has a minimum MSE after 22 min (long-term). Since the experiment only uses the trajectory data set of a single flight route, this conclusion has some limitations (Appendix A: an additional case). However, when training and predicting by different methods, this paper tries to control these variables to guarantee the comparability of the prediction performance among different methods. Therefore, the performance comparison results of the various methods obtained in this paper are still worth referring to after changing the trajectory data set.

(2)

A method for constructing a hybrid trajectory prediction model has been proposed. The hybrid trajectory prediction model is a combination of models, with predictive advantages in different time spans of prediction. This hybrid trajectory prediction model combines the advantages of each sub-model in the different time spans of prediction, so the hybrid model’s performance will at least be no worse than any of the individual models that comprise it.

(3)

By comparing the prediction performance of other prediction models with the performance of flight plan interpolation, the maximum accurate time span of prediction was initially proposed. The performance of the trajectory prediction method beyond this time span of prediction will be close to or even inferior to the performance of flight plan interpolation. In other words, for a specific prediction model and data set, the trajectory prediction beyond this time span of prediction does not offer the significance of complex calculation.