Non-Linear Filtering for Precise Point Positioning GPS / INS integration

This research investigates the performance of non-linear estimation filtering for GPS-PPP/MEMS-based inertial system. Although integrated GPS/INS system involves nonlinear motion state and measurement models, the most common estimation filter employed is extended Kalman filter. In this paper, both unscented Kalman filter and particle filter are developed and compared with extended Kalman filter. Tightly coupled mechanization is adopted, which is developed in the raw measurements domain. Un-differenced ionosphere-free linear combination of pseudorange and carrier-phase measurements is employed. The performance of the proposed non-linear filters is analyzed using real test scenario. The test results indicate that comparable accuracy-level are obtained from the proposed filters compared with extended Kalman filter in positioning, velocity and attitude when the measurement updates from GPS measurements are available.


Introduction
Global satellite navigation systems such as the global positioning system (GPS) are often used for navigation applications.With the lack or complexity of attitude determination and the accuracy degradation during partial outages, GPS can be integrated successfully with independent environmental navigation system such as inertial navigation system (INS).The motivation of the proposed integrated system is the aiding of GPS observations in compensating INS errors which accumulated in integrating process, on the other hand, during complete and partial GPS outages, the navigation parameters can be continue obtained from the INS.Differential GPS technique are commonly used in GPS/INS integration in higher accuracy-level applications.This is mainly due to the high accuracy of differential GPS in comparison with standalone mode.Unfortunately, this involves the deployment of a base station, which limits the range of navigation area and increases the cost and complexity of the system.With the development of precise point positioning PPP (Zumberge et al, 1997 andKouba andHeroux, 2001), which is capable of providing decimeter to centimeter positioning accuracy without the need for a base receiver; it is possible to develop a high accuracy GPS/INS system based on one GPS receiver only.In addition, the great advances in micro-electro-mechanical sensors (MEMS) provide the development of a generation of low cost inertial sensors.MEMS sensors are characterized by small size, light weight and low cost with respect to high-end inertial sensors.Generally, MEMS sensors have poorer performance compared with high-end IMU due to the severe biases and errors affecting inertial sensors.
For the integrated GPS/INS system, the ultimate target is to estimate the navigation state of the moving body at the current time based on sensor readings, given a set of updated measurements (GPS observations) collected at time steps.The Bayesian filtering addresses the estimation problem of the dynamic model by propagating the conditioned probability density function (PDF) of the navigation state variables f x|Z from one epoch to the next epoch, knowing the PDF of the observationsf z|x .This is applied into two steps; prediction step and update step.In the prediction step, the dynamic motion model is considered based on IMU readings.Assuming that the PDF f x |Z is available at time step k-1, the motion model is used to predict the current navigation state variables of the moving body by computing the priori PDF f x |Z through integration process.In the update step, the measurement model is considered to incorporate new information (GPS observations) to update the priori density to obtain the posterior PDF f x |Z according to Bayes theorem (Jekeli, 2001).The prediction and update steps repeated recursively to obtain the estimated navigation states at time steps.The Bayesian filtering yields the optimal estimation solution without approximation for the integrated system states.However, this recursive solution yields multi-dimensional integrals.For nonlinear systems such as our integrated system, it is impossible to evaluate these recursive solutions analytically.
One of the solutions is to transform the nonlinear system to linear through linearization process by applying Taylor series expansion and neglecting the second and higher order terms assuming Gaussian distribution density which is applied in Kalman filter (KF) estimation.Generally, KF gives optimal estimation solution for linear models.In other words, for the integrated GPS/INS system, KF gives optimal solution for the approximate system (linearized system) rather than the original system (nonlinear system).Using first order Taylor linearization may cause divergence of motion models especially during GPS outages due to the impact of neglecting higher order terms especially when low cost MEMS-based IMU is used.The Unscented Kalman filter (UKF) which is first introduced as a non-linear filter by Julier et al, (1995), is considered as a linear regression Kalman filter, Unlike, Extended Kalman filter, UKF linearize the nonlinear systems through a linear regression between n points, drawn from the prior distribution of the random variable rather than applying Taylor series expansion linearization, called sigma points .A set of sigma points with appropriate weights is deterministically chosen propagated through the true non-linear system, captures the posterior mean and covariance (Bergman, 2001).
On the other hand, nonlinear filters such as Particle filter (PF) avoid the linearization of the system model, but instead, obtain approximate estimation solution for the nonlinear model rather than approximate the system model itself.PF has been used for GPS/INS integrations using different approaches (Yi and Brzezinska 2006, Giremus et al, 2005, Caron et al, 2007and Georgy et al, 2010).PF handles the multi-dimensional integrals numerically rather than the exact analytical solution.PF can accommodate arbitrary sensor characteristics, motion dynamics, and noise distributions, PF has no approximations for the system and measurement models the only approximation is the numerical solution by estimating the PDF P x|Z discretely by a set of N samples based on the Monte Carlo method (Ristic et al., 2002).However, as mentioned in (Doucet, 1998), it impossible to sample fromP x|Z , an importance Density is used instead.This method is called the importance sampling (SI) which based on selecting probability density which taking into account the observations to time k.Generally, in GPS/INS integration, the prior probabilistic motion density is used as an importance density.This makes the PF suffer specially when using low cost MEMS-based inertial sensors because this importance density does not produce enough samples in regions where the true PDF is large due to the high drift of these low cost sensors.This research aims to develop a new integrated navigation system based on integrating GPS-based PPP with MEMS accelerometers and fiber optic gyros for precise navigation applications.The proposed integrated system requires rigorous modeling of all errors and biases affecting both inertial and GPS observations.We use un-differenced observations for GPS PPP.Inertial sensor biases are accounted for through calibration, while Gaussian-Markov stochastic process model is used to account for the sensor's random errors.Both UKF and PF are employed to merge GPS and inertial sensors data compared with the traditional EKF.

