Abstract
Aiming to improve the position and velocity precision of the INS/GNSS system during GNSS outages, a novel system that combines unscented Kalman filter (UKF) and nonlinear autoregressive neural networks with external inputs (NARX) is proposed. The NARX-based module is utilized to predict the measurement updates of UKF during GNSS outages. A new offline approach for selecting the optimal inputs of NARX networks is suggested and tested. This approach is based on mutual information (MI) theory for identifying the inputs that influence each of the outputs (the measurement updates of UKF) and lag-space estimation (LSE) for investigating the dependency of these outputs on the past values of the inputs and the outputs. The performance of the proposed system is verified experimentally using a real dataset. The comparison results indicate that the NARX-aided UKF outperforms other methods that use different input configurations for neural networks.
- Artificial neural networks
- Global Navigation Satellite Systems
- Inertial Navigation Systems
- NARX
- unscented Kalman filter
1 INTRODUCTION
In order to overcome the shortcomings associated with the stand-alone operation of Inertial Navigation Systems (INS) and Global Navigation Satellite Systems (GNSS), and to combine advantages of each system, INS and GNSS are often integrated to obtain accurate navigation solution with superior performance in comparison with either a GNSS or an INS stand-alone system. Many fusion algorithms are employed to fuse INS and GNSS data; the traditionally employed fusion techniques are Kalman filters (KF), such as extended Kalman filter (EKF) (Al Bitar & Gavrilov, 2019; Crassidis, 2006; Faruqi & Turner, 2000) and unscented Kalman filter (UKF) (Al Bitar & Gavrilov, 2019; Chang, 2014; Crassidis, 2006). With correct dynamic and stochastic models of GNSS and INS errors, KF can produce very accurate solutions, if there is continuous access to GNSS signals. However, KF does have limitations. The major inadequacy related to the utilization of KF for INS/GNSS integration is the necessity to have accurate stochastic models for each of the sensor errors. The inaccurate description of the system noises, measurement errors, and uncertainty in the dynamic models lead to unreliable estimates and degradation in accuracy, especially during GNSS outages when KF operates in prediction mode based on the predefined state error models, which are not necessarily correct. In addition, there are several significant drawbacks of KF, such as sensor dependency and observability problems (Hong et al., 2005; Klein & Diamant, 2018; Tang et al., 2008).
The limitations of KF have motivated researchers to investigate alternative methods for improving the accuracy of navigation solution during GNSS outages. These methods were predominantly based on artificial intelligence (AI). Much research has been conducted to investigate the use of AI-based techniques to bridge GNSS signal outages in INS/GNSS systems. Researchers have utilized various approaches for combining the AI module(s) with the rest of the INS/GNSS system (Al Bitar et al., 2020).
Chiang et al. (2008) suggested the replacement of KF by AI module using the so-called position update architecture (PUA). The proposed scheme was implemented using a constructive neural network (CNN) to overcome the limitations of conventional techniques that are predominantly based on the KF. The PUA module is used to estimate the INS position error during GNSS signal outages using velocity and azimuth of the INS.
Chiang and Huang (2008) proposed position, velocity and azimuth update architecture (PVAUA) based on multilayer perceptron neural networks (MLPNN). The PVAUA module uses the velocity, azimuth of the INS and time to estimate the INS position error during GNSS signal outages.
El-Sheimy et al. (2004) introduced an alternative INS/GPS integration method using an adaptive neuro-fuzzy inference system (ANFIS). The ANFIS-based module is implemented to predict the error drift of the standalone INS-estimated position during GNSS signal blockage using the INS position and time.
In fact, the replacement of KF by an AI module worked well for navigational-grade INS. However, these techniques showed a very limited success when applied to a MEMS-based INS/GNSS navigation system, due to the high noise level and bias instability of MEMS inertial sensors. As a result, KF is kept as the primary state estimation tool in INS/GNSS integration, and thus the logical step was toward an integration technique that uses both KF and AI module in the same system.
Wang et al. (2007) first proposed the concept of aiding KF. The authors utilized radial basis function neural networks (RBFNN) to predict the position differences between INS and GNSS in three orthogonal directions to form estimates of the measurement update of EKF during GNSS outages. The input parameters of the RBFNNs are the attitude angles and the changes of vehicle velocity and attitude angles in each epoch. However, only position measurements of KF are predicted during GNSS outages. This means that KF is not fully operational as no velocity measurements are predicted.
Chen and Fang (2014) proposed a hybrid prediction method for bridging GNSS outages using RBFNN and time series analysis, which aided EKF by forecasting position and velocity measurement updates. The proposed hybrid prediction method uses RBFNN to predict the six components of position and velocity measurement updates. The inputs of RBFNN are the measurements of accelerometers and gyroscopes. These measurements are first passed through a Wavelet denoising filter in order to lower the noise level. The residual error of training the RBFNN is modelled as a time series. The outputs of the RBFNN and the time series are summed together to form the final prediction of position and velocity measurement updates for EKF during GNSS outages. However, the complexity of the proposed system is not suitable for real-time implementation.
Jingsen et al. (2016) proposed a hybrid prediction method that combines extreme learning machines (ELM) and EKF. ELMs are applied to predict EKF position and velocity observations during GNSS outages. The measurements of gyroscopes and accelerometers are selected as inputs of the ELMs. The use of raw measurements of gyroscopes and accelerometers without a denoising stage complicates the learning process of ELMs, especially in the case of MEMS-based INS, as the level of noise is relatively high.
Yao et al. (2017) proposed a hybrid fusion algorithm to provide a pseudo position information to assist the integrated navigation system during GNSS outages using MLPNN. The proposed MLPNN-based model relates the current and past one-step values of velocities, angular rates and specific force of INS to the increments of the GNSS position. The GNSS position increments are accumulated to achieve the pseudo-GNSS position measurements.
Yao and Xu (2017) proposed robust least squares support vector machine (RLS-SVM)-aided fusion methodology for INS during GNSS outages. The RLS-SVM is used to predict the pseudo-GNSS position during GNSS outages similar to the previously proposed system in Al Bitar et al. (2020). The inputs of the RLS-SVM model are the specific force, velocity, and yaw information.
Wang et al. (2019) proposed a fusion algorithm based on back propagation neural networks (BPNNs) to predict the pseudo-GNSS position during GNSS outages. The proposed BPNN-based model relates the current and past values of velocities, angular rates, specific force of INS and the time elapsed since the beginning of GNSS signal outage to the increments of the GNSS position.
Recently, Fang et al. (2020) proposed an algorithm based on long short-term memory (LSTM) to predict the pseudo-GNSS position during GNSS outages. The inputs of LSTM model are the four-step information of specific force, angular rates, velocity and yaw.
There are some common drawbacks related to the methods mentioned above. It is clear that the selection of inputs of an AI module differs from one method to another without any justification or comparison. In fact, the selection of the inputs of an AI module affects the system. A fewer number of inputs means a simpler internal structure of AI module and, consequently, less training time, while a large number means more complicated structure and thus longer learning time and less real-time capability. Including a wrong input or excluding a right input lead to degradation in prediction accuracy of the AI module. The selection of the measurements of gyroscopes and accelerometers as inputs of an AI module is usually justified by the fact that the outputs of INS (position, velocity and attitude angles) are the result of the mathematical integrating (over time) of these measurements. In other words, the measurement errors (noises) of gyroscopes and accelerometers are embedded in the errors of INS outputs. As a result, the erroneous INS outputs (during GNSS outages) can be used as inputs of an AI module instead of the raw measurements of gyroscopes and accelerometers. One advantage of using erroneous INS outputs instead of the measurements of gyroscopes and accelerometers is that they do not require a denoising stage, as the process of mathematical integration smooths them.
The second drawback of all aforementioned methods is using EKF as the only choice for fusing INS and GNSS data. Other options, such as UKF, for example – which is proven to be less sensitive to the nonlinearities of process and observations models compared to EKF – were not considered.
The third drawback is that these methods use the same covariance matrices of KF during the availability and the outages of GNSS. This is not true, as the measurements provided by the aid of the AI module during GNSS signal have error characteristics that differ from those of GNSS measurements.
This paper considers the problem of aiding KF in INS/GNSS systems using a nonlinear autoregressive neural network with external inputs (NARX) (Siegelmann et al., 1997). The paper also addresses the problems mentioned above. First, UKF is chosen as the integrating filter instead of EKF. Secondly, a new approach for selecting the optimal inputs of an AI module (NARX networks, in our case) is proposed. This approach is based on mutual information (MI) theory (Brown, 2009; Peng et al., 2005) for identifying the inputs that influence each of the outputs (the measurement update for UKF during GNSS outages), and lag-space estimation (LSE) (He & Asada, 1993) for determining the model order (i.e., the dependency of the outputs on the past values of inputs and the outputs themselves). Third, the covariance matrices of UKF are linked to the prediction errors of AI module. The proposed method is shortly called “NARX-aided UKF.” The performance of the NARX-aided UKF is experimentally verified using a real dataset.
The rest of this paper is organized as follows: A detailed explanation of the proposed method is given in Section 2. The performance of the proposed method is presented in Section 3. Conclusions are presented in Section 4.
2 STRUCTURE OF THE PROPOSED SYSTEM
2.1 Principle of operation
The main idea of the proposed system is to employ a NARX network to predict the measurement update (the difference between GNSS outputs and the outputs of INS) of UKF during GNSS outages. The system works in two modes: learning mode when a GNSS signal is available, and prediction during GNSS outages. Two copies of INS are created, INS1 and INS2. In learning mode, INS1 and GNSS are integrated using a loosely coupled scheme. The positions and velocities P𝐺𝑁𝑆𝑆, 𝑉𝐺𝑁𝑆𝑆 provided by GNSS are merged as updates of the INS1 estimates of position and velocity P𝐼𝑁𝑆1,𝑉𝐼𝑁𝑆1 through a UKF (Figure 1). UKF estimates the errors on state of INS1 
. These estimates are added to state of INS1 X𝐼𝑁𝑆1 to form corrected values X𝐶. At the same time, INS2 works autonomously, and its outputs P𝐼𝑁𝑆2, 𝑉𝐼𝑁𝑆2 are corrected periodically every 60 seconds by GNSS measurements, as shown by the dashed arrows. In fact, the 60-s duration is related to real-life scenarios; in real life, the GNSS signal may be lost when moving through tunnels or around obstacles in urban areas. The duration of these outages in most cases is less than 60 s. The position and velocity of INS2 are then subtracted from the ones of the GNSS to form error signals of position 𝛿𝑃𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2 and velocity 𝛿𝑉𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2. These error signals are used as target values for training NARX networks. The estimates of gyro drifts and accelerometer biases are fed back to INS1 and INS2 mechanization equations, as shown by the dotted arrows.
Principle of operation of proposed method for aiding UKF during GNSS signal outages - training mode
The training algorithm of NARX networks starts after collecting a specific amount of training data that represents the 60-s simulated outage duration. This amount is called a window (δPGNSS/INS2,𝛿𝑉GNSS/INS2. The inputs of NARX networks are chosen from the outputs of INS1 (position, velocity and attitude angles), the past values of errors, and the time elapsed since the beginning of GNSS signal outage. The process of selecting the inputs of NARX networks is based on MI theory and LSE, and it is conducted in offline stage, as will be explained later in Subsection 2.5.
When real GNSS signal outage occurs, the system switches to prediction mode, as shown in Figure 2. The NARX module predicts the errors 
. These errors are added to the position and velocity of INS2 P𝐼𝑁𝑆2,𝑉𝐼𝑁𝑆2 to form estimations of position and velocity of GNSS 
. The difference between 
and P𝐼𝑁𝑆1,𝑉𝐼𝑁𝑆1 forms the estimation of measurement update for UKF 
. Using these updates, UKF continues to operate as if no GNSS outage had occurred.
Principle of operation of proposed method for aiding UKF during GNSS signal outages - prediction mode
To train the NARX networks in online mode, a training procedure that utilizes non-overlap moving window technique is used. The non-overlap moving window doesn’t have the disadvantage of redundancy in the information when using sliding window; thus, it doesn’t require a long time for data processing compared to the sliding window technique. The NARX networks are updated (trained) within this window. For real-time purposes, the NARX networks are trained until reaching certain minimum mean-squared error (MSE) or after completing a certain number of training epochs (determined empirically). This procedure is repeated when a new window is acquired, as shown in Figure 3. Whenever a GNSS outage occurs, the NARX networks switch to prediction mode and provide estimates of errors (
).
Online training of NARX networks
2.2 Navigation equations
Taking into consideration the coordinate frames shown in Figure 4, the navigation equations written in N-frame are given as follows (Jekeli, 2012):
Coordinate frames: 𝑂𝑋𝐼𝑌𝐼𝑍𝐼 is the geocentric inertial coordinate system (I-frame); 𝑂𝑋𝐸𝑌𝐸𝑍𝐸 is the geocentric Earth coordinate system (E-frame); 𝑂𝑁𝐸𝐷 is the Local North-East-Down coordinate system (N-frame); 𝑂𝑋𝑌𝑍 is the body fixed coordinate system (B-frame); 𝑎 is the semi-major axis of the ellipsoid, 𝑏 is the semi-minor axis of the ellipsoid, 𝜔𝐸 is the Earth’s angular rate, and 𝑡 is the time
1
2
where 𝑷 = [𝜑 𝜆 ℎ]𝑇 is the object’s position and 𝜑, 𝜆, ℎ are the latitude, longitude, and height of the object’s center of mass. 
is the object’s velocity relative to the E-frame written in N-frame. r𝐸 is the object’s center-of-mass coordinate vector in e-frame, and 
is the matrix of direction cosines from E-frame to N-frame. 
is the vector of specific force in N-frame, where 𝒇𝐵 is the vector of specific force in B-frame (the output signals of the accelerometer triad) and 
is the matrix of direction cosines from B-frame to N-frame. 
are skew-symmetric matrices composed of angular velocities 
, where 
is the angular rate vector of E-frame relative to I-frame written in N-frame, and 
is the angular rate vector of N-frame relative to I-frame written in N-frame. 
is a diagonal matrix 
, where M and 𝑁 are radii of curvature of ellipsoid (the figure of the Earth is described by biaxial ellipsoid (Jekeli, 2012). 
is the gravity vector written in N-frame, where 𝒈𝑁 = [ 0 0 𝑔 ]𝑇 is gravitational vector, and g is given according to World Geodetic System WGS-84 and 
, where 
is the angular rate vector of E-frame relative to I-frame written in E-frame.
The matrix 
can be represented through the Rodrigues–Hamilton parameters (quaternions) 
. The time behavior of quaternion is described by the following differential equation:
3
where
4
where 
is the angular rate vector of the B-frame relative to I-frame (the output signals of the gyroscopes triad). The attitude angles (Euler angles) are calculated from quaternion components using the following equations:
5
where 𝜃, 𝜙, 𝜓 are pitch, roll and yaw angles.
2.3 Measurement model of inertial sensors
The inaccurate measurements of inertial sensors are explained by various reasons, among them, nonorthogonality of the measuring axes of the units of accelerometers and gyroscopes, and biases that can be represented as a sum of systematic and random components. The measurement model of accelerometers and gyroscopes based on MEMS technology can be written in a generalized form as (El-Sheimy et al., 2007; Jafari et al., 2014; Quinchia et al., 2013)
6
7
8
9
10
11
12
13
where 
are the three-dimensional vectors of the output signals of the accelerometers and the gyroscopes, respectively, and 
are their true values. 
are 3×3 coefficient matrices. 𝑰3 is 3× 3 unity matrix. 𝑩𝐴,𝑆,𝑩𝐺,𝑆 are the systematic components of the accelerometer biases and gyro drifts, respectively. 𝑩𝐴𝐶𝐶𝑅𝑊,𝑩𝑅𝑅𝑊 are the acceleration and the rate random walks, respectively, and 𝑾𝐴𝐶𝐶𝑅𝑊,𝑾𝑅𝑅𝑊 are zero-mean white noises. 𝑩𝐴,𝐺𝑀,𝑩𝐺,𝐺𝑀 are first order Gauss-Markov (GM) processes, 𝑻𝐴,𝐶,𝑻𝐺,𝐶 are 3×3 correlation matrices, and 𝑾𝐴,𝐺𝑀,𝑾𝐺,𝐺𝑀are zero-mean white noises. 𝑾𝑉𝑅𝑊,𝑾𝐴𝑅𝑊 are zero-mean white noises that represent the velocity and the angle random walks, respectively, and 𝑺𝐴,𝑺𝐺 are coefficient matrices that include scale factors and other coefficients due to the non-orthogonality of measuring axes of the accelerometers and the gyroscopes blocks.
2.4 Unscented Kalman filter
The UKF uses a deterministic sampling technique known as the unscented transformation to pick a minimal set of 2L+1 sample points (called sigma points) around the mean, where L is the size of state vector. The sigma points are then propagated through the non-linear functions, from which a new mean and covariance estimate are then formed (Al Bitar & Gavrilov, 2019; Crassidis, 2006). In addition, the UKF removes the requirement to calculate the Jacobians, which can be a difficult task for complex functions. Compared to EKF, UKF is less sensitive to the nonlinearities of process and observation models. In order to apply the algorithm of UKF, it is necessary to write the process and measurement equation of INS/GNSS in discrete time. The process equation in discrete time can be written as
14
where 𝑿𝑘 is the state-vector of INS and 𝑾𝑘 is the vector of process noise:
15
𝒇(𝑿𝑘,𝑾𝑘) is nonlinear vector-function that can be written by transforming Equations (1)–(3) and Equations (8)–(13) into discrete time
16
17
where 𝑇𝑆 is the sampling time of INS.
The measurement equation of UKF is given as
18
where 𝒁𝑘 is the vector of GNSS position and velocity measurements:
19
where 𝝊𝑘 is measurement noise, which is assumed to be zero-mean white noise with covariance matrix 𝑹𝑘. The measurement covariance matrix 𝑹𝑘 is given by
20
where 
, 
are the standard deviations of GNSS position and velocity measurements’ noise, respectively.
The UKF can be represented as a “prediction– correction” procedure. To initialize the UKF, weight coefficients are calculated for each of the 2𝐿 + 1 sigma-points according to the following rules:
21
where 𝛼, 𝛽, 𝜅 are parameters that determine the position of sigma points in state-space. 𝛼, 𝜅 regulate the spread of points relative to their expected value. Parameter 𝛽 is used to incorporate prior knowledge of the distribution. For a normal distribution, it is conventional to set the values of these parameters: 𝜅 = 0, 𝛽 = 2, 𝛼 ∈ [10−4,1] (Crassidis, 2006).
At the prediction stage, 2𝐿 + 1 sigma-points are generated. These sigma-points form the 2𝐿 + 1 columns of matrix 𝑆𝑘−1 as follows:
22
where 
are the posterior estimates of statevector and covariance matrix, respectively. 𝑸𝑘 is process noise covariance matrix. The Cholesky decomposition is an effective method for calculating the square root of 
. Further, to denote a column of 𝑺𝑘 we use the index 𝑗 = 0, 1, …, 2𝐿. For example, 𝑺𝑘,𝑗 means a j-th column of matrix 𝑺𝑘. The columns of 𝑺𝑘 are propagated through nonlinear function 𝒇(∗)
23
Then, the a priori estimates of state-vector and covariance matrix are calculated using weighted sums
24
25
At the correction stage, the a priori estimates of statevector and covariance matrix are updated (corrected)
26
27
28
where
29
30
31
32
The initial value of covariance matrix 
is given by
33
where 
denotes the accuracy of the initial alignment of the INS,
,
denote the initial errors in determining constant gyro drifts and accelerometer biases, 𝝈𝑅𝑅𝑊,𝝈𝐴𝐶𝐶𝑅𝑊 are the standard deviations of the noises of the rate and the acceleration random walks, and 𝝈𝐺,𝐺𝑀,𝝈𝐴,𝐺𝑀 are the standard deviations of the noises of GM processes.
2.5 Offline stage
Four essential tasks are performed in offline stage: 1) the selection of the optimal inputs of NARX networks, 2) the design of the internal structure (the number of layers/neurons) of NARX networks, 3) preliminary training of NARX networks, and 4) calculating the new covariance matrices of UKF that will be used during GNSS outages (in online mode).
2.5.1 Offline data
A dataset is acquired from both INS and GNSS during a trip that contains as many maneuvers as possible. The data of INS1 and GNSS are fused by UKF using the loosely coupled scheme, as shown in Figure 5.
The calculation of target values and candidate inputs for NARX networks
The target values (the measurement updates) 𝛿𝑃𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2,𝛿𝑉𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2 are obtained as follows:
34
The target values will have the shape of the signal shown in Figure 5. The INS2 outputs are position 𝑷𝐼𝑁𝑆2, velocity 𝑽𝐼𝑁𝑆2 and attitude angles 𝑨𝐼𝑁𝑆2
35
The goal is to choose a set of inputs from the group 𝐼𝐼𝑁𝑆2 = {𝜑,𝜆,ℎ,𝑣𝑁,𝑣𝐸,𝑣𝐷,𝜙,𝜃,𝜓,𝑡} that have an impact on each of the six error components from the group 𝑇 = {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷}, where t is the time elapsed since the loss of the GNSS signal and varies in interval [0, 60 s]. The measurement updates {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} are predicted using six NARX networks. The inputs of NARX networks are selected based on MI criterion and the LSE. First, the MI criterion ranks the inputs in group 𝐼𝐼𝑁𝑆2, where the inputs with positive rank are selected. Then, LSE is used to determine the model orders, that is, the dependency of each error component in group 𝑂𝑂𝑈𝑇 = {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} on the past values of the error itself and the past values of the inputs chosen by MI.
Next, the theoretical background of MI and LSE is presented (subsections 2.5.2 and 2.5.3). The structure of NARX network is presented in Subsection 2.6.
2.5.2 Mutual information
The relation between the inputs 𝐼𝐼𝑁𝑆2 and outputs 𝑇 is not linear, so the methods based on linear relations (like correlation) are prone to mistakes and will not give accurate results. The MI criterion is a good candidate to solve this problem, as it measures the arbitrary (linear or nonlinear) dependencies between variables. MI is widely used in machine learning for canonical tasks, such as classification, clustering and feature selection (Brown, 2009; Peng et al., 2005). MI is one of the feature selection methods. These methods define a statistical criterion that is used to rank characteristics (or features) according to their usefulness for classification. The characteristics with a high rating are chosen, and characteristics with a low rating can be discarded. Given two random variables 𝑥 and 𝑦, their mutual information is determined through probability density functions 𝑝(𝑥), 𝑝(𝑦), 𝑝(𝑥, 𝑦)
36
where 𝐷𝑥,𝐷𝑦 are the spaces of (𝑥, 𝑦). The concept of MI can be employed for arranging or ranking a number of candidates (features) 𝑥𝑙,𝑙 = 1,…,𝑀 according to their usefulness or influence on target signal 𝑦, based on 𝑁 samples of 𝑥𝑙 and 𝑦. Peng et al. (2005) proposed a solution for this problem by calculating a rank based on the Maximum Relevance Minimum Redundancy (MRMR) criterion. The rank is calculated as follows:
37
The MRMR criterion takes into account the redundancy of candidates 
and subtracts it from 𝐼(𝑥𝑙;𝑦). The MRMR criterion can be summarized as “a set of features should not only be individually relevant, but also should not be redundant in relation to each other.” By calculating 𝐽𝑀𝑅𝑀𝑅(𝑥𝑙), the features can be arranged or ranked; the feature with the largest 𝐽𝑀𝑅𝑀𝑅 has the largest effect or influence on 𝑦 and vice versa. The MRMR criterion is applied to rank the candidate inputs in group 𝐼𝐼𝑁𝑆2 in accordance to their impact on each of the six error components from the group 𝑇. The results of applying MI criterion are provided in Subsection 3.2.
2.5.3 Lag-space estimation
The next task is to investigate the dependency of each error component in group 𝑇 on the past values of the error itself and the past values of the inputs that were chosen by MRMR criterion. This problem is referred to as lag-space estimation, or model order estimation. The use of higher order dependencies in the modeling process leads to possibly over-parameterized, and thus less efficient, models. It is therefore important to estimate the optimal lag-space, that is, to find the primary dependencies. This allows minimizing the number of parameters and optimizing the predictive abilities of the AI module (NARX module, in our case). In system theory, a nonlinear dynamical system is generally described by differential or difference equations that represent input/output relations. However, in many practical situations, it is difficult to write down the accurate state dynamics and observation equations for a continuous or a discrete time system. What are available are the input and output data of the unknown dynamical system, that is, 𝑢(𝑡) and 𝑦(𝑡), which are observed at sampling times 𝑡𝑖 = 𝑖𝑇𝑆,𝑖 = 0,…,𝑁 − 1. It has been shown that under some mild assumptions, the following input/output model (Siegelmann et al., 1997):
38
can represent nonlinear dynamical systems described by differential or difference equations, where 𝑔(∗) is nonlinear function, and parameters 𝑛𝑦 and 𝑛𝑢 are orders of the input/output model that should be determined. In our case, 𝑢(𝑡), is a subset of the group 𝐼𝐼𝑁𝑆2 and 𝑦(𝑡) is one of the six error components in group 𝑇. Authors He and Asada (1993) proposed a method to identify model orders 𝑛𝑦,𝑛𝑢 based on Lipschitz quotients. According to them, the model described by Equation 38 can be written as
39
where 
are the input variables and n is the number of input variables 𝑛 = 𝑛𝑦 +𝑛𝑢 +1. Now, the goal is to reconstruct the function 𝑔(∗) based on the pairs (
). The Lipschitz quotient for 𝑚 input variables can be calculated by
40
Usually, the following index is used to determine the optimal amount of input variables
41
where 𝑞(𝑚)(𝑘) is the k-th largest Lipschitz quotient among all 
with m input variables, and p is a positive parameter (𝑝 = 0.01𝑁 ∼ 0.02𝑁). If n is the optimal number of input variables, then 𝑞(𝑛+1) is very close to 𝑞(𝑛) and 𝑞(𝑛−1) is much larger than 𝑞(𝑛). Moreover, 𝑞(𝑛−2) is much larger than 𝑞(𝑛−1) and 𝑞(𝑛+2) is very close to 𝑞(𝑛+1). Therefore, looking at the curve of 𝑞(𝑚) as a function of 𝑚, we can observe that starting from a certain value 𝑚 = 𝑛, further increase in m will not significantly change the index 𝑞(𝑚) and thus the value of n can be determined (Figure 6). The results of applying this LSE method are provided in Subsection 3.2.
The typical curve of 𝑞(𝑚)
2.5.4 The reconfiguration of the covariance matrices of UKF
During GNSS outages, the proposed system works in prediction mode. The measurements obtained by the aid of NARX networks have error characteristics that differ from the characteristics of GNSS measurements. Therefore, it is necessary to reconfigure the UKF covariance matrices. To do this, the following steps are performed: First, the proposed system is applied using the offline dataset for simulated GNSS outages. Then, the standard deviation of the position and velocity errors with respect to GNSS measurements is calculated as follows:
42
where 
,
are the position and the velocity obtained by applying the proposed system, 𝝁𝑃,𝝁𝑉 are the mean values of the errors, and 𝑁𝑝 is the number of data samples. When a real GNSS outage occurs, the corresponding components of the covariance matrices 𝑷𝐶𝑂𝑉,𝑹 are updated with 
,
.
2.6 NARX neural network
Considering the system model given by Equation (38), a good choice of AI module is NARX, as it obeys the same system equation. The NARX is a recurrent dynamic neural network that can be used to model extensive variety of nonlinear dynamic systems. NARX networks have been applied in various applications, including black-box system identification and time-series modeling (Diaconescu, 2008; Siegelmann et al., 1997). The architecture of the NARX is shown in Figure 7.
The architecture of NARX
Six NARX networks are utilized to predict systems errors 𝛿𝑃𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2 = [ 𝛿𝜑 𝛿𝜆 𝛿ℎ ]𝑇, 𝛿𝑉𝐺𝑁𝑆𝑆/𝐼𝑁𝑆2 = [ 𝛿𝑣𝑁 𝛿𝑣𝐸 𝛿𝑣𝐷 ]𝑇 during GNSS outages.
3 RESULTS
3.1 Experimental setup
Raw experimental data were acquired from a Micro-Electro Mechanical System-Strapdown Inertial Navigation System (MEMS-SINS) (Ekinox-D Inertial Navigation System) with sampling frequency 200 Hz. The characteristics of gyroscopes and accelerometers of this SINS were obtained in (Gonzalez & Dabove, 2019; Gonzalez et al., 2017) using Allan variance method. A Global Navigation Satellite System/Global Positioning System (GLONASS/GPS) receiver was used with sampling frequency of 5 Hz. The accuracy in position is 0.5 m for latitude and longitude and 1 m for altitude. The accuracy in velocity for all components is 0.1 m/s. Both systems were mounted on the roof of the vehicle, as shown in Figure 8a. The experiment was conducted in the city of Turin in Italy (Gonzalez et al., 2017). The duration of the dataset used in this work is 2,300 s. The trajectory of the vehicle is shown in Figure 8b. The first segment of the trajectory (the first 600 s) is magnified. The dataset of the first segment is used for offline stage. The rest of the dataset (1,700 s) is used for online validation of the proposed system. Table 1 shows the specifications of the gyroscopes and the accelerometers of the Ekinox-D INS (Gonzalez & Dabove, 2019).
(a) Experiment hardware; (b) trajectory
The specifications of the gyroscopes and the accelerometers of the Ekinox-D INS
3.2 The results of the offline stage
The dataset of the first segment of the trajectory is used in offline stage. Figure 9a shows the first segment of the trajectory, and Figure 9b shows the speed of the vehicle ||𝑽|| along this segment. For better readability, the speed of the vehicle is shown in (km/h). The segment contains many types of possible maneuvers (accelerating and decelerating, zero velocity, straight lines, turning, etc.).
(a) The first segment of the trajectory of the vehicle; (b) the speed of the vehicle along the first segment of the trajectory
Using the dataset of the first segment, six GNSS outages (six windows) were simulated to form the target signals {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} as shown in Figure 10.
Target signals for offline preprocessing stage
Figure 11 shows the candidate input signals {𝜑,𝜆,ℎ,𝑣𝑁,𝑣𝐸,𝑣𝐷,𝜙,𝜃,𝜓} along the first segment.
The candidate input signals {𝜑,𝜆,ℎ,𝑣𝑁,𝑣𝐸,𝑣𝐷,𝜙,𝜃,𝜓} along the first segment of the trajectory
The first task is to identify the inputs that influence each of the targets {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} using the MRMR criterion explained earlier. The results of applying MRMR criterion are presented in Table 2. The large negative score means that the input has a high redundancy, that is, the information imbedded in this input are found in other inputs, so there is no need to consider this input. The positive score (in bold) means high relevance and low redundancy of the input. The scores close to zero (underlined) reflect an insignificant influence of the inputs on the corresponding output and can be neglected. As a result, the inputs with high positive scores are chosen. It is worth mentioning that the results given in Table 2 are limited to land vehicles, and they cannot be generalized for the case of aerial vehicles because they have different types of movement.
The results of applying MRMR algorithm
The second task is to determine the dependency of the target signals {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} on their past values and the past values of selected inputs, that is, the inputs selected using MRMR algorithm and shown in bold in Table 2, using LSE explained earlier. Figure 12 shows the results of applying LSE to investigate the dependency of the target signals {𝛿𝜑, 𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} on their past values. To determine the proper lag-space, we look at the point where increasing the lag-space (𝑚) will not change the order index 𝑞(𝑚) significantly. As Figure 12 shows, a lag-space of 𝑛𝑦 = 2 is a good choice for all target signals, as the order index doesn’t change significantly for values larger than 2.
Lag-space estimation for the position and velocity errors
Next, the LSE is applied to determine the lag-space for inputs selected using MRMR algorithm. Figure 13 shows the results of applying LSE for the case of target signal 𝛿𝜑 and the selected inputs {𝜑, 𝑣𝑁}. As Figure 13 shows, the order index doesn’t change significantly for values larger than 2. This means that the lag-space is 𝑛𝑢 = 2. For the other cases {𝛿𝜆, 𝛿ℎ, 𝛿𝑣𝑁,𝛿𝑣𝐸,𝛿𝑣𝐷} the same lag-space was found (𝑛𝑢 = 2).
Lag-space estimation for the case of target signal 𝛿𝜑 and the selected inputs {𝜑, 𝑣𝑁}
At this stage, all inputs of the six NARX networks are determined.
Table 3 summarizes the results of applying LSE. The symbol (-) means that there is no dependency found between the output and the input.
The results of applying LSE method
Figure 14 shows the final input/output configurations of the six NARX networks according to Tables 2 and 3.
Final input/output configurations of the six NARX networks
The next task is the design of the internal structure of the NARX networks, that is, the number of hidden layers and neurons. It is proved that an artificial neural network (ANN) with two hidden layers can approximate any function (Gonzalez & Dabove, 2019; Lippmann, 1987). Therefore, six ANN with two hidden layers were utilized. Choosing the right number of neurons is very important. An ANN with a small number of neurons will not be able to learn. A large number of neurons will lead to an increase in the network training time and can also lead to over-fitting. In this article, the number of neurons is derived empirically using rules-of-thumb (Goodfellow et al., 2016; Huang, 2003). It is shown that an ANN with two hidden layers and 
neurons can learn 𝑁𝑆 example with any arbitrary precision, where 𝑁𝑦 = 1 is the number outputs. Here 𝑁𝑆 represents the number of samples in the window 𝑁𝑆 = 𝑊 samples. The selection of this value of window size is based on offline trials. In fact, there is a trade-off in choosing the window size; large window sizes guarantee that more motion dynamics are mimicked, thus providing better accuracy over long GNSS outages. On the other hand, a small window size guarantees fast learning, but the system provides high accuracy of estimation only for short GNSS outages. The window size is chosen as (𝑊 = 300 samples), which represents 60 s record of data, because the sampling frequency of GNSS is 5Hz. Therefore, the number of neurons is 𝑁 = 50. In fact, 𝑁 = 50 is only a starting value for the number of neurons. The final values were 30, 36, 20, 32, 36 and 40 for NARX-1, NARX-2, NARX-3, NARX-4, NARX-5 and NARX-6. These values were achieved after many offline trials. The hyperbolic tangent sigmoid (tan-sigmoid) transfer function is applied as activation function. To train the NARX networks, Levenberg-Marquardt (LM) (Moré, 1978) training algorithm is used. LM algorithm is the most widely used optimization algorithm. It outperforms other methods in a wide variety of problems because of its fast and stable convergence. The LM algorithm is well suited for training small and medium-sized problems, which are the case here. Figure 15 shows the results of preliminary training of the NARX networks based on the offline dataset.
The results of the preliminary training of NARX networks based on offline dataset
3.3 Online validation of the proposed system
In order to test the proposed system, the second part of the dataset (with duration of 1,700 s) is used. Six GNSS outages were simulated, as shown in Figure 16. The outages contain straight lines (outages 1, 2 and 3) and turnings (outages 4, 5 and 6).
The trajectory of vehicle with GNSS outages for online validation
Figure 17 shows the speed of the vehicle ||𝑽|| along the trajectory with the GNSS outages highlighted. It can be seen that the GNSS outage segments contain accelerating and decelerating (all outages), zero velocity (outages 1, 2 and 3), high speed ∼ 50-90 km/h (outages 2, 4, 5 and 6), low speed ∼ 0-20 km/h (outages 1, 2 and 3), and mid speed ∼ 20-50 km/h (outages 1 and 2).
The speed of the vehicle along the trajectory, with the GNSS outages highlighted
Figure 18 shows the GNSS-outage segments of the vehicle’s trajectory.
Performance of different methods (UKF, M1, M2 and NARX-aided UKF) for six GNSS outages
The proposed method (NARX-aided UKF) is compared to UKF and two widely adopted methods that use different input configurations in order to validate the selection of the input configurations of NARX networks. The first method (shortly M1) uses the current information of specific force and angular rates 
for estimating the position and velocity errors as in (Chen & Fang, 2014; Jingsen et al., 2016). The second method (shortly M2) uses the four-step information of specific force, angular rates, velocity and yaw 
as in (Fang et al., 2020). For comparison, the reference GNSS trajectory is also shown in Figure18. The superiority of the proposed method over UKF, M1 and M2 in terms of positioning accuracy is obvious.
Figure 19 shows the horizontal error (with respect to GNSS) in position using UKF, M1, M2 and NARX-aided UKF.
Horizontal error in position using UKF, M1, M2 and NARX-aided UKF for six GNSS outages
Table 4 provides numerical values of horizontal errors using different methods. The proposed method (NARX-aided UKF) improved the positioning accuracy by 82%, 65% and 55% with respect to UKF, M1 and M2 respectively.
Numerical values of horizontal errors (in meters) using UKF, M1, M2 and NARX-aided UKF
To demonstrate the effect of modifying the covariance matrices of UKF during GNSS outages, Figure 20 shows the horizontal errors in position in two cases: the first case, the NARX-aided UKF without modification of covariance matrices, and the second case with modification of covariance matrices. It can be seen that the updating of covariance matrices slightly enhanced the positioning accuracy. The improvement of positioning accuracy is 5% to 8%.
Demonstration of the effect of updating covariance matrices of UKF during GNSS outages
4 CONCLUSIONS
The problem of aiding UKF during GNSS outages in INS/GNSS systems was considered in this paper. A new method, namely NARX-aided UKF, was suggested and tested. The proposed method consists of offline and online stages. The offline stage is essential for selecting the input signals of NARX networks based on MI and LSE, the design of the internal structure of NARX networks, preliminary training of NARX networks and calculating the new covariance matrices of UKF that were used during real GNSS outages. The covariance matrices of UKF during GNSS outages were linked to prediction accuracy of NARX networks. The performance of the proposed method was experimentally verified using real datasets. The results indicated that the proposed method improved the accuracy of navigation systems during GNSS outages. The results also confirmed the superiority of the proposed method over widely adopted methods that use different input configurations for neural networks. The future work will consider the case of aerial vehicles.
HOW TO CITE THIS ARTICLE
Al Bitar N, Gavrilov A. A novel approach for aiding unscented Kalman filter for bridging GNSS outages in integrated navigation systems. NAVIGATION. 2021;68:521–539. https://doi.org/10.1002/navi.435
ACKNOWLEDGMENTS
The authors would like to thank professors Rodrigo Gonzalez (National University of Technology, Mendoza, Argentina) and Paolo Dabove (Politecnico di Torino University, Turin, Italy) for providing the experimental data.
- Received August 27, 2020.
- Revision received March 14, 2021.
- Accepted March 24, 2021.
- © 2021 Institute of Navigation
This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Jump to section
Related Articles
Cited By...
- No citing articles found.