Low-Cost, Triple-Frequency, Multi-GNSS PPP and MEMS IMU Integration for Continuous Navigation in Simulated Urban Environments

In this research, a next-generation, low-cost triple-frequency GNSS, microelec - tromechanical (MEMS) based inertial measurement unit (IMU)

In the past decade, many manufacturers such as u-blox, Broadcom, MediaTek, Qualcomm (European GNSS Agency, 2018;Wang, 2018), and others have retailed single-frequency (SF) GNSS receivers. The accuracy produced by such low-cost SF GNSS systems is at the meter-level. Dual-frequency (DF), as well as recently triple-frequency (TF) GNSS receivers, were limited to high accuracy applications such as agriculture, avionics, crustal deformation, remote sensing, etc., where the accuracy expected is at the centimeter-level. Along with changing trends in the GNSS technology, some manufacturers such as SwiftNav, Unicore communications, u-blox, and others have started showcasing low-cost DF receivers, evolving in the low-cost receiver space. Septentrio and Unicorecomm are the first among the few to release a TF low-cost GNSS receiver. The accuracy achieved by these receivers is at the cm to dm level in Real-time Kinematic (RTK) processing depending on static or kinematic mode.
Precise point positioning (PPP) is a technique wherein a global network of stations are used to broadcast precise satellite clock and orbit products which can be used by standalone receivers (Zumberge et al., 1997). The downside of PPP is that it takes a few minutes to a few tens of minutes for the position solution to converge to the decimeter level. However, there is a lot of attention that is diverting towards PPP-based GNSS navigation, as well as sensor integration with PPP because of the technique's ability to offer decimeter-level accuracy without requiring local reference infrastructure.
Although there are GNSS receivers available in the market that can perform exceptionally well in terms of software processing (positioning solution), the performance of the receiver is largely restricted by the quality of the antenna (Sun et al., 2017) and in turn quality of measurements, the antenna can procure. The low-cost applications targeted by this research require an antenna that can be used as a rover, is more portable, and is less expensive. Therefore, a patch antenna was a natural choice. Patch antennas have been used widely by SF receivers as they closely possess the characteristics to receive L1 band GNSS signals. However, recently companies like Tallysman, Amotech, and others have released DF and TF patch antennas (Panther, 2012) with good multipath rejection capacity compared to the previous generation of low-cost patch antennas.
While GPS/GNSS and sensor integration is a popular topic of research, GNSS PPP and IMU integration is a more recent area of study and application. Work conducted by Elsheikh et al., (2020); Gao, Ge, et al., (2017); and Gao, Zhang, et al., (2017) used high-precision GNSS receivers such as a Trimble NetR9 with MEMS-based IMU with tightly-coupled (TC) integration. The authors claim that the algorithm achieves a decimeter-level accuracy during simulated outages. These studies achieved decimeter-level accuracy for short outages such as 3-10 s and sub-meter-to meter-level accuracy for longer period outages, depending on the number of satellites used for the solution during the outage. The decimeter-level accuracies are expected because of the high-quality hardware used as the primary sensor for navigation. Elsheikh et al. (2019) studied low-cost, SF-PPP and MEMS-IMU integration using a loosely-coupled (LC) architecture in urban environments in real-time. The authors explicated that the LC algorithm sustained sub-meter level accuracy during short outages such as due to road underpasses. The study also considered a long outage by driving in underground parking for 3 min, where the solution performed with a few meters of error. Low-cost DF PPP + IMU using a geodetic antenna was first analyzed by (Vana et al., 2019) and the authors concluded that the performance of the system remains at a decimeter-level even during 30 s of simulated outages.
Research in the PPP + IMU area explored so far is either related to integrating a high-precision DF GPS/GNSS receiver with a low-cost MEMS IMU or low-cost SF PPP integrated with MEMS-IMU. The only research conducted using a low-cost TF GNSS and MEMS IMU was undertaken by Cheng et al. (2018). The research assesses performance using the BeiDou constellation in RTK mode with a short baseline of 5 km. With the emergence of next-generation low-cost, TF receivers in the market, it becomes an interesting part of GNSS PPP + IMU research to assess the performance of a low-cost, TF GNSS receiver with a low-cost MEMS IMU. DF/ TF GNSS receivers are more attractive for PPP processing than SF ones, as ionosphere refraction (one of the major contributors to GNSS signal error) can be mitigated and DF/TF GNSS PPP position solutions typically converge faster than SF GNSS PPP solutions.
This research work is focused on land vehicle systems that require accurate navigation in the areas of low-cost autonomous vehicles, intelligent transportation systems (ITS), connected vehicles, and others. The Society of Automotive Engineers (SAE) international standard J3016 has published various standards of automation. Currently, manufacturers have vehicles that are partially automated (Level 2) and are now extensively researching and testing Level 3 (conditional automation), where the accuracy required for positioning will be at the decimeter rather than meter to sub-meter level (Dorn et al. 2016;Onishi et al. 2017;Bisnath et al. 2018a). The objective of this research work is to 1) design and develop a low-cost TF GNSS PPP and MEMS IMU integrated solution using TC architecture, 2) apply ionosphere constraining (IC) and land vehicle constraining/non-holonomic constraints (NHC) to the GNSS and IMU measurements, respectively, to reduce the convergence as well as reconvergence period and noise level, and 3) assess the performance of TF GNSS PPP + MEMS IMU solution during a 30 s simulated outage period.
The hardware components used in this research include the low-cost TF GNSS receiver, MEMS IMU, and a TF patch antenna. The low-cost GNSS and IMU are tightly integrated. PPP augmentation and IC are applied on the GNSS measurements while the vehicle constraints/NHC are applied to the MEMS IMU measurements. The novelty addressed by this research is the combination of the low-cost hardware and the software constraining that is used together to provide significantly improved continuous and accurate navigation in the urban environment which has not been examined in the PPP + IMU research area so far. Furthermore, analysis of TF GNSS PPP + MEMS IMU accuracy performance with and without introduced GNSS outages is analyzed via rigorous testing.

