Adaptive Filtering on GPS-Aided MEMS-IMU for Optimal Estimation of Ground Vehicle Trajectory

Fusion of the Global Positioning System (GPS) and Inertial Navigation System (INS) for navigation of ground vehicles is an extensively researched topic for military and civilian applications. Micro-electro-mechanical-systems-based inertial measurement units (MEMS-IMU) are being widely used in numerous commercial applications due to their low cost; however, they are characterized by relatively poor accuracy when compared with more expensive counterparts. With a sudden boom in research and development of autonomous navigation technology for consumer vehicles, the need to enhance estimation accuracy and reliability has become critical, while aiming to deliver a cost-effective solution. Optimal fusion of commercially available, low-cost MEMS-IMU and the GPS may provide one such solution. Different variants of the Kalman filter have been proposed and implemented for integration of the GPS and the INS. This paper proposes a framework for the fusion of adaptive Kalman filters, based on Sage-Husa and strong tracking filtering algorithms, implemented on MEMS-IMU and the GPS for the case of a ground vehicle. The error models of the inertial sensors have also been implemented to achieve reliable and accurate estimations. Simulations have been carried out on actual navigation data from a test vehicle. Measurements were obtained using commercially available GPS receiver and MEMS-IMU. The solution was shown to enhance navigation accuracy when compared to conventional Kalman filter.


Introduction
The GPS (Global Positioning System) and the INS (Inertial Navigation System) can provide stand-alone navigation solutions; however, both are prone to intrinsic errors/biases which may limit the accuracy and productivity of each navigation technique. The GPS provides fairly precise positioning information but at a low frequency of 1 Hz. It may also suffer from phenomena such as unavailability of a clear line-of-sight to the satellites, multipath propagation/interference in dense urban environments, signal jamming, and degradation in accuracy or complete denial to service by the US government [1,2]. INS, on the other hand, finds preference being a dead reckoning system

Inertial Navigation Mechanism
Inertial navigation requires determining a vehicle's attitude and position in a three-dimensional frame of reference. This is accomplished by using three perpendicularly mounted gyroscopes to provide rate-of-turn in three different axes, while three orthogonally mounted accelerometers provide accelerations in these three axes as experienced by the vehicle. Measurements provided by the gyroscopes are with respect to the inertial frame of reference. However, accelerometers cannot provide the acceleration with respect to the inertial frame only. Their measurements also include a component of acceleration caused by the gravitational force. Measurement provided by an accelerometer can be written as difference between the true acceleration of a body in space and the acceleration due to gravitational pull; that is, This quantity is known as the specific force, and can be described as the non-gravitational force exerted on per unit mass of the body. The equation can be rearranged to calculate the true acceleration. This equation is termed the navigation equation in literature (and can take different forms depending upon the choice of reference frame). If r is the unit vector representing the vehicle's position from origin, then acceleration of the vehicle in an inertial frame is given by: However, when the navigation process is being carried out in a non-fixed rotating frame of reference, such as the Earth frame, additional forces act on the vehicle as a function of the rotation of the reference frame. This effect is known as the Coriolis effect, and the velocity of vehicle can then be calculated using the inertial velocity as follows: where ω ie = [ 0 0 Ω e ] T represents the Earth's rate of turn with respect to the inertial frame, and the symbol × indicates vector product. The symbol Ω e = 15.041067°/hr is the Earth's angular speed.
The specific force measurements provided by accelerometers are expressed in the vehicle's body frame as: Specific force measurements from body frame represented by the subscript, b, are then transformed to the local geodetic or the navigation frame using the transformation matrix, C n b : where the subscript b and superscript n of the transformation matrix, C n b , denote body and frame references, respectively, and the transformation matrix is given by, cos ψ cos ϕ − sin ψ sin θ sin ϕ − sin ψ cos θ sin ψ cos ϕ + cos ψ sin θ sin ϕ cos ψ cos θ − cos θ sin ϕ sin θ cos ψ sin ϕ + sin ψ sin θ cos ϕ sin ψ sin ϕ − cos ψ sin θ cos ϕ cos θ cos ϕ where ϕ, θ, and ψ are vehicle's roll, pitch, and yaw angles, respectively. However, the three acceleration components expressed in the navigation frame cannot directly provide velocity components in the navigation frame. This is due to three reasons: (1) The Earth's rotation rate, which is expressed as angular velocity vector ω n ie in the local navigation frame. (2) There is a change in orientation of the local navigation frame with respect to the Earth's frame, which is because of the definition of the local level north and vertical directions. The local north is at a tangent to the Earth's meridian at all times, whereas the vertical direction is perpendicular to the Earth's surface. This effect can be expressed as the angular velocity vector ω n en , which is the turn-rate of the navigation frame with respect to the Earth-fixed frame, and is known as the transport rate. (3) the third reason is gravity of the Earth, expressed in terms of the local level gravity vector, g n l . Mathematically, these three factors can be expressed as where L represents vehicle's latitude.
where v N and v E represent vehicle's north and east velocity components, respectively; R = 6, 378, 137 m is the equatorial radius of Earth; and h is the height of vehicle above the Earth's surface.
Taking into account all the factors discussed above, the rate of change of vehicle's north, east and down (NED) velocities in a navigation frame can be described as: The position of the navigating vehicle can be expressed in terms of latitude (L), longitude (l), and height (ḣ) above the Earth's surface as:L The whole principle of strap-down inertial navigation in the local navigation frame is demonstrated by the schematic diagram in Figure 1. The GPS model considered in this research can be found in [17].

