2.1. Traditional predictive model

From a mathematical point of view, the prediction first solves the problem of describing the known elements by establishing an approximate description function, and then solves the extrapolation of this approximate function. Dynamic platform navigation should not only measure parameters such as position, speed, acceleration of current time and past times, but also use this data to predict the state, trend and trajectory of the platform for the next time period. It is necessary to build an approximate forecasting function which can reflect the change of the trajectory by using previous and current trajectory status information (Morzy, Reference Morzy2006; Song and Barabási, Reference Song and Barabási2010; Qiao et al., Reference Qiao, Jin, Han, Tang and Gutierrez2015a). Therefore, the optimal approximate function of the known state and the extrapolated model of the trajectory at the next time period are key to solving the forecasting problem.

2.1.1. Dead reckoning algorithm

The dead reckoning algorithm was originally widely used in navigation and positioning of ships and vehicles (Hwang, Reference Hwang2002; Retscher and Kealy, Reference Retscher and Kealy2006). The platform is considered to be in linear motion, and the next time position of the platform can be predicted when a priori information and the course angle of the current time position are known. The dead reckoning algorithm is illustrated in Figure 1.

Figure 1. A schematic diagram of the dead reckoning algorithm.

Assuming the position $p_{t_{n}} \lpar x_{t_{n}}\comma \; y_{t_{n}}\rpar $ at time t _n, previous velocity $v_{t_{n}}$ and course angle $\theta_{t_{n}}$ are known, the position of $p_{t_{n+1}} \lpar x_{t_{n+1}}\comma \; y_{t_{n+1}}\rpar $ at time t _n+1 can be calculated by the dead reckoning algorithm as Equation (1).