DUAL-AND TRIPLE-FREQUENCY PROCESSING
In this research, an uncombined measurement model is used. An uncombined mode of processing is preferred as it allows easy scalability towards more frequencies and allows for the use of external ionospheric information in the form of global ionospheric maps (GIM). The pseudorange and carrier-phase measurements of the three GNSS frequency measurements can be expressed in the form of the following Equation (1) (Hofmann-Wellenhof et al., 2007).
In Equation (1), s represents per satellite, P s i and ϕ s i are the pseudorange and carrier phase measurements, respectively, i = 1, 2, and 3 frequencies, ρ is the geometric range between the satellite and the user position, c is the vacuum speed of light, dt r is the receiver clock errors, dt s is the satellite clock error, T is troposphere delay, γ i = f f i 1 2 2 / is the ratio of frequencies applied to first frequency ionosphere delay I to estimate the ionosphere delay at frequency i, b i s is the satellite hardware bias at frequency i, λ i i c f = / is the signal wavelength, N i is the integer ambiguity on each frequency, IFB r 3, is the inter-frequency pseudorange bias on the third frequency, and e Pi and e i ϕ represent the receiver noise and multipath, respectively, on the measurements with the respective frequencies. The PPP observations are corrected for the relativity effect, Earth rotation effect, satellite antenna phase center offset (PCO) and variation (PCV). The troposphere delay consists of dry and wet delay. The dry component is corrected by using an a priori model and the wet component is estimated as a random walk process. The latter's components are modeled using the global mapping function (Boehm et al., 2006). While performing TF GNSS PPP + IMU integration, if satellite data is available in DF, such a satellite measurement is still used in DF mode in combination with the other satellites that have the TF signals.

IONOSPHERE CONSTRAINING
One of the bottlenecks that causes longer convergence using the PPP technique is the estimation of ionosphere error. To improve the ionosphere estimation during initialization as well as re-convergence, usage of GIM correction has been explored by various research groups such as ; Banville et al. (2014); Cai et al. (2017). The previous research claims to decrease the convergence time by 40-50%. In this research work, the ionospheric parameters have been constrained using GIM products. Typically, in PPP processing, the ionospheric estimates contain a receiver differential code bias (DCB)/code hardware delay lumped into it. In order to be able to perform the constraining with GIM, the receiver DCB needs to be estimated. Therefore, ionosphere delay can be expressed as: In Equation (2), I is the ionosphere delay combined with receiver DCB, I slant is the slant ionosphere delay, and r dcb is the receiver DCB. The constraining has been achieved by forming pseudo-observations between the estimated ionosphere delay and the GIM products and the estimated r dcb . The misclosure for the pseudo-observable is:

INERTIAL MECHANIZATION
The triad of accelerometers and gyroscopes in the IMU package measure specific force and turn rates. The measurement model for accelerometer and gyro measurements with errors is given in Equations (3) and (4) (Farrell, 2008).  f a is the actual accelerometer measurement that accounts for all measurement errors, δ SF a is the accelerometer scale factor error, δ b a is the accelerometer bias, δ nl a is the accelerometer non-linearity error, v a is the measurement noise,  ω ip g is the actual gyro measurement accounting for all measurement errors, δ SF g is the gyro scale factor, δ b g is gyro bias, δ k g represents gyro g-sensitivity, and v g is measurement noise.
 f a and ω ip g  are obtained in the body-frame.
IMU mechanization is the technique of converting the raw IMU measurements into position, velocity, and attitude information. The input to the mechanization process is the body frame specific force and turn rates from the IMU, initial position, velocity, and attitude. The IMU mechanization in the Earth-centered, Earth-fixed (ECEF) frame is pictorially represented in Figure 1. The following steps give a summary of the process (Groves, 2013): 1. Attitude update 2. Conversion of the specific force from the body frame to the ECEF frame 3. Velocity update which includes converting the specific force into acceleration using a gravity model and compensating for the Coriolis effect 4. Position update In this research work, the ECEF frame has been chosen because it is easy to perform the integration for a tightly coupled architecture, where the measurements are in terms of ranges. The mathematical representation of the mechanization process in the ECEF frame is given in Equation (5) (Farrell, 2008).

FIGURE 1 Inertial mechanization in ECEF frame
The dot on the quantity represents the time derivative of the respective quantity. r e is the 3D position, v e is the 3D velocity, and Ω is the skew-symmetric matrix with gyroscope measurements ω.
R Ω and local gravity vector is g G r r e e ie e ie e e = − ( ) Ω Ω . Initial conditions for the initial position ( ( ) r e 0 ), velocity ( ), v e ( ) 0 and attitude (R b e (0)) need to be supplied (Farrell, 2008). The mechanization process in different frames along with more details can be found in, e.g., Farrell (2008); Grewal & Andrews (2001); Groves (2013).

GNSS PPP AND IMU INTEGRATED FILTER
Three methods of integration are mainly followed to integrate GNSS and IMU (Falco et al., 2012;Godha, 2006). The first technique is called loosely-coupled integration, wherein processed position and velocity from GNSS and IMU are used in the misclosure. LC integration involves a smaller number of states. The second technique is tightly-coupled integration, where raw measurements from GNSS (pseudorange and carrier phase) and IMU (specific force/turn rates or the predicted pseudorange/carrier phase) are used in the misclosure of EKF for processing, which means that more states are involved in the estimation process. The third type of integration is deep-coupling, where the IMU is calibrated by the GNSS updates and, in turn, the IMU assists the receiver tracking loops when there is interference or low signal strength . Detailed explanations of these architectures can be found in reference texts (Farrell, 2008;Groves, 2013).
In this paper, the TC approach has been used because this architecture can estimate position and velocity with better accuracy when there are less than four satellite measurements available. Figure 2 is the block diagram of the software architecture used in the research to integrate GNSS PPP with low-cost MEMS IMU. The Septentrio Mosaic-x5 GNSS receiver is used to collect the GNSS measurements and the antenna used in the setup is a Tallysman patch antenna. Precise clock and orbit corrections, along with corrections such as Earth rotation, atmosphere, relativistic, and phase wind up are applied to the measurements. and form the measurement input to the EKF. The outputs from the EKF are corrected IMU position δ P, velocity δ v, attitude δ , bias for accelerometers b a , and bias for turn rates b g . The error states δ δ δ δ δ P v b and b a g , , ,  are fed back to the inertial mechanization to correct for the errors. The corrected position, velocity, and attitude of the IMU is P V and A IMU e IMU e IMU e , , . Based on the knowledge of the car motion, the zero velocity update, zero angular rate update, land vehicle constraint, and height constraining aid in minimizing the unexpected behavior from the IMU measurements while stationary as well as in motion. The details and benefits of the constraints are explained in Vana et al. (2020). Due to the inherent drifting nature of the IMU sensor, constant calibration is necessary.
] is the error in the position represented in the ECEF frame.
] is the error in velocity represented in the ECEF frame. δε δε δε δε = [ ] r p y is the attitude error in roll, pitch, and yaw. d tropo is troposphere wet delay, b a is the accelerometer bias, b g is gyroscope bias, N i is ambiguity for satellite i, δ tc . is GNSS receiver clock drift error, I 1 is the estimated ionosphere per satellite, and r dcb is the receiver DCB per constellation (Park, 2004).
As explained in the previous section, the residuals between the corrected GNSS measurements (pseudorange and carrier phase) and the predicted IMU measurements (pseudorange and carrier phase), and the GIM constrained ionosphere measurements are the inputs to the EKF as represented in Equation (7).
The design matrix A is as shown in Equation (8). It consists of the design model for GNSS-related states.
The first three columns are the partial derivates to the 3D position coordinates with respect to each satellite followed by 1x3 velocity and attitude entries. The next column that has an entry 1 is for clock bias followed by clock drift. The troposphere dry component is modeled and the wet component is estimated according to Boehm et al. (2006) to reduce the impact of errors in the estimation process. λ is a partial derivative for the ambiguity terms Nλ. The last three columns are ionosphere, inter frequency bias (IFB), and receiver DCB related entries in the design matrix.

