Real-time Kinematic Positioning of INS Tightly Aided Multi-GNSS Ionospheric Constrained PPP

Real-time Precise Point Positioning (PPP) technique is being widely applied for providing precise positioning services with the significant improvement on satellite precise products accuracy. With the rapid development of the multi-constellation Global Navigation Satellite Systems (multi-GNSS), currently, about 80 navigation satellites are operational in orbit. Obviously, PPP performance is dramatically improved with all satellites compared to that of GPS-only PPP. However, the performance of PPP could be evidently affected by unexpected and unavoidable severe observing environments, especially in the dynamic applications. Consequently, we apply Inertial Navigation System (INS) to the Ionospheric-Constrained (IC) PPP to overcome such drawbacks. The INS tightly aided multi-GNSS IC-PPP model can make full use of GNSS and INS observations to improve the PPP performance in terms of accuracy, availability, continuity, and convergence speed. Then, a set of airborne data is analyzed to evaluate and validate the improvement of multi-GNSS and INS on the performance of IC-PPP.


Method
In the INS tightly aided IC-PPP model, the GNSS raw observations (pseudo-range, carrier-phase, and Doppler) and the INS solutions (positions, velocities, and attitudes) are combined in one filter 22,23 . Simply, the observation model can be expressed as where Z k is the innovation vector at epoch k calculated by making difference between GNSS raw observations (O GNSS ) and INS predicted GNSS values (C INS ); H k denotes the designed matrix of the state parameter vector (X k );  k represents the observation noise with the apriori covariance R k . Generally, Z k for INS aided IC-PPP consists of the innovations of pseudo-range (P), carrier-phase (L), and Doppler (D) from two frequencies (f = 1,2), and the other virtual observations including ionospheric delays (I), zenith tropospheric delay (T), and receiver DCBs (d), which can be written as follows    relativity effect, earth rotation effect, antenna Phase Center Offset (PCO) and Variation (PCV) of GNSS receiver and satellite etc. 11 , all of these errors will be corrected by the classic models 11 ;  r denotes the position lever arm correction between GNSS receiver phase center and IMU center; d r j s , and d j s are DCB [14][15][16]23 for GNSS receiver and satellite on frequency j.
Analogously, the innovation functions for carrier-phase on frequency f 1 and f 2 can be expressed as where λ and N represent carrier-phase wavelength and integer ambiguity, respectively; u r j s , and u j s are un-calibrated phase delay 14 of receiver and satellite on j frequency which can be absorbed by float ambiguity ) in PPP parameter estimation, and the random constant model is used to describe the float ambiguities when no cycle slip happens; ∆ρ L r s , stands for other corrections in carrier-phase. In this paper, the ionospheric delays at f 1 frequency (I s 1 ) in line of sight is estimated as parameters 15,16 instead of ionosphere-free combination adopted by convention PPP to eliminate the first order ionospheric delay due to its frequency dependent character, and the ionospheric delays at f 2 frequency can be expressed as The ionospheric delay can be modeled as where q I 2 is the Power Spectral Density (PSD) of the ionospheric dynamic noise (ε k−1 ); Δ t is the interval between two adjacent epochs.
Since there are too many parameters that should be estimated in IC-PPP mode and there are high correlation between the ionospheric delays and the DCB 16 , the PPP solutions are rather weak. To enhance IC-PPP solution, Global Ionosphere Mapping (GIM) data provided by IGS are adopted as virtual observations to constrain the ionospheric delay for each satellite-user pair. The mathematic model can be expressed as 16 : is the ionospheric delay calculated by GIM model with a priori variance σ GIM 2 for its noise (ε GIM ); VTEC and z denote vertical total electron content and zenith angle at ionosphere pierce point (IPP). Considering the temporal and spatial characteristics of ionospheric delay, σ GIM 2 can be further written as cos ( ) cos 14 12 /sin ( ), are the variance of the ionospheric zenith delay and the variance of ionospheric delay variation along with latitude and local time.
As a result, the innovation function for ionospheric delays can be written as is a predicted value in INS aided IC-PPP. Besides, according to the previous research about PPP, the DCBs can be absorbed by receiver clock and satellite clock in convention PPP mode 16 . However, such DCBs will not be absorbed by the single clock in IC-PPP mode, which will impact obviously on the initial position accuracy and the convergence time of IC-PPP 16 . Therefore, the satellite DCBs are eliminated using International GNSS Service (IGS) DCB products and the receiver DCB is estimated as parameter by modeling as random walk process in IC-PPP mode. The model for receiver DCB can be expressed as where q d 2 is the PSD of the DCB dynamic noise (ε k−1 ). And its innovation function can be expressed as d GNSS The tropospheric delay consists of dry component and wet component, and both of the two can be expressed as a product of Zenith Tropospheric Delay (ZTD) and mapping function 34 . The dry delay can be modelled accurately and can be corrected precisely by classical model, but the wet one cannot be corrected well. To weaken the effect of the residual of the wet ZTD on IC-PPP performance, it can be modeled as a random process and estimated as parameter in IC-PPP, expressed as is the PSD of the wet ZTD dynamic noise (ε k−1 ). And the innovation function is defined as . The Doppler data can accelerate IMU sensor errors estimation and provide velocity solution for INS tightly aided IC-PPP. The innovation functions for Doppler at frequency f 1 can be written as According to the mathematical functions above, the designed matrix H k can be obtained by differential of eqs (3,4 and 11). The final form of state parameter vector X k for the INS tightly aided IC-PPP can be written as where δp n , δv n , and δθ are the correction vectors of position, velocity, and attitude, respectively; n denotes local frame (North-East-Down); δB and δS are the correction vectors of biases and scale factors of IMU sensors; Then, the extend Kalmam filter (EKF) is utilized for parameters estimation, and the corresponding state function can be expressed as  where ω en n is the rotation angular rate of navigation frame (n-frame, North-East-Down) with respect to earth centered earth frame (e-frame) projected in n-frame; '× ' denotes cross-product; ω ie n is the earth rotation rate of e-frame with respect to inertial frame (i-frame) projected in n-frame; C b n represents the attitude direction cosine matrix for transforming from body frame (b-frame, Forward-Right-Down) to n-frame; δf b , ω δ ib b , and δg n are offset of accelerometer's specific force, offset of gyroscope's angular rate, and gravity offset.
Although there are many hardware errors contained in INS outputs 35 , just biases and scale factors of gyroscopes and accelerometers should be considered for tactical IMU. Usually they can be modelled as 1 st Gauss-Markov procedure 36 .
where τ is the correlation time; σ is the variance of the driving noise (ξ) depending on the instability of the IMU biases and scale factors. Meanwhile, for one multi-GNSS receiver there exit Inter System Biases (ISB) between each pair of two individual GNSS systems due to different frequencies and signal structures adopted by each system. Besides, since the Frequency Division Multiple Access (FDMA) technique is used by GLONASS, it leads to the Inter Frequency Bias (IFB) for every pair of two GLONASS satellites. In order to avoid such errors' impact on IC-PPP, we estimate 3 independent receiver clock offsets (t t t r G r R r B ) to eliminate the ISBs and 2 + m receiver DCBs to remove the IFBs . To estimate the receiver clock offset and drift accurately, the model proposed by Brown (1992) is applied, which can be expressed as where h 0 , h 2 and K are the instabilities of the GNSS receiver crystal oscillator and the experiential amplification factor. The dynamic behavior of the ionosphereic delay variation, receiver DCBs, and the residual of the wet component ZTD can be described by random walk procedure as shown by eqs (5,9 and 10), respectively. Then, the state transfer matrix can be obtained from eqs (5,9,10,14~16). Finally, according to the mathematical algorithm as described above, the closed loop extend Kalman filter 36 can be employed to update the state parameters and the corresponding variance in the INS tightly aided IC-PPP.