Filtering
Development of a reliable aided inertial navigation system depends on selecting an appropriate estimation method. The quality of parameter estimates depends chiefly on the type and accuracy of inertial sensors. Thus, the information from sensors should adequately be formulated with careful consideration of error characteristics in each sensor. Based upon the type of coupling, there are two basic GPS/INS integration modes; namely, loose and tight coupling modes [18]. The loosely coupled integration is the simplest and most robust technique of coupling. Navigation solutions are generated by both the INS and the GPS receiver, independently. A third navigation solution is then derived by using an estimator to merge the information provided by the two independent solutions. The Kalman filter including its derivatives is the most widely used filter to realize this fusion [19][20][21]. The fusion Kalman filter calculates error estimates between the two solutions, which are fed back to provide corrected attitude, velocity and position. The robustness of the loosely coupled scheme is due to the fact that if one system fails, the other sensor can still provide navigation, although with reduced accuracy. This research work utilizes the loose coupling strategy for INS-GPS fusion due to its simplicity, robustness, and lesser computational load.

Conventional Kalman Filter
The Kalman filter is a mathematical estimator used for predicting a system's behavior from noisy observations. It is a recursive algorithm which uses a series of steps, involving prediction and updating, to obtain an optimal minimum variance estimate of a system or state vector. It is commonly used for estimation in the field of navigation, particularly for autonomous applications [22,23]. It combines noisy sensor outputs (such as those from GPS, accelerometers, and gyroscopes measurements) and estimates the dynamic state of a system (which may include a vehicle's attitude, acceleration, velocity, position, etc.). For a discrete-time system, the Kalman filter algorithm is given by the following expressions: Prediction step:x Measurement update step: where Φ is the n × n state transition matrix;x k|k andx k+1|k represent the corrected(or a posteriori) value and the predicted (or a priori) value of the estimated state vector x, respectively; P k+1|k+1 and P k+1|k are the a posteriori and a priori state covariance matrices of order n × n, respectively, at time t k ; and the matrix K k+1 of order n × m is the Kalman gain. The term [z k+1 − H k+1xk+1|k ] is called the residual or measurement innovation, and H k is the measurements matrix of order m × n. The matrices Q (of order n × n) and R (of order m × m) represent the system and measurement noise covariances, respectively.
Assuming that system noise (w k ) and measurement noise (v k ) are uncorrelated and have characteristics of white noise, we can write the mean matrices as: And covariance matrices can be written as:

Adaptive Kalman Filtering
Due to its reliance on a priori information (of dynamic model and noise parameters), the Kalman filter may suffer from degradation or even divergence in situations where this a priori information is time-varying (or unknown) [24]. The matrix R usually contains information provided in a sensor's data sheet, while the matrix Q is attained analytically. Thus, any errors in Q and/or R may lead to sub-optimal estimation. Although usually taken as constants, Q and R for most practical situations must be adjusted to incorporate time-varying states and sensor inaccuracies which also vary under different dynamic conditions [25,26]. Similarly, in situations where system states are changing abruptly, robustness of the filter also becomes very critical, so that the filter can keep up with sharp changes in observations/measurements. This is achieved by adjusting Kalman gain K, or process covariance matrix P [27]. In order to achieve strong convergence, filter accuracy may have to be suppressed in such cases. This type of filtering is known as adaptive Kalman filtering