FIELD TESTS AND RESULTS
Field tests were conducted by collecting kinematic data in the vicinity of the York University campus in Toronto, Canada. The GNSS sensors used in the setup are Septentrio Mosaic-x5 GNSS modules -multi-band, multi-constellation, and triple frequency. The MEMS IMU used in the experiment is the XSens MTi-7, which consists of an SF u-blox GNSS receiver, magnetometer, and a MEMS IMU. The specifications of the XSens MTi-7 are provided in Table 1. Based on these specifications and the specifications in (VectorNav, 2021), one can consider IMU to be an industrial grade MEMS IMU. The antenna used in the setup is a triple-band patch antenna from Tallysman. The patch antenna is capable of covering GPS L1/L2/ L5, GLONASS G1/G2/G3, Galileo E1/E5/E6, and BeiDou B1/B2 and has a good multipath rejection capacity. Figure 3 is a pictorial representation of the setup used for the data collection. The Mosaic-x5 and XSens MTi-7 were placed on the car roof rack along with the Tallysman antenna. In order to have an RTK reference solution, a base station and a rover were set up using Novatel PowerPak7 (enclosure including OEM7 receiver and Epson IMU with low measurement noise) receiver connected to a geodetic antenna. The base station was setup at a known location while the rover high precision system was set up on the car roof beside the low-cost sensors. The frequency of GNSS data collection is 5 Hz, whereas the IMU data were collected at 100 Hz. The raw GNSS and MEMS IMU data were post-processed with the York-PPP+IMU software.
In Figure 4, the experiment set up on the car roof is shown. The box on the car roof rack contains the low-cost sensors as well as the Novatel PowerPak7 (an enclosure including OEM7 receiver and Epson IMU with low measurement noise).
The geodetic antenna for rover reference as well as the low-cost patch antenna were placed on top of the box. The lever arm between the antenna and the IMUs were measured accurately. Data were collected on January 02, 2021, DOY 002. Constellations used for processing are GPS, GLONASS, Galileo, and BeiDou. At the times of data collection, there were no BeiDou-2 satellites visible. An elevation have been used for processing. Data collection was performed for about 40 min on surface roads and highways with some occasional stops at traffic signals in a suburban environment. The distance between the base station and rover was maintained at less than 6 km. Figure 5 shows the vehicle track. Novatel's Inertial Explorer software was used to post-process the data in RTK smoothed mode to provide a reference solution.