GPS-PPP/MEMS-based IMU tightlycoupled algorithm
Tightly coupled (TC) architecture is implemented adopting a central filter to process GPS raw pseudorange, phase and Doppler measurements and the INS-derived observations to produce estimates of the state vector including position, velocity and attitude.
The mathematical modeling of INS is commonly described in the framework of linear dynamic systems.The dynamic behavior of such systems can be described using a state-space representation.For this purpose, a system of non-linear firstorder differential equations (mechanization equations) is employed.In the local-level frame (n-frame), this state vector for INS mechanization is represented by (Jekeli, 2001): Where; r is position vector; latitude, longitude and altitude.V is the velocity in the ENU frame, consisting of the three East, North and Up components, V is the kinematic acceleration in the ENU frame, Ω .V is the effect of the motion of the ENU frame with respect to the ECEF frame, 2Ω .V is the Coriolis acceleration.g is the gravity vector, including the Gravitation term and the centripetal term related to the Earth rotation and f is the specific force vector in the body frame and is measured by the accelerometers.The matrix Ω is the skew-symmetric matrix of rotation rate vector of Earth expressed in the ENU frame as; The matrix Ω is the skew-symmetric matrix of the rotation rate vector of the ENU frame with respect to ECEF frame, expressed in the ENU frame as: The matrix Ω is the skew-symmetric matrix of the rotation rate vector of the body frame with respect to the ECI frame ω , expressed in the body reference, which is measured by the gyros.The matrix Ω is the skew-symmetric matrix of the rotation rate of the navigation frame with respect to inertial frame ω expressed in the body frame which is computed combining ω and ω transforming the result in the Body frame as follows; The bias drift is modeled as a first-order Gauss-Markov process (Petovello, 2003) represented by the equations: where The subscript "i" indicates the axis, and are the correlation times for the accelerometers and gyros respectively, and and are the Gauss-Markov process driving noises, whose spectral densities are and .
The clock errors unique to the GNSS measurements such as clock offset and drift are modeled and written as (Brown and Hwang 1997): and w are the clock offset and drift driving noise with spectral density q and q respectively.The clock error spectral densities are computed as Brown and Hwang (1997).The complete state vector can be written as; The measurement model of the GPS/INS filter in the TC architecture has the typical form: (3) Where y(k) is the corrected un-difference ionosphere-free GPS measurements, h(x) is the nonlinear observation GPS satellites equation function of INS positioning (predicted measurements) and ε is the Gaussian white noise with zero mean and covariance matrix R.
The un-differenced ionosphere-free combination of code pseudorange and carrier phase mathematical model can be written as; where P and P are code measurements on L and L , respectively; Φ and Φ are the carrier phase measurements on L and L , respectively; dt and dt are the clock errors for receiver and satellite, respectively; T is the tropospheric delay; d and d are frequency-dependent code hardware delay for receiver and satellite, respectively; δ and δ are frequencydependent carrier phase hardware delay for receiver and satellite, respectively; e , e are relevant system noise and un-modeled residual errors for the un-differenced ionospherefree combination of the code and carrier-phase measurements, respectively; and λ are the wavelengths for un-differenced ionosphere-free carrier frequencies; N is un-differenced ionosphere-free ambiguity bias; c is the speed of light in vacuum; and ρ is the true geometric range from the antenna phase centre of the receiver at reception time to the antenna phase centre of the satellite at transmission time.
Knowing the precise GPS Satellites ephemeris, the outputs of position and velocity from the INS mechanization are used to predict the pseudorange, phase and Doppler measurements through the non-linear observation equations.The corrected pseudorange, phase and Doppler measurements from GPS are differenced with the INS-predicted measurements.Then the integration filter directly processes those residuals to estimate the INS state vector.Finally, the obtained INS state estimates are feed backed to the INS mechanization using the closed loop approach.Figure 1 shows the integration mechanism.