Sage-Husa Adaptive Kalman Filter
The Sage-Husa Adaptive Kalman filter [28] is based on a maximum a priori estimation technique. Since the old measurements can influence the estimation process, the filter performance can be optimized by utilizing a finite set of recent measurements. This is achieved by assigning weights to Q and R matrices, which causes the filter to adjust the gain K adaptively. So, in situations where noise statistics vary with time, increasing the weights of latest innovations will lead to enhanced estimation accuracy. For systems with time varying noise characteristics, we can write the mean matrices as Filtering functions for Sage-Husa-Kalman filter are: Measurement update step: The parameter d is known as the attenuation factor, given by where b is called the forgetting factor. Generally, b is selected using trial and error method but a value between 0.95 and 0.99 has been found to work well [29] to converge the filter to optimal Q and R.

Strong Tracking Robust Kalman Filter
When dealing with inaccurate linear or non-linear dynamic system models, or in situations where state changes are more sudden, the strong tracking filter may provide another effective method for adaptive state estimation [27]. Although lacking in extensive theoretical analysis, the filter has widely been used in practical applications. It makes use of a time varying fading factor which is computed from the most recent measurement innovation, and dynamically corrects the error covariance for state prediction, providing strong tracking of a maneuvering target with abrupt state changes.
The strong tracking Kalman filter algorithm updates the system covariance matrix P by suppressing the filter accuracy, but making it more robust for real time situations, as generalized below: where the parameter λ is called the fading factor and is calculated upon the arrival of new measurements, as proposed in [30]. where and C k+1 is defined as and v o(k+1) is the measurement innovation given by The value of parameter ρ is generally taken to be 0.95-0.99. The filter is active only when the system is experiencing abrupt state change (λ > 1); that is, when

Proposed Scheme
A hybrid filter based on Sage-Husa and strong tracking Kalman filters is presented which offers balance between accuracy and robustness in practical real-time situations where variation in noise statistics is also accompanied by sudden changes in system states. When state change is not abrupt (λ ≤ 1), the filter will prioritize estimation by Sage-Husa algorithm; otherwise, when λ > 1, estimation will be carried out through strong tracking filter. The flow diagram demonstrates the filter operation. The filter has been implemented using a loosely coupled scheme due to its simplicity and robustness. The proposed algorithm is depicted in Figure 2.

Mathematical Model
In state-space representation for a linear dynamic system, a system of n-order differential equations is transformed into a set of n-coupled, first-order differential equations. Similarly, for the case of INS-GPS integration, the differential equations representing INS errors can be expressed in a state-space form, which will then be used as system model for Kalman filter. The state-space model consists of two equations; i.e., a dynamic equation and a measurement equation, which can be written respectively as:ẋ = Fx + w (46) where x is state vector of order n × 1, F is called system dynamic matrix with order n × n; w, an n × 1 vector, represents the system noise; z and v are called the observation vector and the measurement noise vector, respectively, both of order m × 1. H is the measurement matrix of order m × n, with m measurements and n states. The above equations can be written in discrete form as: where x k is the state vector of order n × 1 at time t k ; Φ k+1 is the n × n state transition matrix at time t k+1 ; w k is the system noise of order n × 1; z k is the observation vector of order m × 1; H k is the measurements matrix of order m × n; and v k is the m × 1 measurement noise vector at time t k . For a discrete-time system, the state transition matrix can be written as: Taylor series expansion of this equation over a small time interval t yields the following expression (ignoring the higher terms): where I is the identity matrix. The system matrix F is a 9 × 9 matrix formed by incorporating the attitude, position and velocity error models of the inertial sensors, which are explained by Titterton in [31]. The system matrix F is given by where F 1 through F 9 are all 3 × 3 matrices given as where f N , f E , and f D are the north, east, and down components of specific force, whereas, v N , v E , and v D are the north, east, and down components of vehicle's velocity, respectively. Ω = 15.041067 • /hr represents the Earth's angular speed. r N = R N + h and r E = R E + h, where R N and R E are known as Earth's meridian radius of curvature and transverse radius of curvature, respectively, when an ellipsoidal model is considered for Earth. S L , C L , and T L represent sin, cos, and tan of the vehicle's latitude (L), respectively.
where e is the eccentricity of the Earth, with a constant value of e = 0.08181919082423, and R e is known as the equatorial radius of the earth, with value R e = 6378137 m. The state vector x is given by where δα, δβ, and δγ are the errors in roll, pitch, and yaw; δv N , δv E , and δv D are the errors in north, east, and down velocities; and δL, δl, and δh are the errors in latitude, longitude, and altitude, respectively. The system noise matrix Q consists of variance errors in its principal diagonal and is formed directly from the IMU sensor data sheet, whereas the uncertainties in the GPS measurements are used to derive the measurement noise covariance matrix R.
The difference between GPS and IMU velocity and position measurements forms the measurement vector z. (67) The latitude and longitude measurements are in radians with very small values, which can make the measurement matrix H singular. This can cause numerical instabilities in Kalman gain K calculation. So the measurement vector can be modified as: The measurement matrix (H) is written as: All navigation calculations are being carried out in NED frame of reference.