Performance -No Outages
Using the data collected on Jan 2, 2021, DF and TF GNSS PPP with MEMS integration has been assessed. First, the integrated solution with single, dual, and triple The horizontal error plots comparing DF PPP + IMU and TF PPP + IMU with and without IC are shown in Figure 6. The convergence threshold assumed in this case is the time taken for the horizontal rms to reach an accuracy of 10 cm. The car was stationary for the initial 15 min and then started moving. Table 2 gives an indication of the horizontal rms of the integrated solutions. The benefit of IC is a reduction in the initial rms in both the DF and TF PPP integrated with IMU. The total rms of DF PPP + IMU and TF PPP +IMU with and without IC are different only by a couple of cm. However, the convergence time taken by the DF and TF GNSS PPP + IMU integrated solutions is slightly longer by a few minutes when IC has not been applied. From Table 2, the accuracy of SF PPP + IMU + IC is 1.2 m which is many-fold higher than the required accuracy targeted by this research as well as DF/TF PPP +IMU. The SF PPP + IMU + IC does not converge to 10 cm accuracy; therefore, the convergence period has not been mentioned. NHC has been applied to IMU measurements for all the solution combinations mentioned in Table 2.
From Figure 6, it can be observed that the initial error reduces significantly in the case of DF PPP + IMU by applying IC when compared to TF PPP + IMU by applying IC. The reason for this behavior is that processing TF GNSS has an increased number of measurements (one pseudorange and one carrier-phase for the third frequency per available satellite) compared to dual-frequency processing. The increased number of measurements in turn increases the estimation degrees of freedom and aids in a more precise estimation of ionosphere states, as well as faster solution convergence. Therefore, after applying IC to TF measurements, the impact of initial horizontal error improvement is not as significant as when compared to applying IC to DF GNSS measurements as the TF IC plays an important role only during initialization to quicken the convergence and to decrease initial position error. In the initial period of convergence, GIM values support estimation of the ionosphere states. However, GIM corrections are not very accurate, and they cannot be relied on for precise ionospheric estimation. Given that GIM corrections are only accurate to the same noise level of pseudorange measurements Banville et al., 2014;Cai et al., 2017), a similar weight is given to both pseurodange measurements and GIM corrections. As a consequence, similar to the pseudorange measurements, GIM corrections only have an effect on the PPP solution during the initial stages of the processing. Once the ambiguities and ionosphere are estimated more precisely, phase measurements are given higher weight which reduces the effect of GIM corrections on ionospheric estimation and position estimates. Figure 7 illustrates the availability of the total number of satellites (as well as per constellation) during the data collection period. The average number of satellites during the data collection was 18.
The accuracy of a GNSS PPP + IMU integrated solution depends on the quality of IMU used. A low-cost MEMS IMU was used in this work and the performance  (6), the number of unknowns are: 3D position states, one troposphere wet delay component, one receiver clock, one receiver clock drift, IFB, receiver DCB, three ambiguity states, and ionosphere state per carrier phase frequency per TF satellite. Therefore, the number of unknowns can be given as: In Equation (9), N is the number of satellites. The number of satellite-dependent states are 4 per satellite. At least 4 satellites are required to estimate all the unknowns in a given epoch for TF processing without any a priori information. While processing dual-frequency the number of unknowns is given by: In Equation (10), N is the number of satellites. The number of satellite-dependent states are three per satellite. The number of states that are not dependent on the satellite is reduced by one in the DF case, because the associated IFB need not be estimated. During DF processing, at least 7 satellites are required. While three frequency measurements are not available for a satellite, DF processing has been applied. In Figure 6, the horizontal error increases between 22 and 22.1 hr for a brief period of about a minute while the vehicle makes a turn to exit from a parking lot. The data were collected in an open sky parking lot and the exit is surrounded by trees towards the north and south directions, causing fluctuation in the number of satellites as seen in Figure 7. The sky view distribution of the satellites during the period is shown in Figure 8. The distribution of satellites just before exiting the parking lot is shown in Figure 8(a). While exiting the parking lot, the trees cause a poor satellite geometry for a few seconds as shown in Figure 8(b). Figure 8(b) has FIGURE 7 Plot of satellites available been highlighted in red to show that all satellites except one lie in the east-west direction. The brief period of poor geometry causes the horizontal position error to increase, and it takes a few tens of seconds to converge back to 10 cm accuracy.