(1)$$\left\{\matrix{x_{t_{n+1}}=v_{t_n}\lpar t_{n+1}-t_n\rpar \cos\lpar \theta_{t_n}\rpar +x_{t_n} \cr y_{t_{n+1}}=v_{t_n}\lpar t_{n+1}-t_n\rpar \sin\lpar \theta_{t_n}\rpar +y_{t_n}}\right. $$

In actual dead reckoning applications, information such as the velocity $v_{t_{n}}$, course angle $\theta_{t_{n}}$, etc, may be provided by the on board sensors such as a speedometer, magnetic compass and GNSS receiver.

2.1.2 Least squares fitting method

In the fitting process of the navigation trajectory, the least squares fitting method is used to solve the platform trajectory fitting and then the next time position data is obtained by extrapolation on the fitting curve by the inertia principle (Wang et al., Reference Wang, Wang, Xiong and He2016), as shown in Figure 2. In Figure 2, according to the six known coordinate points, a curve is fitted by the least squares criterion, and then the position of the next moment is extrapolated on the curve.

Figure 2. Least squares fitting-extrapolation algorithm illustration.

Given the position coordinates $\lpar x_{t_{1}}\comma \; y_{t_{1}}\rpar \comma \; \lpar x_{t_{2}}\comma \; y_{t_{2}}\rpar \comma \; \lpar x_{t_{3}}\comma \; y_{t_{3}}\rpar \comma \; \ldots\comma \; \lpar x_{t_{n}}\comma \; y_{t_{n}}\rpar $ at time t ₁, t ₂, t ₃, …, t _n, the coordinates at time t _n+1 can be extrapolated through the least squares fitting curve. The coordinates x and y can be calculated by Equations (2) and (3):

(2)$$\left\{\matrix{\min X = \sum\limits_{i=1}^n \, \lpar \, f\lpar t_i\rpar -x_{t_i}\rpar^2 \cr \min Y = \sum\limits_{i=1}^n \lpar \varphi \lpar t_i\rpar -y_{t_i}\rpar^2}\right. $$

(3)$$\left\{\matrix{x_{t_{n+1}}=f\lpar t_{n+1}\rpar \cr y_{t_{n+1}}=\varphi \lpar t_{n+1}\rpar }\right. $$

where f(t _i), φ(t _i) are the polynomial fittings of the x-axis and the y-axis. The two fitting curves are used to extrapolate the coordinate $\lpar x_{t_{n+1}}\comma \; y_{t_{n+1}}\rpar $ in the next time period.

The least squares extrapolated model has been widely used, but in some applications, such as when a dynamic platform moves at high speed or when turning through large angles, it is difficult for the extrapolated model to reflect the movement of the platform. In order to solve this problem, a new approximate model and a new solution are proposed in this paper. The key idea is to take advantage of the information such as the angle of the platform at the current time to predict the future trajectory of the dynamic platform, so it can reflect dynamic platform trends better.

2.2. Generalised extension extrapolation mathematical model

The trajectory prediction model of a dynamic platform is divided into two categories: fitting and interpolation. The representative mathematical models include the least squares fitting polynomial function, the power polynomial interpolation function, the spline function containing the derivative value and the Kalman filtering method. A generalised extension extrapolation model that combines interpolation and fitting two approximation elements is proposed in this paper. The method introduces the angle information which is effective for prediction and it can improve the accuracy and performance of the prediction.

The generalised extension extrapolated model proposed in this paper is different from the traditional extrapolated model. It takes the advantages of dead reckoning and least squares fitting extrapolation and takes into account the a priori information of the real-time platform position information such as angle information to fit the extrapolation. At the same time, real-time angle information is added as a constraint condition, which makes the new model more robust against external changes and it can respond quickly to the change of the latest data points, especially when adapting to a large change in the angle of the vehicle's trajectory. In this model, the random error of the motion parameters can be filtered by fitting, and the constraints of interpolation can make full use of the important motion parameters to improve the accuracy of extrapolation, so it can improve the accuracy and performance of prediction.

The specific mathematical formula of the generalised extension extrapolated model (Shang et al., Reference Shang, Cheng, Sheng, Shi and Yue2016) using a known trajectory in the model is fitted with the quadratic curve under the minimum error criterion, and then the angle information of the current moment is added as the constraint condition.

(4)$$\left\{\matrix{\min I = \sum\limits_{i=m}^{n} \, \lpar ax_i^2+bx_i+c-y_i\rpar^2 \hfill \cr subject\, \, to\, \, 2ax_n+b=\tan\alpha_n \hfill \cr}\right. $$

where I is the objective optimisation function; x _i, y _i are the position of measurement point i; (x _m, y _m), (x _m+1, y _m+1), …, (x _n, y _n) are all fitting points in the model; α _n is the angle value (0 ≤ α ≤ 2; π) between velocity direction of the point n and positive x-axis and a, b and c are the coefficients of the extrapolated polynomial.

The optimal solution of the parameters a, b and c can be obtained using Equation (4). Using the fitting curve and the velocity information, the coordinates of the prediction point (x _n+1, y _n+1) can be obtained using Equation (5):

(5)$$\left\{\matrix{y_{n+1} = ax_{n+1}^2+bx_{n+1}+c \hfill \cr \lpar x_{n+1}+x_n\rpar^2+\lpar y_{n+1}+y_n\rpar^2=\lpar v_nT\rpar^2}\right. $$

where v _n is the velocity value and T is the time interval.

2.3. The motion model for trajectory prediction

From Equation (4), the approximate polynomial of the generalised extension extrapolated model is a quadratic polynomial and the quadratic polynomial function is limited: in the curve shown in Figure 3, the expression of the quadratic function is fitted easily because the curve has obvious analytical expressions, but in Figure 4, the curve does not have the form of an analytic function; it is more difficult to fit the curve when the quadratic function is used. To solve the problem, a coordinate transformation is used to make the platform trajectory meet the approximate quadratic polynomial. The coordinate transformation makes full use of the symmetry of the quadratic function about the pole and avoids the extreme solution problem when the tangent function value is infinite in the mathematical model.

Figure 3. Quadratic function basic graphics.

Figure 4. The quadratic polynomial cannot approach the situation.

The method of coordinate transformation is as follows:

Take the current position O ^′ of the platform as the new coordinate origin, define the tangent direction of the arc that is formed by the prior motion trajectory and motion direction current time as the x ^′-axis. Define the angle of the bisector formed by the prior point and the direction of motion as the y ^′-axis direction. See Figure 5.

Figure 5. Coordinate transformation diagram.

In Figure 5, ABO ^′C is the platform trajectory curve. Point O ^′ is the current point. The line O ^′C is the current motion direction. A and B are a priori points. Define the tangent direction of curve BO ^′C as the x ^′-axis. Define the angle bisector direction of curve ∠BO ^′C as the y ^′-axis.

The position of the next point C can be extrapolated from the points A, B, O ^′. Assuming that the coordinates of A, B, O ^′ in the original coordinates are (x _n−2, y _n−2), (x _n−1, y _n−1), (x _n, y _n). The angle between BO ^′ and the x-axis is α ₁. The angle between O ^′C and the x-axis is α ₂ and the angle between the x ^′-axis and the positive direction of the x-axis is α _x.

When α ₂ − α ₁ > 0, the angle between the x ^′-axis and the x-axis positive direction is:

(6)$$\alpha_x=\displaystyle{{1}\over{2}}\lpar \alpha_2-\alpha_1\rpar +\alpha_{1} $$

When α ₂ − α ₁ < 0, the angle between the x ^′-axis and the x-axis positive direction is:

(7)$$\alpha_x=\displaystyle{{1}\over{2}} \lsqb \lpar \alpha_2-\alpha_1\rpar +2\pi\rsqb +\alpha_1 $$

The coordinates of A, B, O ^′ in the new coordinate system are obtained (Ye, Reference Ye2006) by the coordinate transformation formula:

(8)$$\left[\matrix{x^{\prime}_{n-2} \cr y^{\prime}_{n-2}}\right]= \left[\matrix{x_{n-2}-x_n & y_{n-2}-y_n \cr y_{n-2}-y_n & x_{n}-x_{n-2}}\right]\left[\matrix{\cos\alpha_x \cr \sin\alpha_x}\right]$$

(9)$$\left[\matrix{x^{\prime}_{n-1} \cr y^{\prime}_{n-1}}\right]= \left[\matrix{x_{n-1}-x_n & y_{n-1}-y_n \cr y_{n-1}-y_n & x_{n}-x_{n-1}}\right]\left[\matrix{\cos\alpha_x \cr \sin\alpha_x}\right]$$

The point O ^′ is the new coordinate origin, so the coordinate of the point O ^′ is $\lpar x^{\prime}_{n}\comma \; y^{\prime}_{n}\rpar =\lpar 0\comma \; 0\rpar $, then substitute $\lpar x^{\prime}_{n-2}\comma \; y^{\prime}_{n-2}\rpar $, $\lpar x^{\prime}_{n-1}\comma \; y^{\prime}_{n-1}\rpar $, (x ^′_n, y ^′_n) into Equation (4) to solve the approximate quadratic polynomial. The position of the next time point C is extrapolated using the fitted polynomial and the velocity information from the current time, and the algorithm uses Equation (10) to solve the coordinates $\lpar x^{\prime}_{n+1}\comma \; y^{\prime}_{n+1}\rpar $ of point C.

(10)$$\left\{\matrix{y^{\prime}_{n+1}=a x_{n+1}^{\prime 2}+bx^{\prime}_{n+1}+c \hfill \cr \lpar x^{\prime}_{n+1}-x^{\prime}_n\rpar^2+\lpar y^{\prime}_{n+1}-y^{\prime}_n\rpar^2=\lpar v_nT\rpar^2}\right. $$

In Equation (10), v _n is velocity at time n and T is the time interval between time n and n + 1. Then the coordinates of point C in the original coordinate system can be obtained using the coordinate inverse transformation, as shown in Equation (11).

(11)$$\left[\matrix{x_{n+1} \cr y_{n+1}}\right]=\left[\matrix{x^{\prime}_{n+1} & -y^{\prime}_{n+1} \cr y^{\prime}_{n+1} & x^{\prime}_{n+1}}\right]\left[\matrix{\cos\alpha_x \cr \sin\alpha_x}\right]+\left[\matrix{x_n \cr y_n}\right]$$

3.1. Algorithm simulation

The simulation was based on the MATLAB-2012a simulation platform. First, a theoretical trajectory was defined and a random error was added to the theoretical trajectory to give initial positioning points. Epoch data was then extrapolated by using the generalised continuation extrapolated model, least squares fitting model and dead reckoning algorithm. The simulation results are shown in Figures 6 and 7. In the figures, the red “*” sign is the positioning coordinate. The dotted red line “ + ” is the positioning trace curve of the generalised extended prediction model, the straight blue line “o” is the positioning trace curve of the least squares fitting extrapolation model and the straight black line “⋄” is the positioning trace curve of the dead reckoning model.

Figure 6. Sine trajectory simulation results.

Figure 7. Simulation results of “” shape trajectory.

In Figures 6 and 7, the generalised extension extrapolated model proposed in this paper, the least squares extrapolated model and dead reckoning model are used to predict the next future position. According to Figures 6 and 7, using the least squares fitting model, the extrapolation data using the principle of inertia has an obvious overshoot phenomenon when a rapid change in platform motion direction occurs. The dead reckoning algorithm simplified the platform motion model to linear motion, so the prediction accuracy depends entirely on the positioning accuracy of the current point and time interval between two epochs. The positioning error of the current point is totally inherited in the prediction. The generalised extension extrapolated model can not only respond extremely quickly to the change of platform motion, but also makes full use of prior data, so it has a better accuracy for extrapolation.

In order to further analyse the performance of the algorithm, the error curve of prediction values from different algorithms of two motion trajectories are shown in Figures 8 and 9. In the figures, the dotted red line “ + ” is the positioning error of the generalised extended prediction model, the straight blue line “o” is the positioning error of the least squares fitting extrapolation model, and the straight black line “⋄” is the positioning error of the dead reckoning model.

Figure 8. Sine track prediction trace error curve.

Figure 9. Prediction error of “” shape trajectory.

For the above two trajectories, the mean errors of the different models are shown in Table 1.

Table 1. Different extrapolated model prediction error.

Table 1 shows that the prediction error of the generalised extension extrapolated model is less than 0.2 m, which is remarkably smaller than the prediction error of the least squares extrapolated model and the dead reckoning model. Therefore, the forecasting function can achieve a high precision. The extrapolated trajectory of each epoch under the sinusoidal trajectory is simulated and shown in Figure 10. Simultaneously the details of the prediction trajectories of the different algorithms under partial data are amplified, the results are shown in Figure 11. Figure 11(b) is an enlarged section of the red frame in Figure 11(a) to highlight the details of the positioning effect.

Figure 10. The trajectory of each epoch.

Figure 11. Extrapolated trajectory comparison chart.

In Figure 11(b), positioning points P ₁, P ₂, P ₃, P ₄ are known. The extrapolated predictive trajectories of the generalised extension extrapolation, the least squares extrapolation and the dead reckoning extrapolation are given, respectively. The theoretical position trajectory at the time of simulation and the predicted position trajectory under the three methods are given, and a partial magnification is given in the lower right corner. It can be seen from the partial enlargement diagram that the predicted value of the generalised extension extrapolated model is the closest to the theoretical position, the accuracy is higher than the initial positioning point at the next time, and the extrapolated curve is smoother. Therefore, the feasibility of the algorithm has been verified by simulation.

3.2. Experimental application in real environment

A GNSS receiver carried on chip ATGM332D from Zhongke microelectronics was selected to test car navigation positioning in the Beijing Olympic Park of the Chinese Academy of Sciences (CAS) to verify the practical application effect. The hardware of the receiver is shown in Figure 12. The receiver's mean positioning accuracy was 2·5 m (in Wide-Open Field). The tracking sensitivity was −162 dBm. The cold-start capture sensitivity was −148 dBm and the output frequency was 1 Hz.

Figure 12. GNSS receiver hardware device.

The detailed testing procedures are as follows.

1. 1. Obtain some prior data through vehicular GNSS, such as P _t−1(x _t−1, y _t−1), P _t(x _t, y _t).

2. 2. Use the prior data to extrapolate the position at the next time $P^{\prime}_{t+1}\lpar x^{\prime}_{t+1}\comma \; y^{\prime}_{t+1}\rpar $ by generalised extension extrapolated model.

3. 3. At time t + 1, obtain the position P _t+1(x _t+1, y _t+1) at time t + 1 through GNSS, which can be taken as prior data so that the position $P^{\prime}_{t+2}\lpar x^{\prime}_{t+2}\comma \; y^{\prime}_{t+2}\rpar $ of the next time can be extrapolated.

4. 4. Obtain the location coordinate P _t+2(x _t+2, y _t+2) at time t + 2, which can be taken as prior data so that the position $P^{\prime}_{t+3}\lpar x^{\prime}_{t+3}\comma \; y^{\prime}_{t+3}\rpar $ of the next time can be extrapolated; repeat the above until time n.

The test result is shown in Figure 13. The original location points P _t−1, P _t, P _t+1, …, P _n are marked with blue “*” in Figure 13. The extrapolated trajectory points $P^{\prime}_{t-1}\comma \; P^{\prime}_{t}\comma \; \ldots\comma \; P^{\prime}_{n}$ are marked with the red curve in Figure 13. The velocity of the platform and course angle were provided by data from the $GPVTG message in GNSS positioning information.

Figure 13. Measured data extrapolated effect graph.

To further demonstrate the extrapolated effect of this method, Figure 14 gives a detailed view of the corner. The original location points are marked with blue “*”, the extrapolated trajectory points are marked with a red curve. From the enlargement diagram we can see that the method can also have good prediction effect when the vehicular vector angle is large and can provide real-time prediction of the moving trajectory at the next time with high extrapolation precision. Therefore, the practical feasibility of the algorithm is verified by the experiment in a real practical environment.

Figure 14. Measured data details.

Figure 15 is the accumulation curve of the positioning trajectory error from Figure 14. It can be seen that the confidence probability distributions of the positioning errors within 1 m, 1.5 m and 2 m are 65.3%, 89.8% and 96.4%. The average positioning accuracy is 0.72 m and the positioning error of 90% probability is less than 1.5 m. Therefore, the proposed algorithm can achieve high-precision trajectory prediction for land vehicle navigation applications.

Figure 15. Positioning error cumulative distribution curve.