INS/GNSS Tightly-Coupled Integration Using Quaternion-Based AUPF for USV

This paper addresses the problem of integration of Inertial Navigation System (INS) and Global Navigation Satellite System (GNSS) for the purpose of developing a low-cost, robust and highly accurate navigation system for unmanned surface vehicles (USVs). A tightly-coupled integration approach is one of the most promising architectures to fuse the GNSS data with INS measurements. However, the resulting system and measurement models turn out to be nonlinear, and the sensor stochastic measurement errors are non-Gaussian and distributed in a practical system. Particle filter (PF), one of the most theoretical attractive non-linear/non-Gaussian estimation methods, is becoming more and more attractive in navigation applications. However, the large computation burden limits its practical usage. For the purpose of reducing the computational burden without degrading the system estimation accuracy, a quaternion-based adaptive unscented particle filter (AUPF), which combines the adaptive unscented Kalman filter (AUKF) with PF, has been proposed in this paper. The unscented Kalman filter (UKF) is used in the algorithm to improve the proposal distribution and generate a posterior estimates, which specify the PF importance density function for generating particles more intelligently. In addition, the computational complexity of the filter is reduced with the avoidance of the re-sampling step. Furthermore, a residual-based covariance matching technique is used to adapt the measurement error covariance. A trajectory simulator based on a dynamic model of USV is used to test the proposed algorithm. Results show that quaternion-based AUPF can significantly improve the overall navigation accuracy and reliability.


Introduction
The development of unmanned surface vehicles (USVs) for scientific, military and commercial purpose in applications such as oil and gas exploration, oceanographic data collection, hydrographic, oceanographic and environmental survey, mine counter measure, surveillance and reconnaissance, anti-submarine warfare and fast inshore attack craft for combat training require the corresponding development of navigation systems [1]. The need of the robustness, accuracy and reliability is a guarantee of long duration, unmanned operations for USVs. The sensor units on board, which serve as environment perceptions, are of critical importance for automation operation. A reliable and accurate navigation system is very important for a USV.
The integration of Global Navigation Satellite System (GNSS) and Inertial Navigation System (INS) is widely used for positioning and attitude determination for vehicles. The development of micro-electromechanical system (MEMS) technology has brought low-cost INS/GNSS integration approaches into practice [2,3]. A growing number of research groups are developing integrated navigation systems utilizing INS and GNSS due to the complementary nature of INS and GNSS principles. The INS/GNSS integration can be classified into loosely-coupled, tightly-coupled and deeply-coupled [4]. For the loosely-coupled manner, independent redundant solutions are available from the GNSS receiver and INS process. However, the disadvantage is that typically four satellites have to be in view to obtain position and velocity solutions from the GNSS receiver. In addition, cascaded filtering problems may occur when one local Kalman filter (KF) is used in GPS data processing and the other is used in integration. These problems can be easily solved in tightly-coupled integration. In tightly-coupled integration, a centralized KF is employed. The pseudorange and delta pseudorange (or Doppler) measurements are used as observations to update the navigation filter. Furthermore, the system does not need a full GNSS solution to assist the INS. This means that the error correlation of inertial measurement unit (IMU) measurements can maintain the update of the estimated solution by the INS even if the number of the tracked satellites is less than four. However, the resulting system observation model turns out to be nonlinear, which should be carefully treated in the design of the integration KF.
The EKF (extended Kalman filter) is known as state-of-the-art for fusion INS and GNSS data in tightly-coupled integration. The linearization should be implemented in both the nonlinear system model and the observation model first in order to apply EKF. This approximation will result in large errors when EKF calculates a posterior mean and covariance, which may lead to suboptimal performance or even divergence of the filter. UKF (unscented Kalman filter) is a recursive MMSE (Minimum Mean Square Error) estimator based on optimal Gaussian approximation Kalman filter framework. Unlike the EKF, the UKF uses the true nonlinear model and approximates the distribution of the state random variables [5]. The state distribution in the UKF is specified using a minimal set of deterministically chosen sample points to capture the posterior mean and covariance accurately to 2nd order for any nonlinearity. It is based on the assumption that it is much easier to approximate a Gaussian distribution than to simulate an arbitrary nonlinear function [6]. Recently, with the development of the computer technology, particle filter (PF) turns out to be more attractive for nonlinear and non-Gaussian applications, and has been successfully used in [7] to recursively update the posterior distribution by sequential importance sampling and resampling. However, the large computational burden impeded the practical use of PF. In addition, the sample impoverishments accompanying the degeneracy of the system performance are primary disadvantages of the basic PF [8]. To overcome these problems, in [3], the strategy of combining the UKF with PF was proposed.
To achieve better performance from the KF framework, the stochastic information provided to the filter must be as accurate as possible. It is therefore necessary to adapt the stochastic model to accommodate for changes in vehicle dynamics and environment conditions. Insufficient or incorrect knowledge about statistics may lead to degradation of the system performance or even divergence of the filter. In [9], a filter innovation sequence based on adaptive Kalman filtering (AKF) was introduced, which showed a major improvement in adjusting process noise error covariance and sensor measurement error covariance adaptively. In this paper, an adaptive unscented particles filter (AUPF) algorithm based on quaternion is proposed using a residual-based covariance matching technique.
In the remainder of this paper, the content is organized as follows. In Section 2, a trajectory simulator with corresponding sensors measurements allocated in a USV is developed. In Section 3, quaternion-based propagation and an observation model are introduced. Furthermore, a quaternion-based GPS/INS integration algorithm using AUPF is proposed. Simulation results are analyzed and compared to illustrate the performance of the proposed AUPF algorithm.