Performance -Simulated Outages
The primary advantage of integrating an IMU with GNSS is that the IMU sensor aids in providing an accurate and continuous navigation solution during partial or complete GNSS outages. In this section, first a case study of introduced outages with varying number of satellites (0 to 4) during a single simulated outage is studied and then results from a rigorous study of multiple simulated outages is presented with varying number of satellites (0 to 4) in several portions of the data. Figure 9 is a plot of horizontal accuracy of position when a 30-s outage is introduced at 22.14 hours. DF PPP + IMU and TF PPP + IMU solutions are analyzed in Figure 9 to compare the performance when the number of satellites available during the outage is varied from 4 down to 0. Statistics for the horizontal accuracy are summarized in Table 3.
The performance of both DF PPP + IMU and TF PPP + IMU is very close when there are 4 satellites available during the outage, as 4 satellites are still enough to maintain a reasonable decimeter-level PPP solution for half a minute. However, as the number of available satellites decreases from 4 down to 0, the GNSS update starts to negatively impact the solution during the outage, and the position error during and after the outage grows as the number of satellites decreases. However, the addition of the third frequency improves overall performance, as the overall rms of the integrated solution with the third frequency is 5-20% better than DF PPP+IMU.
From Table 3 and Figure 9, it can be noticed that overall rms for the 2-satellite case is better for DF PPP + IMU by 10 cm as compared to TF PPP + IMU. However, the rms after the outage in the 2-satellite case is better in the case of TF PPP + IMU by 6 cm as compared to DF PPP + IMU. Such behavior prevents a conclusive remark on the overall accuracy between DF PPP + IMU and TF PPP + IMU just by looking at one case study, so multiple simulations were conducted in different areas of the dataset. Figure 10(a) depicts a plot of the overall rms accuracy of the data involving the DF PPP + IMU and TF PPP + IMU during the 30 s of introduced outages. Several simulations including a combination of different constellations and frequencies (DF and TF) were carried out for different categories of satellite availability (0-4) Time series of horizontal error to depict the performance of DF PPP + IMU and TF PPP + IMU when the number of satellites is varied from 4 down to 0 during an outage at 21.14 hours and the results were averaged. NHC has been applied to the IMU measurements for all the solution combinations analyzed. The represented plot in Figure 10 are averaged values of horizontal rms from the multiple simulations. From Figure 10(a) it is evident that the TF GNSS PPP + IMU performs better than the DF GNSS PPP + IMU as the number of satellites decreases. When the number of satellites available during the outage were 4, the horizontal rms accuracy remained close to a decimeter with DF PPP + IMU and by adding the third frequency, the accuracy goes down to cm level. When the number of satellites decreases to 1, the TF GNSS PPP + IMU solution is better than the DF GNSS PPP + IMU solution by about 3 dm. It is also very interesting to notice that the overall rms while using a TF PPP + IMU with 2 and 3 satellites available during the introduced outages impressively remains within 2-3 dm. Several decimeters of horizontal accuracy using low-cost equipment is valuable and significant for low-cost applications. Although one would expect that the performance of the navigation system would be the same when there are no satellites available during the complete 30 s of outage period, the results from Figure 10(a) indicate that TF PPP + IMU performs better than DF PPP + IMU at the sub-meter level. The first point to note is that the results in Figure 10 are a representation of overall rms and rms after an outage and not rms during an outage. Secondly, superior TF PPP + IMU performance is due to the additional GNSS measurements that enter into the estimation process and aid in higher precision estimation before and after the outage. Given that the experiments conducted are introduced outages, the DOP between DF PPP +IMU and TF PPP + IMU is the same but the additional measurements from TF PPP + IMU will help in quicker convergence after an outage as the accuracy of estimation of ionosphere is better (Guo et al., 2016;Naciri & Bisnath, 2021;Yi et al., 2021). Once a complete GNSS outage (0 satellites tracked) occurs, the estimation process is re-initialized. In this case, having the third frequency measurements helped in reducing the initial error and convergence time after an outage, resulting in a better rms compared to DF PPP + IMU. The test results in Figure 10(b) show the performance accuracy of the integrated solution after the introduced outage. The rms computed in Figure 10(b) includes the rest of the data set after the introduced outage. In the PPP technique, slow reconvergence of the solution is one of the drawbacks. While the horizontal rms difference between TF GNSS PPP + IMU and DF GNSS PPP + IMU with 4 satellites availability is only different by a few centimeters, the difference is clearly noticeable when there are 3 to 0 satellites available. In all four cases (0 to 3 satellites are available), the TF PPP +IMU is better than DF PPP +IMU by a few decimeters. The re-convergence of the solution with good accuracy after an outage is important as the partial GNSS signal availability in an urban environment is quite common. Therefore, one can infer from the results shown in Figure 10(b) that the higher accuracy during the reconvergence period helps in ensuring that overall rms is maintained at a decimeter level. From Figure 10, TF PPP + IMU performs better than DF PPP + IMU in the case of lower number of satellite availability (less than 4 satellites).
The impressive accuracy performance with the addition of the third frequency is due to the additional measurements that are available before and during the outage which increase the degrees of freedom in the estimation process. The degrees of freedom and the number of satellites required to make the system determined or overdetermined is also explained in Equations (9) and (10). The increased degrees of freedom improve the accuracy of the GNSS as well as IMU states such as the attitude error, accelerometer, and gyroscope biases. An estimation with greater precision of the IMU states before an outage helps the IMU to sustain positioning performance during an outage by keeping the error levels lower, which in-turn, aid in re-convergence after an outage.