Estimation Filtering
Nonlinear estimation filtering techniques are employed to estimate the proposed integrated PPP-GPS/INS system state vector (described above).In this section, the filter algorithms of the UKF and the PF are briefly reviewed.The general integrated system can be considered as; Where f and h are the motion and measurement nonlinear models, w and v are the motion and measurement noises respectively.

Unscented Kalman Filter (UKF)
In UKF, a set of scaled sigma points with appropriate weights is deterministically chosen so as to capture the mean and covariance of this random vector up to a third order accuracy.Consider a random vector x with mean x and covariance matrix P .The scaled sigma points and the corresponding weights can be defined according to Bergman, (2001) as follows: Initialize with (k=0); Where i=1:n, are the sigma points and n is the dimension of the state vector.The parameter λ is a scaling parameter.Time update step; This contribution has been peer-reviewed.doi:10.5194/isprsarchives-XL-2-127-2014 where ̅ and are the initial state vector and variancecovariance matrix respectively.and are the state and observation vectors for the corresponding sigma points, f and h are the non-linear motion and observation models respectively., , , and , are the time prediction state vector, observation vector and variance-covariance matrix respectively., and are the time update state vector variance-covariance matrix respectively.

Particle Filter PF
In the PF, the probability density function is described and evaluated by the set of samples called particles of x as the numerical approximation based on the Monte Carlo method (Doucet et al., 2000).The filter recursively updates the particles and corresponding weights of the particles as follows; Initialize with (k=0); For i =1 …N, the filter particles are drawn for x from prior P(x .In this research, x and P(x are state vector and variance-covariance matrix solution output of EKF algorithm. Importance sampling (k=1: ∞);

Time update
The prior probabilistic motion density is used as an importance density by passing the state vector samples through the nonlinear mechanization equations

Measurement update
In measurements updating step, the time updating samples are passing through the non-linear measurements system to create the observation probability density; For i =1 …N, the importance weight is evaluated as follow; Normalize the importance weights; Estimation the mean state vector;

Resampling step
In this step, the samples with high weights are selected and redistributed with neglecting the low weights samples.The multinomial distribution resampling technique is applied as pointed out in (Sanjeev Arulampalam et al., 2002).The positioning results of integrated system show that decimeter-level accuracy for all filters while few centimeters difference between the filters estimation as illustrated in Figure 4 for altitude, as an example.The results clarify that the effect of non-linearity in positioning determination can be marginal.In addition there are no consistent improvements in the positioning when non-linear filters are used which indicates that using EKF for high accurate integrated system is sufficient.

Trajectory test and Result analysis
(   For the velocity determination, Figure 5 shows the velocities errors in east, north and up using EKF as a central filter.The results show centimeter-level accuracy can be achieved using single receiver compared with the differential mode.While as shown in Figure 6 millimeter-level difference is produced in east component as an example when using non-linear filters compared with EKF.This small difference comes from the high accuracy achieved originally using EKF.For attitude estimation, due to the lack of determination of the attitude using GPS with multi-antenna, the attitude estimation accuracy will depend mainly on the velocities estimation especially in roll and pitch due to the their strong coupling with the horizontal velocities (east and north), and the quality of the gyros used especially for the azimuth determination.While the velocities estimation using the integrated system as shown before achieved centimeter-level difference with the differential mode integrated system, the attitude determination for the integrated system achieved comparable results with the differential mode integrated system as shown in Figure 7.The same criteria are applied for the estimation of the attitude using non-linear filters compared with EKF.As small velocity estimation difference experienced using non-linear filters compared with EKF, the expected difference between the attitude estimation is small as shown in Figure 8 for roll angle as an example.

Conclusion
This paper investigated the performance of two non-linear filters, namely PF and UKF, in comparison with the traditional EKF for the proposed GPS-PPP/MEMS-IMU integrated system.Un-differenced ionosphere-free linear combinations of carrier-phase and code measurements were processed.Tropospheric delay, satellite clock, ocean loading, Earth tide, carrier-phase windup, relativity, and satellite and receiver antenna phase-center variations were accounted for using rigorous modeling.Tightly coupled mechanism was adopted, which was carried out in the raw measurements domain.Field experiments were executed, which cover different scenarios in land-based navigation.
The test results indicate that insignificant position, velocity and attitude differences (centimeter-level position differences, millimeters-level velocity differences and less than 0.02 degree differences in orientations) are obtained when the measurement updates from GPS carrier phase measurements are available.The future work will include the investigation of the obtainable accuracy of the proposed estimation filters during GPS partial outages.
Vehicular test is conducted to evaluate the performance of the developed integrated GPS-PPP/MEMS-based IMU system.The vehicular test was carried out in downtown Kingston (Canada) and shows very difficult scenarios for satellite navigation, with frequent partial GPS outages of several seconds.The equipments used are Novatel SPAN-CPT system and Trimble R10 receiver.The SPAN-CPT system consists of a Novatel OEM4 receiver and an MEMS-IMU, which contains of three MEMS-based accelerometers and three fiber optic gyros.Carrier phase-based GPS/MEMS-based IMU solution is used to provide the reference solution.In order to create this reference solution, a Trimble R7 receiver was setup on a bench mark station to act as the base station.The raw GPS pseudorange, carrier phase and Doppler measurements in dual-frequency were logged at 1 Hz for the GPS receivers while the IMU raw data is logged at 100 Hz.The duration of the trajectory test was set for about 55 minutes.Figure2shows the trajectory test area.

Figure 2 .
Figure 2. Trajectory areaFigure3shows the positioning errors for latitude, longitude and altitude compared with the carrier phase-based GPS/MEMSbased IMU solution.The positioning results of integrated system show that decimeter-level accuracy for all filters while few centimeters difference between the filters estimation as illustrated in Figure4for altitude, as an example.The results clarify that the effect of non-linearity in positioning determination can be marginal.In addition there are no consistent improvements in the positioning when non-linear filters are used which indicates that using EKF for high accurate integrated system is sufficient.

Figure 4 .
Figure 4. Altitude estimation results difference from EKF results

Figure 6 .
Figure 6.East velocity estimation results for PF and UKF differenced from EKF results