Sensor Measurement Simulation
In this part, a USV trajectory simulator is designed. We briefly summarize the different component parts of the USV platform. These include the mathematical dynamic state-space models used for the simulator, the sensors subsystem and the controller design. The core components of such a simulator are schematically depicted in Figure 1, which include a control system, a guidance system and a navigation system (GNC). The GNS system takes the noisy GPS and MEMS-IMU measurements as inputs, and then fuses them from the dynamics model of USV to estimate the optimal vehicle navigation state solutions. There state estimations together with desired trajectories from guidance systems are then adopted by the controller, which generates an optimal control law to drive the thruster of the USV. optimal vehicle navigation state solutions. There state estimations together with desired trajectories from guidance systems are then adopted by the controller, which generates an optimal control law to drive the thruster of the USV.

Trajectory Simulator of the USV
As in [10], the marine craft equations of motion in six degrees of freedom can be written in vectorial setting form: where ∈ denotes the position and the orientation vector. ν ∈ denotes the linear and angular velocity vectors that are decomposed in the body-fixed reference frame. τ ∈ describes the forces and moments acting on the craft in the body-fixed fame. The generalized position, velocity and force vectors have the form that represented as Equation (3): D v represent the inertial, Coriolis-Centripetal and damping matrices, respectively.   g η represents the restoring forces and moments.
A sliding mode trajectory tracking controller is designed to track the reference trajectory for USV as in [11]. All relevant position, orientation, linear and angular velocities, acceleration and forces describing the USV's trajectory are calculated.

GPS Date Simulation
There are three kinds of measurements from GPS receiver, i.e., pesudorange, Doppler and carrier phase. In this subsection, we will introduce the measurement mode of pesudorange and Doppler. The first step is to develop the GPS constellation model, which can be used to generate the position of the satellites in the simulation. GPS archive data are available from the website of the International GNSS service (IGS). The next step is to model the signal transmission from the satellites. The main GPS measurement errors include ionospheric errors, tropospheric errors, etc. Further details can be found in [12]. Considering all of these aspects, the biased pseudorange measurement from one satellite vehicle at one instance is formulated as:

Trajectory Simulator of the USV
As in [10], the marine craft equations of motion in six degrees of freedom can be written in vectorial setting form: where η P R 6 denotes the position and the orientation vector. ν P R 6 denotes the linear and angular velocity vectors that are decomposed in the body-fixed reference frame. τ P R 6 describes the forces and moments acting on the craft in the body-fixed fame. The generalized position, velocity and force vectors have the form that represented as Equation (3): The same notation as in [10] is used, p n b{n P R 3 is the position expressed in the north east down (NED) frame, {N}. Θ nb P S 3 represents the Euler angles. υ b b{n and ω b b{n represent the linear and angular velocity of body frame expressed in{B} frame with respect to {N} frame. f b P R 3 and m b P R 3 are the forces and moments acting on the vehicle, respectively. R 3 and S 3 denote the Euclidean space of dimension three and the sphere, respectively. J Θ pηq is the transformation matrix. M, C pvq, and D pvq represent the inertial, Coriolis-Centripetal and damping matrices, respectively. g pηq represents the restoring forces and moments.
A sliding mode trajectory tracking controller is designed to track the reference trajectory for USV as in [11]. All relevant position, orientation, linear and angular velocities, acceleration and forces describing the USV's trajectory are calculated.

GPS Date Simulation
There are three kinds of measurements from GPS receiver, i.e., pesudorange, Doppler and carrier phase. In this subsection, we will introduce the measurement mode of pesudorange and Doppler. The first step is to develop the GPS constellation model, which can be used to generate the position of the satellites in the simulation. GPS archive data are available from the website of the International GNSS service (IGS). The next step is to model the signal transmission from the satellites. The main GPS measurement errors include ionospheric errors, tropospheric errors, etc. Further details can be found in [12]. Considering all of these aspects, the biased pseudorange measurement from one satellite vehicle at one instance is formulated as: where r ρ denotes a measured range of user to satellite, ρ is a true range of user to satellite, t u and t s denote receiver clock bias and satellite clock bias, respectively, T iono and T tropo denote ionospheric delay and tropospheric delay, respectively, ε ρ denotes other un-modeled errors, i.e., multipath delay.
The relative motion between the satellite and the receiver results in change of the observed frequency of satellite signal. The Doppler can be used to estimate the user velocities from the satellite velocities. As in [3], the Doppler shift can be written as a projection of the relative velocity vector on the satellite line-of-sight vector: where I s u is the use-to-satellite line-of-sight unit vector, v s and v u represent the satellite and receiver velocity respectively. I s u can be expressed as: where p s " px s , y s , z s q and p u " px u , y u , z u q denote the position vector of satellite and receiver expressed in Earth-Central Earth Fixed (ECEF) coordinate.

IMU Date Simulation
For the purpose of simulating IMU data, raw measurements of accelerometer and gyroscope are needed. The trajectory simulator data are used as the basis in simulating the sensor data. Moreover, the velocity and angular rate data for the USVs are used.
The IMU measurements are provided by extracting the acceleration and angular velocity from the simulator model of USVs, which can be modeled as: where g n represents gravity expressed in a navigation frame. ε b acc and ε b gyro denote the zero mean Gaussian distribute noise. f bias b and ω bias b are bias errors of the specific forces and angular rate measurements, respectively. For the low cost MEMS-IMU, the sensor errors have the non-Gaussian characteristics. The bias errors in the gyroscope and accelerometer measurements in body frame will be transformed to be the position and velocity drifts in the navigation frame.

Quaternion-Based Propagration and Observation Models
Define the vector form of quaternions as q " " q 0 q T ı T with one real part q 0 and three imaginary parts given by the vector q " Quaternions are represented as a complex number with four bases and are used to compute the rotation from navigation frame to body frame. Based on Euler's theorem, every change in the relative orientation of two rigid bodies or reference frame can be produced by means of a simple rotation from one frame to another along fixed axes. Given the invariant axis (rotation axis) and rotation angle (with magnitude ||ϕ||), the quaternion vector can be represented as: where u denotes the unit vector along the invariant axis. Apparently, q has the normality property, that is to say ||q|| " 1. It can be concluded that the quaternion vector has only three degrees of freedom, although q has four elements. Given a rotation vector ϕ " rϕ 1 , ϕ 2 , ϕ 3 s T , the quaternion vector can be computed as: Conventionally, with a sufficiently small time interval, the quaternion vector is updated using vectors added together in discrete time domain. However, it should be noticed that the unit sphere defined by quaternion is not a Euclidean vector space. That is to say, the common definition of addition and scaling cannot be applied directly. In the AUPF algorithm, we will apply the quaternion product rule to update the quaternion vector. Thus, the system propagation model in discrete time domain can be expressed as: (11) where b denotes the quaternion product. p n,k and v n,k denote the position and velocity in the navigation frame at epoch k. g n represents gravity expressed in navigation fame, which is assumed to be constant for local navigation. c∆ . t k denotes the receiver clock drift error, which is modeled as a random walk process. c∆t k represents the range equivalent of the receiver clock bias, which is the integration of clock drift error. T is the system propagation time interval. w f , w ω , w cbk and w cdk are Gaussian nose terms. R n b pqq denotes the rotational transformation matrix from the body frame to the navigation frame, which can be formulated using the parameters of quaternion as: ffi fl (12) ∆q k represents the quaternion rotation of body frame during the time interval: where θ k denotes the integral of the body frame angular rate measurements. Equation (11) forms the propagation model of the system. In the tightly-coupled integration, the system observation model is formulated as Equation (14). Where j denotes the number of satellites in view, I j represents the estimated line-of-sight unit vector pointing from the initial estimate of the user position to the j-th satellite: whereρ j,k is the predicted pseudorange measurement from the j-th satellite, and: wherex n u,k ,ŷ n u,k ,ẑ n u,k are the USV position estimates expressed in NED frame. x n j,k , y n j,k , z n j,k are the j-th satellite position coordinates expressed in NED frame.
The predicted pseudorange-rate measurements are calculated as: .
where I n,n j,k "´x n j,k´x n u,k¯{d j,k , I n,e j,k "´y n j,k´ŷ n u,k¯{d j,k , I n,d j,k "´z n j,k´ẑ n u,k¯{d j,k d j,k " c´x n j,k´x n u,k¯2`´y n j,k´ŷ n u,k¯2`´z n j,k´ẑ n u,k¯2 whereˆ. ρ j,k is the predicted delta range measurement from the j-th satellite.ˆ. Due to introducing the quaternion, there exists dimension mismatch among the state vector and the state error covariance matrix. The reason is that the degree of freedom of a quaternion vector is three rather than four. When the state vector includes quaternion vector elements, the dimension of the state vector is 18ˆ1; however, the dimension of the state error covariance matrix is 17ˆ17. In order to solve the problem of dimension mismatch, we introduce the rotation vector in the rotation space, which is transformed from the corresponding quaternion vector error.

Quaternion-Based Integration Using UPF
The idea of UKF comes from the fact that it is much easier to approximate a Gaussian distribution, rather than to simulate a nonlinear function. The particles for UPF are generated using the UKF a posterior estimates. The state vector and its associated errors are defined as: whereφk is the rotation vector corresponding to qk b pq k q´1. The state error covariance Pk can be formulated as: The sigma-points are generated according to Equation (19) where 2n pn " 17q equally weighted sigma-points are generated. The calculated sigma-points has the form of:xk´1 Define the delta rotation angle rate in the quaternion form as: The a priori mean and the state error covariance matrix are computed through Equations (22) and (23): The same as in Equation (19) xk "xḱ`K k pr y k´ŷk q Pk " Pḱ´K k P yy For the UPF, the a posterior estimates ofxk and Pk from Equation (27) are used to form the importance density distribution for generating particles:

/ -
where m denotes the observer dimension which varies during the time and is based on the number of satellites being tracked by the receiver. The a posteriori mean and covariance estimates are computed as: ϕk ,i is the rotation vector corresponding toqk ,i b`qk˘´1. Due to the fact that the unit quaternion is not mathematically closed for addition and scalar multiplication, a normalization procedure is necessary to ensure the constraint q T q " 1 is satisfied in the presence of measurement noise and numerical round-off errors. The normalization can be conducted by replacingqk withqk {||qk ||.

Residual-Based AUPF
For the a posterior state estimations, a reliable resolution is proposed which greatly depends on prior statistics of the system process and measurement noise. However, these statistics are hardly known exactly in practice because they are based on the types of applications and process dynamics [13]. In addition, the estimation environment regarding INS/GNSS kinematic application is not always fixed but subject to change [14]. In order to maintain the estimation accuracy, a residual-based covariance matching technique is used to adaptively adjust the algorithm parameters and to track the changes in the noise source. For the INS/GNSS integration, the IMU measurements are used in system model, and the product manual specifies the sensor turn-on bias, temperature related variations in biases and noises. In this paper, we assume that the process noise error covariance Q is approximately known, and we focus on the adaptive estimation of sensor measurement error covarianceR.
Different from [2], which proposed an adaptive EKF algorithm, for the UKF, modified measurement noise error covariance matrices are adaptively updated as: where P yyk is calculated as: whereŷk " 1 2n 2n ř i"1 h k´xk ,i¯. Thexk ,i denotes the sigma-points drawn from the UKF a posterior estimates. The sigma-points are generated as: The calculated sigma-points have the form of Equation (35). The delta rotation angle rate in the quaternion form can be formulated as: The flowchart of the proposed AUPF is illustrated in Figure 2.
The flowchart of the proposed AUPF is illustrated in Figure 2.

Simulation Results
A simulator with the module of USV trajectory simulation and the raw measurements simulation for IMU and GPS was developed. The trajectory used in this simulation scenario is illustrated in Figure 3. We assume that the USV is operating at sea level. The sliding mode controller in [15] drives the USV fallowing the designed trajectory. To represent a real MEMS-IMU as closely as possible, different error types including the Gaussian model and random bias were induced in the simulated IMU. Raw GPS pseudorange and Doppler measurements were simulated at a 1 Hz rate, while the IMU raw data were simulated at a 100 Hz rate. The sensor error characteristics of MEMS-IMU are given in Table 1. The program is implemented in MATLAB R2013a using Intel Xeon E31270 CPU and 8 GB RAM.

Simulation Results
A simulator with the module of USV trajectory simulation and the raw measurements simulation for IMU and GPS was developed. The trajectory used in this simulation scenario is illustrated in Figure 3. We assume that the USV is operating at sea level. The sliding mode controller in [15] drives the USV fallowing the designed trajectory. To represent a real MEMS-IMU as closely as possible, different error types including the Gaussian model and random bias were induced in the simulated IMU. Raw GPS pseudorange and Doppler measurements were simulated at a 1 Hz rate, while the IMU raw data were simulated at a 100 Hz rate. The sensor error characteristics of MEMS-IMU are given in Table 1. The program is implemented in MATLAB R2013a using Intel Xeon E31270 CPU and 8 GB RAM.    Figure 3 shows the estimation results of trajectory of USV in horizontal plane. Although the low cost MEMS-IMU is used, the agreement with the actual trajectory is generally good. The estimated position and velocity errors can be found in Figure 4. After the initial transient, the estimation error in the north direction remains less than 3 m approximately, 0.5 m in the east direction and −4.2 m in the down direction; the errors of velocity remain within a region of approximately ±0.05 m/s. In Figure 5, the attitudes estimated by the AUPF algorithm (using quaternion) are plotted. The quaternion are converted to their corresponding Euler angles. Except the initial transient, the roll and pitch errors remain mostly in the region ±0.1 deg, but the yaw error tends to be a little larger. Accelerometer bias estimation results are illustrated in Figure 6. Gyroscope bias estimation results are plotted in Figure 7. Due to the fact that the yaw angle is the least observable state in practical vehicle motion, we can see form Figure 7 that the bias error on the z-axis takes more time to be stable.   Figure 5, the attitudes estimated by the AUPF algorithm (using quaternion) are plotted. The quaternion are converted to their corresponding Euler angles. Except the initial transient, the roll and pitch errors remain mostly in the region˘0.1 deg, but the yaw error tends to be a little larger. Accelerometer bias estimation results are illustrated in Figure 6. Gyroscope bias estimation results are plotted in Figure 7. Due to the fact that the yaw angle is the least observable state in practical vehicle motion, we can see form Figure 7 that the bias error on the z-axis takes more time to be stable.              In order to illustrate the performance of the proposed AUPF algorithm, several groups of tests adopting the AEKF, AUKF, PF and AUPF are conducted, respectively, in this part. For comparing the navigation accuracy among the algorithms mentioned above, we choose the norm of the position (||∆p||) and velocity (||∆v||) estimation errors for analysis, which are calculated as As shown in Figure 8, AUPF presents best estimation accuracy compared with others.
of the norm of the position and velocity estimation errors. The time represents the algorithm processing time, which depends on the computation power of the computer and only for reference. We conduct one simulation run for AEKF and AUKF because the same prior statistical parameters and the same group of sensor data result in unchanged performance. For the PF based algorithm, the particles are generated randomly, so the estimation results have slight differences after each simulation. Table 2 shows the average values of 100 runs of the algorithms. We can see that AUPF (100 particles) can significantly improve the accuracy compared with PF (300 particles), and with computational efficiency at the same time.
Although the validation of the proposed AUPF algorithm is conducted from simulation in this paper, it has a theoretical foundation with a stringent mathematical representation. However, a real data experiment is an intuitively clear way to validate the performance of the presented AUPF algorithm, and it is of great practical interest to validate and compare it with other filters in experiments. Furthermore, with the same group of GPS and IMU data, the estimation accuracy and computation time of the algorithms are illustrated in Table 2. M(||∆p||) and M(||∆v||) denote the mean of the norm of the position and velocity estimation errors. V(||∆p||) and V(||∆v||) denote the variance of the norm of the position and velocity estimation errors. The time represents the algorithm processing time, which depends on the computation power of the computer and only for reference. We conduct one simulation run for AEKF and AUKF because the same prior statistical parameters and the same group of sensor data result in unchanged performance. For the PF based algorithm, the particles are generated randomly, so the estimation results have slight differences after each simulation. Table 2 shows the average values of 100 runs of the algorithms. We can see that AUPF (100 particles) can significantly improve the accuracy compared with PF (300 particles), and with computational efficiency at the same time.
Although the validation of the proposed AUPF algorithm is conducted from simulation in this paper, it has a theoretical foundation with a stringent mathematical representation. However, a real data experiment is an intuitively clear way to validate the performance of the presented AUPF algorithm, and it is of great practical interest to validate and compare it with other filters in experiments. Adaptive extended Kalman filter (AEKF); Adaptive unscented Kalman filter (AUKF); Particle filter (PF); Adaptive unscented particle filter (AUPF).

Conclusions
In this paper, we have studied the problem of navigation system design for low-cost USVs, in which we mainly focus on the filtering algorithm. A quaternion-based tightly-coupled integration approach is chosen to fuse the low-cost MEMS-INS measurements and GPS data. Due to the nonlinear and non-Gaussian properties of the INS/GNSS tightly-coupled integration problem, an APUF algorithm is proposed, which combines UKF with PF. A residual-base covariance matching technique is used to adaptively adjust the parameters and to track the changes in observation data in an adaptive manner. The resulting AUPF algorithm can reduce the computational burden without degrading the system estimation accuracy. A USV trajectory simulator is designed to generate the raw senor measurements of IMU and GNSS. The simulation results verify the robustness and accuracy of the AUPF algorithm in comparison to the AEKF and AUKF. In addition, the comparison with PF underlines the computational advantage of AUPF.