Multi-Frequency GNSS PPP + IMU With Ionosphere Constraining Performance -Simulated Outages
The benefits of applying IC have been studied by various groups for different grades of GNSS. As part of this research, the impact of IC on the TF PPP + IMU solution is studied. DF PPP + IMU and SF PPP + IMU are also visualized in this context in order to compare the performance of TF PPP + IMU and to have a full picture of all the three frequency bands available on the market.
In this section, first a case study of introduced outages with two satellites during the outage is presented and then results from a thorough study with multiple simulations are discussed.
As part of the case study, a half-minute partial outage was created with only two satellites available for processing just after 21.14 hours. Figure 11 is a representation of horizontal accuracy performance of SF PPP + IMU with IC, DF PPP +IMU with and without IC, and TF PPP + IMU with and without IC during the partial outage. The SF PPP + IMU + IC solution is accurate at the meter-level just before the outage, but the remaining solutions perform at the decimeter-level of accuracy. Table 4 provides the horizontal positioning accuracy statistics of the case study for all discussed solution combinations, where overall rms and rms after an outage are examined. Rms after an outage is specified as re-convergence plays an important role in PPP processing, which in turn contributes to the overall rms and performance. Table 4, the SF PPP + IMU + IC solution performs with lower accuracy, while the TF PPP + IMU + IC solution has the best accuracy performance in both overall rms, as well as rms after the outage categories. One may expect that applying IC would help in all conditions and situations for initial convergence as well as re-convergence. However, in this case study, the performance of DF PPP + IMU with IC performs with an rms of 70 cm after the outage, while the DF PPP + IMU performs with a 50 cm rms. This behavior prompted more rigorous tests to be conducted to draw a conclusion for this comparative study.
As the case study was conducted with two satellites partial outage in only one area of the dataset, a more rigorous study was conducted involving simulated outages in different portions of the data. The study not only includes creation of outages in different areas of the data, but also the number of satellites available during the outage were varied from 4 down to 0. Figure 12 is a representation of averaged results of multiple simulations to compare the performance of DF PPP + IMU and TF PPP + IMU with and without IC applied to the integrated solution. NHC has been applied to the IMU measurements for all the solution combinations. Figure 12(a) is the overall rms of the integrated solution averaged over various rigorous simulations performed with 30 seconds of introduced outages. Multiple simulations were performed for each of the following test cases: 1) a combination of satellites from different constellations; 2) satellites from the same constellation; and 3) combination of frequencies (DF and TF) available from each satellite for different levels of satellite  Figure 12(a) is the total absence of satellites. During this period, SF PPP +IMU + IC still performs with the least accuracy as expected. However, it is very impressive to note that the low-cost hardware still maintains overall meter-level accuracy. When there are no satellites available, the TF PPP + IMU + IC has 22% improved accuracy when compared to DF PPP + IMU, which is significant for such low-cost equipment. SF PPP + IMU + IC performs with the least accuracy compared to DF/TF PPP + IMU in all cases since mitigation of ionosphere and fast convergence capability of SF GNSS PPP is limited. It is evident from Figure 12 that the accuracy improvement between SF PPP with IMU and TF PPP with IMU is huge, i.e., in terms of meters. Figure  12(b) illustrates the rms of the integrated solution after the introduced outages. The rms computed in Figure 12(b) includes the entire data set after the introduced outage. IC aids in quickening the reconvergence, and hence the analysis of accuracy during reconvergence has been performed. From Figure 12(b) it is evident that the IC does not significantly help in improving the accuracy when the satellites available move from 4 to 3. However, as the number of satellites drop below 3, the accuracy improvement by applying IC to TF PPP + IMU is obvious. For example, when there are two satellites visible, DF PPP + IMU and DF PPP + IMU + IC performance is within cm, however, TF PPP + IMU + IC shows a decimeter improvement over TF PPP + IMU and a couple of decimeters improvement over DF PPP + IMU. As the visible satellites further reduce to 1, TF PPP + IMU + IC performs at a 25% better accuracy than TF PPP + IMU and 34% better accuracy than DF PPP + IMU. While there are no satellites available for the complete 30-s outage period, the performance improvement achieved by TF PPP + IMU+IC is sub-meter better than DF PPP + IMU. The accuracy improvements during reconvergence are promising results for low-cost applications used in urban canyons where the satellite geometry varies very often.
The improvement exhibited by the TF PPP + IMU + IC algorithm over the DF PPP + IMU algorithm is represented in terms of percentages in Figure 13. The comparison between DF PPP + IMU and TF PPP + IMU + IC has been presented as a summary of the detailed results from Figure 10 and Figure 12. Also, Figure 13 is a comparison of the DF PPP + IMU that is being used as a default combination, e.g., the automotive industry and the proposed solution TF PPP + IMU + IC. Overall rms of TF PPP + IMU + IC is always better than the DF PPP + IMU when the satellite availability varies from 4 to 0. The accuracy enhancement by using TF PPP + IMU + IC is 10% when there are 4 satellites. However, the accuracy improvement jumps to 40% as the number of satellites are less than 4. The trend of the rms after an outage by TF PPP + IMU +IC exhibits a betterment within the range of 18-30% when compared to DF PPP + IMU. These results reinforce that the TF PPP + IMU + IC is definitely achieving decimeter-level accuracy with significant improvements over a DF PPP + IMU navigation system. The results explained are significant for the modern low-cost application that demands decimeter-level accuracy even in an urban canyon.