Results
As is well known, the GNSS and INS are usually used to provide the related location information including position, velocity, and orientation for many dynamic applications such as aerial photogrammetry and mobile mapping system. In order to evaluate the performance of the INS tightly aided single-and multi-GNSS IC-PPP in real-time applications, a set of airborne GNSS/INS data from an aerial photogrammetry mission collected in Taiyuan, China on 25 April, 2015 was analyzed. In this experiment, the Trimble R9 multi-GNSS receiver was used to collect GPS, GLONASS, and BeiDou observations (pseudo-range, carrier-phase, and Doppler), and a tactical IMU (POS310) from Wuhan MAP Space Time Navigation Technology Company was used to provide INS data (increments of velocity and angular). The Trimble R9 is a geodetic multi-GNSS receiver with fixing its sampling rate to 1 Hz and the POS310 consists of three fiber-optic gyroscopes and three servo accelerometers with setting the sampling rate to 200 Hz in this test. The whole mission takes about 3.5 hours, and its trajectory is shown in left sub-figure of Fig. 1 with the distances about 41.5 km along North-South direction and 26.9 km in East-West direction, respectively. Shown in the right sub-figure of Fig. 1 is the height of the whole test, and it takes about 2.7 hours in the mission region with the flight height about 1.8 km. The sky plots of available GNSS satellites are displayed in Fig. 2 with the available GPS satellites (2a), the available BeiDou satellites (2b), the available GLONASS satellites (2c), and the three-system available satellites (2d). It is evident that the satellites-users spatial structure in East-West direction is better than that of in North-South direction, which may impact on the positioning accuracy of such two components.
The platform route is arranged into 'S' shape with the main directions along West-North to East-South which are shown visibly in Fig. 1. When the airborne-platform goes to the next route from current route periodically, the velocity and the attitude will change uniformly as shown in Figs 3 and 4, respectively. Accordingly, the velocities are about ± 55 m/s in horizontal and ± 5 m/s in vertical. Meanwhile, the velocity in north and east components (Fig. 3) and attitudes in roll, pitch, and heading directions (Fig. 4) change periodically and significantly due to the 'S' route (moved along West-North to East-South). Here, the changes in roll component are the bank angles with the values about ± 30 degrees for this mission.
In the data evaluation phase, the GNSS/INS data is processed in the IC-PPP mode and the INS tightly aided IC-PPP mode by using raw pseudo-range and carrier-phase of dual-frequency GNSS data with the ionosphere and receiver DCB constraint. Both of the two stages are processed in kinematic mode with the platform   Fig. 6. Accordingly, the averaged number of available satellites and the averaged PDOP are 8.2 and 1.9 while using GPS only. The values are enhanced to 14.3 and 1.4 while using GPS and GLONASS together, and such values are further ameliorated to 22.5 and 1.2, respectively while using the three-system data. It is significant that there are visible improvements in both available satellites and PDOP when multi-GNSS data are used. Besides, we also evaluate the impacts of the number of satellites and PDOP on IC-PPP performance by setting the satellite cut-off elevation angles from 10° to 35° with a stepsize of 5°. According to the results as shown by Figs 7 and 8, it is clear that fewer available satellite number (dropping from 22.5 to 14.2) and larger PDOP (increasing from 1.2 to 3.1) will be obtained along with the satellite cut-off angles increasing. It seems while setting the satellite cut-off elevation angle bigger than 25° the PDOP of G + R + B (Fig. 7) is larger than that of GPS only as shown in Fig. 6, even the number of available satellites is about 2 times larger than that of the GPS only. That is obviously due to the fact that PDOP is not only influenced by the number of available satellites but also constrained by the spatial geometric distribution of these satellites.
In order to validate the performance of the real time kinematic IC-PPP and INS tightly aided IC-PPP, the data is processed in GPS (G), GPS + GLONASS (G + R), and GPS + GLONAS + BeiDou (G + R + B) IC-PPP mode and the corresponding INS aided IC-PPP mode. In the IC-PPP mode, the Global Ionosphere Mapping     improvements when using GPS and GLONASS data together compared to GPS-only solutions. Besides, it seems the solutions of G + R + B IC-PPP are slight worse than these of G + R IC-PPP in north component. It may be caused by the special constellation of BeiDou satellite system at present (5 GEOs over the equator along East-West direction). Generally, the reason for the visible improvements of IC-PPP is that the multi-GNSS can provide more available GNSS satellites and better PDOP. Shown in Fig. 11 are the position differences between the INS tightly aided real time kinematic IC-PPP solutions and the references values in North-East-Up system, and the corresponding statistics are listed in Fig. 10.    As is well known, the number of available GPS satellites in dynamic applications may decrease seriously due to the severe user environments, which will degrade GPS position performance. To analyze the multi-GNSS performance in such complex conditions, we did simulations by setting different cut-off elevation angles in the data processing. The varieties of the number of available satellites and PDOP along with different satellite cut-off elevation angles have been assessed in Figs 7 and 8. In order to evaluate its impact on the position accuracy of real time kinematic multi-GNSS PPP, the experimental data are also processed in G + R + B IC-PPP mode and the corresponding INS tightly aided IC-PPP mode under different cut-off elevation angles ranging from 10° to 35° in a stepsize of 5°, and the corresponding position time series are depicted in Figs 13 and 14, respectively. Evidently, we see that the position RMSs of IC-PPP mode (Fig. 15) decrease dramatically from 4.9 cm, 4.5 cm, and 15.5 cm to 5.9 cm, 11.2 cm, and 33.2 cm in North, East, and Up components, respectively when the satellite cut-off elevation angles increase from 10° to 35°, especially in the vertical direction. In contrast, the corresponding solutions of the INS tightly aided IC-PPP exhibit much better accuracy in both RMS and stability. Compared to the RMSs of G + R + B IC-PPP mode, the average improvements of the position RMS using the INS tightly aided IC-PPP are about 6.9%, 19.3%, and 51.0% in North, East, and Up directions. And the improvement in vertical is much more visible than in horizontal. The RMSs of the INS tightly aided G + R + B IC-PPP are dropped from 4.6 cm, 3.6 cm, and 9.7 cm to 5.9 cm, 7.5 cm, and 16.8 cm along with the satellite cut-off elevation angles increasing from 10° to 35° (Fig. 16). What is more, the RMS of the INS aided G + R + B IC-PPP using 35° cut-off elevation angle is better than that of the GPS-only IC-PPP with 10° cut-off elevation angle (see in Fig. 11). Such simulation indicates that it is possible to obtain high accuracy position in complex environment by multi-GNSS (G + R + B) IC-PPP     low elevation angle will be rejected and fewer GNSS satellites are used in positioning calculation. It can shorten the data processing time with little position accuracy loss as proved in this paper.
In order to get the relative relationship between PPP and the INS tightly aided PPP (both IC-PPP and LC-PPP), we made the difference between PPP solutions under different satellite cut-off elevation angles and that of the INS aided PPP, and the corresponding position differences are also transformed into local coordinate system and the corresponding statistics are listed in Table 1    East directions. It means that INS really can enhance the position accuracy of both IC-PPP and LC-PPP, and such enhancements are more obvious in Up component than in horizontal direction. It may be due to the fact that the availability and sensibility of IMU are much better in vertical direction than in other two directions. Besides, as shown in Figs 16 and 18, the position accuracy of the INS tightly aided PPP degrades slightly along with the increasing cut-off angle, which is much different from that of PPP as shown in Fig. 15. It indicates that INS can provide effective improvement for PPP especially when PPP works not well in worse situations (less available satellites number). Such advantage is very important for the real time dynamic precise positioning applications.
Convergence time is also one of the key factors for real-time IC-PPP. According to the research from Bisnath and Gao (2009) 37 , INS is one of the key tools to shorten the PPP convergence time. Gao et al. (2015Gao et al. ( , 2016 23    For some applications (mapping, navigation, unmanned aerial vehicles etc.), the attitudes determination is also very important. Since no higher precision IMU sensor was used in the test, we just have an assessment about the relative accuracy of the attitudes determination of the INS tightly aided IC-PPP using single-and multi-GNSS data under different satellite cut-off elevation angles. The results are evaluated by comparing with the reference attitudes obtained from the loose integration of RTK and INS. The attitudes accuracy of the RTK/INS (tactical IMU) loosely integration has been evaluated in http://www.novatel.com/assets/Documents/Papers/FSAS.pdf. According to the statistics of attitude differences as depicted in Fig. 22, the RMSs of the attitudes calculated by the eight schemes are extremely close to each other with the maximum discrepancy within ± 0.0006°, and the average RMSs of the attitudes are 0.017°, 0.016°, and 0.104° in roll, pitch, and heading directions, respectively. Besides, the accuracy in heading seems worse than the other two, which may be owing to the worse availability of the gyroscope in heading component.

Discussion
In order to enhance the performance of real-time PPP, we integrate tightly the multi constellation GNSS IC-PPP with INS to take full advantages of the current available GNSS and INS observations. With a set of GNSS and INS dynamic data from an airborne mission, we evaluate the number of available satellites and PDOP of GPS only, GPS + GLONASS, and GPS + GLONASS + BeiDou under different cut-off elevation angles. The multi-GNSS can significantly increase the number of available satellites and improve the spatial geometry in term of PDOP. Then, the performances of the real-time single-and multi-GNSS IC-PPP and the corresponding INS tightly aided IC-PPP are evaluated. The results indicate that both the multi-GNSS and INS can improve the performance of IC-PPP. In details, the multi-GNSS not only make PPP position time series much more stable, but also provide better positions with average improvements of 20%, 50%, and 30% in North, East, and Up components compared to GPS only solutions. The INS tightly aided IC-PPP mode, can further provide another 6% to 40% position improvements in term of RMS for horizontal and vertical components.
Meanwhile, the position accuracy of both GPS + GLONASS + BeiDou IC-PPP and the corresponding INS tightly aided IC-PPP shows a trend decreasing along with the increasing satellite cut-off elevation angles, especially for IC-PPP mode. The position RMS of the multi-GNSS IC-PPP drops visibly in vertical direction when the cut-off elevation angle increases from 10° up to 35°, but the position accuracy loss in horizontal is hardly to see. In comparison, the position RMS of the INS aided multi-GNSS IC-PPP is about centimeter level even setting the satellite cut-off elevation angle to 35°, and no more than 5 cm position accuracy is lost for all of the three components