Simulation, Data Collection and Results
The field test drive for data collection each lasting for fifteen minutes was carried out or Murree road (34.156024 • N, 73.222560 • E), Abbottabad city, Pakistan in January 2019. A total of twenty GPS and IMU measurement data sets were gathered over a period of three days. Two of the data sets were discarded due to erroneous data as a result of a hardware malfunction. Due to mountainous terrain of the city, satellite coverage was moderate and the average number of satellites available was equal to six. In order to avoid any timing errors between the GPS and IMU data, software written in MATLAB was used to time stamp the MEMS IMU data.
The YH-7000 commercial grade MEMS IMU was used for inertial measurements, which include angular rates and specific force data provided by three dimensional gyros and accelerometers. IMU data was obtained at 200 Hz. Simulations have been carried out on real vehicle test drive data in MATLAB R2014A platform (3.1 GHz, Intel, Core i5 CPU). In order to better visualize the filter's performance for highly dynamic systems, Gaussian noise was introduced to the measurements during the interval between 300 s and 500 s. The actual GPS data was used to compare the estimation results. A total of 250 Monte Carlo simulations were carried out to omit any bias introduced by adding the noise.
Novatel's OEMV series GPS receiver was used to obtain raw data which was converted to position and velocity measurements using the on-board software at 1 Hz. The GPS measurements included vehicle's position in terms of latitude, longitude, and altitude, and vehicle's velocity in terms of north, east and down velocities. Typical error of commercial GPS data is about 30 m in latitude and longitude, while error in altitude is about 1.5 times greater than latitude/longitude error. The YH-7000 commercial grade MEMS IMU Sensor specifications are given in Table 1: The results for the conventional linear Kalman filter and the hybrid adaptive Kalman filter for INS-GPS integrated navigation are now compared, as shown in the following simulation results with actual sensor data. Vehicle position in terms of latitude, longitude, and altitude is shown first in  The proposed hybrid filter produces better estimates for vehicle's position and velocity. Superior performance of the hybrid filter is more evident during the interval between 300 s and 500 s, when random noise was added to the actual sensor measurements. It converges more rapidly compared to the conventional Kalman filter. Similarly, when there is abrupt change in vehicle's latitude and longitude measurements (near time = 650 s), the hybrid adaptive filter clearly excels at convergence (Figures 3 and 4). This behavior is more pronounced for vehicle's velocity estimation, which is understandable, given the fact that measurements have been taken for a ground vehicle, which generally experiences more abrupt changes in its velocity/acceleration rather than its position. The moving mean square error (MMSE) plots provided in Figures 9-15 offer a better insight into the performance of proposed hybrid adaptive filter compared to conventional Kalman filter. These curves clearly illustrate superior estimation accuracy of the proposed filter scheme. A quantitative measure of the filter's accuracy has also been provided in terms of accumulative error plots (Figures 15 and 16), which show overall improvements made by the proposed hybrid scheme in vehicle's position and velocity estimates.
These results prove the efficiency of the proposed scheme for parameter estimation in navigation, especially for the case when measurements contain large uncertainties/errors, or when states are changing more abruptly, such as the vehicle's velocity.   Switching decision between the Sage-Husa filter and the robust tracking filter is made based on the check provided in Equation (45). When an anomaly in the measurements is detected, like a large error, or abrupt state change, the condition in Equation (45) is satisfied, and the filter switches to robust mode where more weight is assigned to latest measurements, thus, providing strong tracking of the target.

Conclusions
A novel algorithm based on a hybrid of the Sage-Husa filter and the strong tracking filter using a loosely coupled scheme is proposed for GPS-INS integration. The proposed algorithm has generated superior navigation for a vehicle's position and velocity estimates. The algorithm offers a balance between accuracy and strong convergence for optimum estimation for highly dynamic systems. It inherently performs much better when implemented for more dynamic systems, as proven by the results. The same hybrid filter will, thus, produce even better navigation estimates for target tracking applications or in case of aerial vehicles, where a conventional Kalman filter will certainly lag behind due to its reliance on a priori information, and consequently, its inability to keep track of sudden changes in system states. The hybrid filter based on Sage-Husa-Kalman filter and strong tracking Kalman filter can, thus, be acknowledged as an optimal filter when compared to the conventional Kalman filter.