CONCLUSIONS AND FUTURE RESEARCH
Providing a continuous positioning solution in a suburban and urban environment using low-cost sensors has been a challenging problem. The trade-off between the accuracy of positioning and cost of the navigation system has existed for some time until the emergence of a new generation of multi-frequency low-cost GNSS receivers as well as MEMS IMUs entered the market. In this research, a complete FIGURE 13 Plot of percentage improvement between DF GNSS PPP+IMU and TF GNSS PPP+IMU +IC low-cost TF GNSS PPP, MEMS IMU, and patch antenna hardware were used to perform TC integration. As the low-cost sensors possess inherent errors, PPP augmentation, IC, and the third frequency are used to minimize the GNSS errors while land vehicle constraining/NHC are used to calibrate the IMU errors. The unique combination of the next-generation low-cost hardware and software (PPP augmentation, IC, and NHC) provides a novel resolution to the continuous navigation solution in the GNSS-denied environment assuring a decimeter-level accuracy even when the number of satellites available decreases from 4 to 2.
The TF and DF PPP integrated with IMU along with IC and NHC perform at a similar accuracy of less than a decimeter when there is good GNSS availability. 30 s of GNSS outages were introduced in different parts of the dataset and rigorous testing was performed by including multiple test cases involving multiple simulations for each of the following test cases: 1) a combination of satellites from different constellations; 2) satellites from the same constellation; and 3) combination of frequencies (DF and TF) available from each satellite for different levels of satellite availability (ranging from 0 up to 4) and the results were compared. During the outages, the overall rms and rms after the outage indicate that including the third frequency helps improve the integrated solution and reduce the impact of the outage, as the rms decreases with TF integrated solution compared to DF integrated solution. As the number of satellites available decreases, the overall rms of the integrated TF PPP with IMU along with IC and NHC performed better than the DF PPP integrated with IMU by 10-40% while the improvement in rms after the outage is in the range of 18-30%. Applying IC and NHC to TF PPP with IMU helps in improving the accuracy of the solution during reconvergence after an outage when the number satellites available is 0-2. These numbers indicate that including the third frequency on low-cost GNSS chips, combined with constraining the ionosphere, certainly helps in improving the integrated solution when there are fewer satellites available. The superior performance of TF PPP + IMU when compared to the DF PPP + IMU comes from the fact that TF PPP provides additional measurements through the third frequency and aids in improved estimation of the IMU states before, during, and after an outage.
The emergence of low-cost TF-GNSS receivers in the market, along with low-cost MEMS IMUs marks a new generation and beginning for more research in integrated navigation. Usage of lower-cost sensors helps in lowering the cost of modern and next-generation applications but with limited compromise in accuracy, precision, and continuity in solution. The outcomes of this research work are reassuring to achieve a continuous navigation solution with decimeter-level accuracy for low-cost autonomous applications, virtual reality, augmented reality, and other next-generation applications.
Future work involves testing the algorithm in an urban environment where the change in the satellite geometry is quite dynamic. To improve the low-cost solution, the vehicle odometer will be integrated.

a c k n o w l e d g e m e n t s
The authors would like to thank Nacer Naciri from the GNSS Laboratory, York University for helping with numerous trials in data collection and York-PPP software debugging issues. The authors would also like to thank the Natural Sciences and Engineering Research Council (NSERC) and York University for providing funding for this research. The German Research Center for Geosciences (GFZ), International GNSS Services (IGS), and Centre National d'Etudes Spatiales (CNES) are acknowledged for correction data.