Ultra-Low-Cost Tightly Coupled Triple-Constellation GNSS PPP/MEMS-Based INS Integration for Land Vehicular Applications

: The rapid rise of ultra-low-cost dual-frequency GNSS chipsets and micro-electronic-mechanical-system (MEMS) inertial sensors makes it possible to develop low-cost navigation systems, which meet the requirements for many applications, including self-driving cars. This study proposes the use of a dual-frequency u-blox F9P GNSS receiver with xsens MTi670 industrial-grade MEMS IMU to develop an ultra-low-cost tightly coupled (TC) triple-constellation GNSS PPP/INS integrated system for precise land vehicular applications. The performance of the proposed system is assessed through comparison with three different TC GNSS PPP/INS integrated systems. The ﬁrst system uses the Trimble R9s geodetic-grade receiver with the tactical-grade Stim300 IMU, the second system uses the u-blox F9P receiver with the Stim300 IMU, while the third system uses the Trimble R9s receiver with the xsens MTi670 IMU. An improved robust adaptive Kalman ﬁlter is adopted and used in this study due to its ability to reduce the effect of measurement outliers and dynamic model errors on the obtained positioning and attitude accuracy. Real-time precise ephemeris and clock products from the Centre National d’Etudes Spatials (CNES) are used to mitigate the effects of orbital and satellite clock errors. Three land vehicular ﬁeld trials were carried out to assess the performance of the proposed system under both open-sky and challenging environments. It is shown that the tracking capability of the GNSS receiver is the dominant factor that limits the positioning accuracy, while the IMU grade represents the dominant factor for the attitude accuracy. The proposed TC triple-constellation GNSS PPP/INS integrated system achieves sub-meter-level positioning accuracy in both of the north and up directions, while it achieves meter-level positioning accuracy in the east direction. Sub-meter-level positioning accuracy is achieved when the Stim300 IMU is used with the u-blox F9P GNSS receiver. In contrast, decimeter-level positioning accuracy is consistently achieved through TC GNSS PPP/INS integration when a geodetic-grade GNSS receiver is used, regardless of whether a tactical- or an industrial-grade IMU is used. The root mean square (RMS) errors of the proposed system’s attitude are about 0.878 ◦ , 0.804 ◦ , and 2.905 ◦ for the pitch, roll, and azimuth angles, respectively. The RMS errors of the attitude are signiﬁcantly improved to reach about 0.034 ◦ , 0.038 ◦ , and 0.280 ◦ for the pitch, roll, and azimuth angles, respectively, when a tactical-grade IMU is used, regardless of whether a geodetic- or low-cost GNSS receiver is used.


Introduction
GNSS precise point positioning (PPP) is capable of providing precise positioning solutions without any additional GNSS base stations [1]. A costly geodetic grade GNSS device is required for PPP to achieve precise positioning, which limits the use of PPP in a wide range of commercial applications. The release of dual-frequency (DF) GNSS smartphones such as Xiaomi mi 8, which supports dual-frequency (DF) GNSS measurements, allowed users to remove the ionospheric delay, which in turn led to an enhanced PPP solution [2]. It has been shown in [3,4] that DF PPP through smartphones can achieve The IF linear combinations of pseudorange and carrier phase measurements, after accounting for the previously mentioned errors for GPS, GLONASS, and Galileo, can be written as: where G, R, and E refer to the GPS, GLONASS, and Galileo satellite systems, respectively; P G IF , Φ G IF , P R IF , Φ R IF P E IF , and Φ E IF refer to the GPS, GLONASS, and Galileo IF linear combinations; ρ G ,ρ R ,ρ E are the geometric range between the receiver's antenna phase center and the satellite antenna phase center; b G r is the receiver clock error; zwd is the zenith wet delay; m w is the wet mapping function; ISB R and ISB E are the inter-system bias between the GLONASS and Galileo satellite systems, respectively, and the GPS satellite system; and ε PIF and ε ϕIF refer to noise and multipath effect of pseudorange and carrier-phase IF linear combinations, respectively. Differential code bias (DCB), which refers to the difference in signal travel time for two signals of a particular GNSS constellation, is required for dual-frequency users if the used pair of signals to formulate the IF linear combinations are different from the ones used to generate the GNSS precise clock products [23]. The GPS precise clock products are generated using the P-code on L1 (C1W) and L2 (C2W) frequencies [24], while the code measurements on E1 (C1X) and E5a (C5X) are used for Galileo precise clock products [25]. As a result, the GPS C1C pseudorange measurements have to be corrected by DCB P1-C1 to be consistent with the precise clock products when the Trimble R9s receiver is used. On the other hand, the GPS C1C pseudorange measurements have to be corrected by DCB P1-C1, and Galileo C7Q pseudorange measurements have to be corrected by DCB 7Q-5X, when the u-blox-F9P receiver is used. The zenith dry component of the tropospheric delay is accounted for using the Saastamoinen model [26]. The dry and wet mapping functions are determined using the Vienna mapping function (VMF) [27]. The effects of relativity, sagnac delay, phase center offset and variation, Earth tides, ocean loading, and phase wind up are modeled as described in [28]. In the adopted GNSS PPP model, the pre-saved real-time orbit and clock products for GPS, GLONASS, and Galileo are obtained from the Centre National d'Etudes Spatials (CNES) analysis center (available at: http://www.ppp-wizard.net/products/REAL_TIME/, accessed on 20 June 2020).

Full IMU Mechanization
The IMU measurements are provided in the vehicle body frame, in which the x-axis refers to the transversal direction, the y-axis refers to the vehicle's forward direction, and the z-axis completing a right-handed system as adopted in this study. The IMU mechanization outputs are produced in the local-level frame (LLF). The LLF axes point towards east, north, and up (ENU) as adopted in this study [29]. In the LLF frame, the mechanization outputs include the position (latitude ϕ, longitude λ, and altitude h), velocity (east velocity v e , north direction v n , and up velocity v u ), and attitude angles (roll r, pitch p, and yaw y) of the IMU. Firstly, the measured accelerations and angular rotations by the IMU are converted to incremental velocities and incremental angles. Then, the accelerometer and gyro biases are then used to correct the incremental velocities and incremental angles as [8]: where ∆v b is the corrected velocity increment; ∆θ b ib is the corrected angular increment; f b and w b ib are the measured accelerations and angular rotations in the body-frame (b-frame); i refers to the inertial frame; ∆t is the sampling interval; b a is the accelerometer bias; and b g is the gyro bias. The initial values of r and p angles are determined through the initial alignment process using static IMU data [30]. The initial y angle is determined in this study using GNSS observations from two antennas mounted on the roof of the test vehicle. The initial values of the attitude angles are used to set the initial transformation matrix between the b-frame (b) and LLF (l) as follows: cos y cos r − sin y sin p sin r − sin y cos p cos y sin r + sin y sin p cos r sin y cos r + cos y sin p sin r cos y cos p sin y sin r − cos y sin p cos r − cos p sin r sin p cos p cos r   (9) Subsequently, the initial quaternion (q) parameters are estimated as [8]: The incremental angle (∆θ) between the b-frame and LLF can be then determined, after accounting for the Coriolis effect, as [30]: where w l il = −v n M+h w e cos ϕ + v e N+h w e sin ϕ + v e tan ϕ N+h T ; R b l is the transformation matrix from LLF to b-frame, which equals the transpose of R l b ; and N is the prime vertical meridian radius of the curvature of the Earth. The updated quaternion can be then derived using ∆θ b lb as follows [8]: where k and k − 1 are two consecutive epochs; Ω is the skew-symmetric matrix representation of ∆θ b lb and is formulated as follows: The updated transformation matrix is constructed using the updated quaternion parameters as: Using the updated transformation matrix R l b , the attitude angles (r, p, y) are determined as follows: The velocity increment is then corrected for Coriolis and gravity effects using the updated R l b as: where Skew symmetric matrices' representation of the angular rotations between the ECEF (e) and inertial frame (i) Ω l ie and between the LLF and ECEF frame Ω l el can be represented as in [30]. The gravity in the LLF g l can be presented as: The velocity can be updated as: The position update, namely latitude, longitude, and altitude at epoch k using the corresponding values for the previous epoch k − 1, can be presented as follows: where M and N are the meridian and prime vertical meridian radii of curvature of the Earth, respectively.

Tightly Coupled GNSS PPP/INS Integration
In this study, an improved robust adaptive Kalman filter (IRKF) is adopted and used as the estimation filter. Raw GNSS carrier-phase and pseudorange measurements are first obtained from the GNSS receiver and corrected using the error models as described in Section 2. Both of the full mechanization position outputs and the orbit and clock products are then used to predict the carrier phase and pseudorange measurements. The difference between the corrected GNSS measurements and the predicted IMU-based counterparts is fed to the IRKF. The IRKF uses the measurement differences to update the predicted position, velocity, and attitude errors from the system model. The updated errors are used to correct the IMU mechanization outputs to provide the TC GNSS PPP/INS integrated position, velocity, and attitude solutions. The estimated accelerometer and gyro bias errors are then fed back to the IMU mechanization through the closed-loop scheme. Moreover, the estimated receiver clock error, zenith wet delay, ISB, and ambiguity parameters are fed back to the PPP errors correction stage. The proposed TC GNSS PPP/INS integration is implemented as shown in Figure 1.

Tightly Coupled GNSS PPP/INS Integration
In this study, an improved robust adaptive Kalman filter (IRKF) is adopted and used as the estimation filter. Raw GNSS carrier-phase and pseudorange measurements are first obtained from the GNSS receiver and corrected using the error models as described in Section 2. Both of the full mechanization position outputs and the orbit and clock products are then used to predict the carrier phase and pseudorange measurements. The difference between the corrected GNSS measurements and the predicted IMU-based counterparts is fed to the IRKF. The IRKF uses the measurement differences to update the predicted position, velocity, and attitude errors from the system model. The updated errors are used to correct the IMU mechanization outputs to provide the TC GNSS PPP/INS integrated position, velocity, and attitude solutions. The estimated accelerometer and gyro bias errors are then fed back to the IMU mechanization through the closed-loop scheme. Moreover, the estimated receiver clock error, zenith wet delay, ISB, and ambiguity parameters are fed back to the PPP errors correction stage. The proposed TC GNSS PPP/INS integration is implemented as shown in Figure 1.

System Model
The discrete form of the system model is given as [8]:

System Model
The discrete form of the system model is given as [8]: where δx k (−) is the predicted state error vector at epoch k; Φ k−1,k = F(t)∆t + I; Φ k−1,k is the transition matrix; F(t) represents the system dynamic matrix; ∆t is the sampling interval; I is the identity matrix; δx k−1 (+) is the estimated state error vector at epoch k − 1; G is the noise coupling matrix; and w refers to the system noise with Q covariance matrix. For the TC GNSS/INS integration, the state vector is given by: where δr L refer to the position errors; δv L refer to the velocity errors; δε L refer to the attitude angles errors; δb a refer to the accelerometer bias errors; δb g refer to the gyro bias errors; δb G r refer to the receiver clock error; δzwd refer to the wet zenith delay error; δISB refer to the ISB error; and δN IF refer to the ambiguity parameter errors. The dynamic matrix can be adopted as in Equation (A1). The first-order Gauss Markov process is used to model the accelerometers and gyro bias as described in Equations (A10) and (A11) in Appendix A [10]. The noise coupling matrix G is in Equation (A12) in Appendix A [31,32]. In the prediction stage, the system transition matrix as in Equation (22) is used to predict the state vector δx k (−), and its a priori covariance matrix P k (−) can be estimated as: where P k−1 (+) is the state vector posteriori covariance matrix at epoch k − 1. It is worth mentioning that the a priori covariance matrix is assumed to be a diagonal matrix at the initial epoch. The Q matrix is diagonal, and it is set based on the system noise covariances with each state vector parameter. The system noise includes the accelerometer noise vector, gyro noise vector, Gaussian noise associated in the GM process of accelerometer bias, Gaussian noise associated in the GM process of gyro bias, driving noise for receiver clock error, driving noise for zenith wet delay error, and driving noise for ISB. The receiver clock error, zenith wet delay, and ISB are modeled as random walk processes. The stochastic parameters of both accelerometers and gyros provided in the IMU's data sheets are used to construct the initial form of the process noise covariance matrix. However, Q requires further tuning to achieve the required performance for the application under consideration.

Measurements Model
The measurement model for the TC GNSS PPP/INS integration can be written as: where υ is the observations noise with R covariance matrix; H is the coefficient matrix; and δz refers to the measurement error vector. The measurement error consists of the differences between the GPS, GLONASS, and Galileo IF linear combinations of carrier phase and pseudorange measurements and the predicted IMU-based counterparts, and can be written as: Geomatics 2021, 1

265
The coefficient matrix for the TC GNSS PPP/INS integration can be constructed as: where U r is an Nx3 matrix that contains a set of unit vectors u r ; N is the number of tracked satellites for each satellite's system, namely N G , N R , and N E ; u r is a 1 × 3 unit vector between the receiver and the corresponding satellite; and N G , N R , and N E refer to the GPS, GLONASS, and Galileo visible satellites. The u r matrix is multiplied by a transformation matrix [8] to convert the coordinate corrections from the ECEF frame δx, δy, δz to the curvilinear coordinate corrections δϕ, δλ, δh. In the update stage, the Kalman gain K can be written as: In the IRKF, an adaptive factor is adopted to balance the contribution of the predicted state vector and the robust estimation using the new measurements to compensate for the system model errors [33,34]. The adaptive factor can be determined as [35]: where ∆ X k is the discrepancy in the predicted state; c • and c 1 are constants and can be set as c • = 1 → 1.5 ; c 1 = 3 → 4.5 ; and α is the adaptive factor, which can take values in the range 0 to 1. ∆ X k can be determined as [35]: where δ x k is the state vector estimate using the GNSS measurements at epoch k. As shown in Equation (29), the coefficient α takes the value of 1 when the ∆ X k is less than c • , which indicates a good estimate of the predicted state vector. On the other hand, the coefficient α takes the value of 0 when ∆ X k is higher than c 1 , which indicates that the contribution of the predicted state to the update stage has to be minimized, as shown in Equation (31). The values of c • and c 1 need to be further tuned to achieve the required performance for the application under consideration. In this study, the values of 1.5 and 4.5 are adopted for c • and c 1 , respectively. To control the predicted state contribution, a new formula of the Kalman gain is presented as: It is worth mentioning that the measurement covariance matrix R is assumed to be a diagonal matrix. The diagonal elements are the variances of the IF linear combinations of the pseudorange measurements σ P,IF and carrier phase measurements ϕ P,IF , which are given by [36]: where el refers to the elevation angle of the satellite; ; f 1 and f 2 are GPS L1/L2 frequencies, GLONASS L1/L2 frequencies, and Galileo E1/E5a or E1/E5b frequencies; and σ P and σ ϕ refer to the pseudorange and carrier phase standard deviations. In this study, two different GNSS receivers are used, namely Trimble R9s and u-blox F9P. The standard deviation of the carrier phase measurements is scaled by 100 times less than the pseudorange due to its high precision [37]. The adopted pseudorange and carrier phase standard deviations are 0.3 m and 0.003 m when the Trimble GNSS receiver is used. Due to the difference in the quality of the GNSS observations between the geodetic-grade and low-cost GNSS receivers, the adopted pseudorange and carrier phase standard deviations are 3 m and 0.03 m when the u-blox GNSS receiver is used. The a posteriori covariance matrix of the updated state vector can be estimated through [13]: The innovation term dx k , which is used to update the predicted state vector, can be estimated as [10]: In this study, we are using the closed-loop scheme, which means that δx k (−) = 0. Therefore, the innovation term can be re-written as: In the IRKF, a modified covariance matrix for the GNSS measurements is introduced to compensate for the effect of measurement outliers. The modified measurement covariance matrix R, which is assumed to be a diagonal matrix, is constructed for GPS, GLONASS, and Galileo pseudorange and carrier phase measurements as [38]: The measurements' residuals are firstly estimated to construct R as: The standardized residuals are then estimated as: The standardized residuals are then used to determine the robust classification factor P for each specific measurement type individually. The classification factor for the ith Geomatics 2021, 1 267 number of measurements in the jth type of measurements P j i can be estimated using the t-test statistics as [38]: where T j i is the t-test statistic for the ith measurement in the jth type as described in [38]; t • and t 1 are the t-test critical values at α • and α 1 significance levels, and ν j = n j − 1 is the degree of freedom for the jth measurement type; and n j is the number of jth type measurements. The significance levels need to be further tuned to achieve the required performance for the application under consideration. In this study, the values α • andα 1 are adopted to be 0.15 and 0.01, respectively. The classification factor is then used to scale the variance of each measurement individually as: where j refers to the observation type; i refers to the ith observation of the jth type. In this study, we have six types of measurements as described in Equation (37). Therefore, the j symbol takes values that range from one to six, and the i symbol depends on the number of the available measurements on each type. After this step, the GNSS measurements, including outliers, are down-weighted in the modified covariance matrix R. The Kalman gain and a posteriori covariance matrix are then estimated as follows: Finally, the updated state vector is estimated as: In the closed-loop scheme, Equation (44) can be simplified as: The updated state vector parameters are used to correct the IMU full mechanization outputs, resulting in the TC GNSS PPP/INS integrated solution.

GNSS Receiver and IMUs
In this study, three GNSS receivers were used. The first is the u-blox F9P GNSS receiver, while the second is the Trimble R9s geodetic receiver. The used GNSS measurements from both receivers are summarized in Section 2. As mentioned earlier, the BeiDou measurements were excluded in this study due to the low number of tracked BeiDou satellites (two or less) during the field trials. The third one is the NovAtel PwrPak7 geodetic receiver, which was used to generate the reference solutions. Two IMUs were used in the field trials, namely the Stim300 tactical-grade IMU and the xsens MTi 670 IMU industrial-grade. The specifications of both IMUs are summarized in Table 2.

Reference Solutions
The NovAtel PwrPak7 GNSS receiver and the Stim300 were used to generate the reference solutions for the three field trials. The raw GNSS measurements from PwrPak7 and the raw measurements from Stim300 were utilized along with GNSS measurements from a base station to generate a TC carrier phase-based differential GNSS (DGNSS)/INS integrated solution for each field trial. The TC DGNSS/INS solution was generated in post-processing mode using NovAtel's waypoint inertial explorer (IE) software. Three virtual reference stations (VRS) were constructed using Cansel's GNSS permanent network (CAN-NET) in the Toronto area to provide raw GNSS measurements, which were used to generate the reference solutions for the three field trials.

Processing Scheme
To thoroughly evaluate the performance of the proposed integrated system, four different TC GNSS PPP/INS integrated solutions were developed as follows: Ublox-Xsens: The u-blox GPS, GLONASS, and Galileo measurements were tightlycoupled integrated with the xsens MTi670 accelerometers and gyro measurements. This represents the main system that was developed in this study to provide an ultra-low-cost integrated solution.
Trimble-Xsens: In this solution, Trimble GNSS measurements and xsens MTi670 accelerations and angular rotation measurements were integrated in a tightly-coupled mode. This system examined the performance of the integrated system when an industrialgrade IMU is used.
Ublox-Stim: This solution took advantage of the tactical grade of the Stim300 to provide better accuracy, especially for the attitude solution.
Trimble-Stim: This system was adopted in this study to provide the most precise autonomous positioning and attitude solution. The performance of the above systems was assessed in comparison with the Trimble-Stim system.
It should be pointed out that the lever arm between each GNSS antenna and the corresponding IMU is precisely measured and compensated for during the processing process. It should be pointed that cost of the Ublox-Xsens system is at least 10 times lower than the Trimble-Stim system.

Sensors Setup
All sensors were rigidly fixed on a wooden panel and firmly mounted on the vehicle's rooftop, as shown in Figure 2. The u-blox F9P GNSS receiver was connected to a u-blox patch antenna, the Trimble GNSS receiver was connected to Trimble Zephyr 3 antenna, and the PwrPak7 was connected to an OmniStar antenna. The axes of both IMUs were aligned with the vehicle's body frame axes, in which the y-axis refers to the forward direction, the x-axis is the transversal direction, and the z-axis is completing a right-handed system.
All sensors were rigidly fixed on a wooden panel and firmly mounted on the vehicle's rooftop, as shown in Figure 2. The u-blox F9P GNSS receiver was connected to a ublox patch antenna, the Trimble GNSS receiver was connected to Trimble Zephyr 3 antenna, and the PwrPak7 was connected to an OmniStar antenna. The axes of both IMUs were aligned with the vehicle's body frame axes, in which the y-axis refers to the forward direction, the x-axis is the transversal direction, and the z-axis is completing a righthanded system.

Result Analysis
Three land vehicular field trials were conduct ed in Toronto, Canada, to assess the performance of the proposed ultra-low-cost TC triple-constellation GNSS PPP/INS integrated system. The first field trial aimed to assess the positioning solution in a suburban area under both an open-sky environment and a partial GNSS outage. Both of the second and the third field trials were carried out in combined suburban and urban areas, where the GNSS signals were partially or fully blocked.

First Field Trial
The field trial was on 22 June 2020 and lasted for about 30 min. Figure 3 shows the full trajectory of the field trial. As shown in Figure 4, the number of tracked GPS, GLONASS, and Galileo satellites ranges from 4 to 13 when the u-blox GNSS receiver is used, while the number of tracked satellites ranges from 9 to 16 when the Trimble GNSS

Result Analysis
Three land vehicular field trials were conduct ed in Toronto, Canada, to assess the performance of the proposed ultra-low-cost TC triple-constellation GNSS PPP/INS integrated system. The first field trial aimed to assess the positioning solution in a suburban area under both an open-sky environment and a partial GNSS outage. Both of the second and the third field trials were carried out in combined suburban and urban areas, where the GNSS signals were partially or fully blocked.

First Field Trial
The field trial was on 22 June 2020 and lasted for about 30 min. Figure 3 shows the full trajectory of the field trial. As shown in Figure 4, the number of tracked GPS, GLONASS, and Galileo satellites ranges from 4 to 13 when the u-blox GNSS receiver is used, while the number of tracked satellites ranges from 9 to 16 when the Trimble GNSS receiver is used. The reason for the low number of tracked satellites for the u-blox receiver is that the u-blox supports only civil signals, L1 C/A and L2C, for GPS constellation. Unfortunately, however, the L2C signals are only transmitted by modernized GPS satellites. This means that fewer GPS dual-frequency measurements were available at the time of data collection to form the conventional IF linear combination. The tracking capability of the u-blox GNSS receiver significantly affects the accuracy of the positioning solution, as shown in positioning solution due to the high variability of the number of tracked satellites when the u-blox receiver is used. The Ublox-Stim solution smoothed most of the jumps that are present during the epochs when a low number of satellites are tracked, compared to the Ublox-Xsens solution. The Trimble-based positioning solution outperforms the u-bloxbased counterpart, proving that the GNSS receiver tracking capability is the dominant factor for the accuracy of the positioning solution. Additionally, both of the Trimble-Xsens and Trimble-Stim solutions show similar positioning performance.   Table 3 shows the root-mean-square (RMS) of positioning errors of the Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions in the east, north, and up directions. It should be pointed out that the first three minutes, which represent the solution convergence, are excluded from the results. As can be seen, the quality of the IMU grade is the dominant factor that degrades the positioning accuracy for the u-bloxbased solutions. The Ublox-Xsens solution achieves sub-meter positioning accuracy in all Ublox-Xsens solution. The Trimble-based positioning solution outperforms the u-blox based counterpart, proving that the GNSS receiver tracking capability is the dominan factor for the accuracy of the positioning solution. Additionally, both of the Trimble-Xsen and Trimble-Stim solutions show similar positioning performance.   Table 3 shows the root-mean-square (RMS) of positioning errors of the Ublox-Xsens Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions in the east, north, and up directions. It should be pointed out that the first three minutes, which represent th solution convergence, are excluded from the results. As can be seen, the quality of th IMU grade is the dominant factor that degrades the positioning accuracy for the u-blox based solutions. The Ublox-Xsens solution achieves sub-meter positioning accuracy in al            Table 3 shows the root-mean-square (RMS) of positioning errors of the Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions in the east, north, and up directions. It should be pointed out that the first three minutes, which represent the solution convergence, are excluded from the results. As can be seen, the quality of the IMU grade is the dominant factor that degrades the positioning accuracy for the u-bloxbased solutions. The Ublox-Xsens solution achieves sub-meter positioning accuracy in all directions. On the other hand, the positioning accuracy is improved to reach decimeterlevel in the east, north, and up directions when the Stim300 is utilized with the u-blox F9P GNSS receiver. This represents an improvement of 15%, 59%, and 58% in the east, north, and up directions, respectively. The Trimble-based solutions achieve better positioning accuracy compared to the u-blox-based counterparts. The   Figure 8 shows the attitude errors for the pitch, roll, and azimuth angles of the xsensbased integrated solutions, while attitude errors of the Stim-based integrated solutions are presented in Figure 9. It can be seen that the attitude solution performance is significantly improved when the Stim300 is used, regardless of whether the Trimble or u-blox GNSS receiver is used.   Figure 8 shows the attitude errors for the pitch, roll, and azimuth angles of the xsensbased integrated solutions, while attitude errors of the Stim-based integrated solutions are presented in Figure 9. It can be seen that the attitude solution performance is significantly improved when the Stim300 is used, regardless of whether the Trimble or u-blox GNSS receiver is used.   The RMS and the mean of the attitude errors of the Ublox-Xsens, Trimble-Xsens Ublox-Stim, and Trimble-Stim integrated solutions are presented in Table 4. The RMS o the attitude errors of the Ublox-Xsens solution is 0.293°, 0.331°, and 3.470° for pitch, roll and azimuth, respectively, compared to 0.292°, 0.363°, and 3.475°, respectively, for the Trimble-Xsens solution. As can be seen in Table 4, the IMU grade is the dominant facto that affects the attitude solution accuracy, regardless of whether a geodetic or low-cos  Table 4. The RMS of the attitude errors of the Ublox-Xsens solution is 0.293 • , 0.331 • , and 3.470 • for pitch, roll, and azimuth, respectively, compared to 0.292 • , 0.363 • , and 3.475 • , respectively, for the Trimble-Xsens solution. As can be seen in Table 4, the IMU grade is the dominant factor that affects the attitude solution accuracy, regardless of whether a geodetic or low-cost GNSS receiver is used. For example, the RMS of the attitude errors of the Ublox-Stim is improved by 83%, 86%, and 90% for pitch, roll, and azimuth angles, respectively, in comparison with the Ublox-Xsens solution. The same improvement in the attitude solution accuracy is achieved when the Stim300 is utilized with the Trimble GNSS receiver instead of the xsens MTi670.

Second Field Trial
The second field trial was conducted in an urban area with a challenging environment. The trial lasted for about 45 min, and the full trajectory is shown in Figure 10. The number of the tracked satellites ranged from 1 to 17 when the u-blox GNSS receiver was used, while the number of the tracked satellites ranged from 7 to 19 when the Trimble GNSS receiver was used, as shown in Figure 11.
Geomatics 2021, 1, FOR PEER REVIEW 18 accuracy is achieved when the Stim300 is utilized with the Trimble GNSS receiver instead of the xsens MTi670.

Second Field Trial
The second field trial was conducted in an urban area with a challenging environment. The trial lasted for about 45 min, and the full trajectory is shown in Figure 10. The number of the tracked satellites ranged from 1 to 17 when the u-blox GNSS receiver was used, while the number of the tracked satellites ranged from 7 to 19 when the Trimble GNSS receiver was used, as shown in Figure 11.   eomatics 2021, 1, FOR PEER REVIEW 1 accuracy is achieved when the Stim300 is utilized with the Trimble GNSS receiver instead of the xsens MTi670.

Second Field Trial
The second field trial was conducted in an urban area with a challenging environ ment. The trial lasted for about 45 min, and the full trajectory is shown in Figure 10. The number of the tracked satellites ranged from 1 to 17 when the u-blox GNSS receiver was used, while the number of the tracked satellites ranged from 7 to 19 when the Trimble GNSS receiver was used, as shown in Figure 11.                Table 5 shows the root-mean-square (RMS) of positioning errors of the Ublox-Xsens Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions in the east, north, and up directions. It should be pointed out that the first three minutes, which represent th solution convergence, are excluded from the RMS values. The Ublox-Xsens solution     Figure 15 shows the attitude errors for the pitch, roll, and azimuth angles of the Xsens-based integrated solutions, while the attitude errors of the Stim-based solutions are presented in Figure 16.
For the Trimble-Xsens and Ublox-Xsens solutions, the error in the pitch and roll angles is typically less than 2 • , while it reaches 6 • for the azimuth. It can be seen that the attitude solution is significantly improved when the Stim300 is used, regardless of whether the Trimble or u-blox GNSS receiver is used. For both of the Trimble-Stim and Ublox-Stim solutions, the error in the pitch and roll angles is typically less than 0.1 • , while it reaches 0.5 • for the azimuth.
The RMS and the mean of the attitude errors of Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions are presented in Table 6. The RMS values of the attitude errors of the Ublox-Xsens solution are 1.623 • , 1.732 • , and 2.254 • for the pitch, roll, and azimuth, respectively, compared to 1.626 • , 1.730 • , and 2.270 • , respectively, for the Trimble-Xsens solution. It can be seen that the IMU grade is the dominant factor that affects the accuracy of the attitude solution, regardless of whether a geodetic or a low-cost GNSS receiver is used, as concluded from the results of the first trial. For example, the RMS values of the attitude errors of the Ublox-Xsens solution are improved from 1.623 • , 1.732 • , and 2.254 • for pitch, roll, and azimuth angles, respectively, to 0.043 • , 0.051 • , and 0.373 • for the Ublox-Stim solution. The attitude accuracy is improved by 97%, 96%, and 87% for the pitch, roll, and azimuth angles, respectively, when the Stim300 is used with the Trimble GNSS receiver instead of the xsens MTi670.  Figure 15 shows the attitude errors for the pitch, roll, and azimuth angles of the Xsens-based integrated solutions, while the attitude errors of the Stim-based solutions are presented in Figure 16.     The RMS and the mean of the attitude errors of Ublox-Xsens, Trimble-Xsens, Ublo Stim, and Trimble-Stim integrated solutions are presented in Table 6. The RMS values the attitude errors of the Ublox-Xsens solution are 1.623°, 1.732°, and 2.254° for the pitc roll, and azimuth, respectively, compared to 1.626°, 1.730°, and 2.270°, respectively, f the Trimble-Xsens solution. It can be seen that the IMU grade is the dominant factor th affects the accuracy of the attitude solution, regardless of whether a geodetic or a low-co

Third Field Trial
The third field trial was conducted in a combined urban and suburban area with a challenging environment that included tall buildings and urban foliage (about 3 m). The trial lasted for about 34 min and the full trajectory of the field trial is shown in Figure 17. The number of tracked satellites ranged from 0 to 19 when the u-blox GNSS receiver was used, while the number of tracked satellites ranged from 4 to 19 when the Trimble GNSS receiver was used, as shown in Figure 18.

Third Field Trial
The third field trial was conducted in a combined urban and suburban area with a challenging environment that included tall buildings and urban foliage (about 3 m). The trial lasted for about 34 min and the full trajectory of the field trial is shown in Figure 17. The number of tracked satellites ranged from 0 to 19 when the u-blox GNSS receiver was used, while the number of tracked satellites ranged from 4 to 19 when the Trimble GNSS receiver was used, as shown in Figure 18.  As shown in Figures 19-21, the Ublox-Stim solution shows better performance than the Ublox-Xsens solution, especially in the east direction. Additionally, the Ublox-Stim

Third Field Trial
The third field trial was conducted in a combined urban and suburban area with a challenging environment that included tall buildings and urban foliage (about 3 m). The trial lasted for about 34 min and the full trajectory of the field trial is shown in Figure 17. The number of tracked satellites ranged from 0 to 19 when the u-blox GNSS receiver was used, while the number of tracked satellites ranged from 4 to 19 when the Trimble GNSS receiver was used, as shown in Figure 18.  As shown in Figures 19-21, the Ublox-Stim solution shows better performance than the Ublox-Xsens solution, especially in the east direction. Additionally, the Ublox-Stim As shown in Figures 19-21, the Ublox-Stim solution shows better performance than the Ublox-Xsens solution, especially in the east direction. Additionally, the Ublox-Stim solution smoothed the majority of the jumps that are present during the epochs when a low number of satellites are tracked, while some jumps reach more than a meter for the Ublox-Xsens solution. It is worth mentioning that at some points in time, the positioning errors in the east and north directions of both of the Ublox-Xsens and Ublox-Stim solutions exceed one meter, which is mainly due to the high variability of the tracked satellites. The Trimble-based positioning solutions outperform the u-blox-based solutions, which is mainly due to the tracking capability of the GNSS receiver.
Geomatics 2021, 1, FOR PEER REVIEW 23 solution smoothed the majority of the jumps that are present during the epochs when a low number of satellites are tracked, while some jumps reach more than a meter for the Ublox-Xsens solution. It is worth mentioning that at some points in time, the positioning errors in the east and north directions of both of the Ublox-Xsens and Ublox-Stim solutions exceed one meter, which is mainly due to the high variability of the tracked satellites. The Trimble-based positioning solutions outperform the u-blox-based solutions, which is mainly due to the tracking capability of the GNSS receiver.          The root-mean-square (RMS) and mean of positioning errors of the Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions in the east, north, and up directions are presented in Table 7. The RMS values of the Ublox-Xsens positioning solution are 1.029 m, 0.687 m, and 0.466 m in the east, north, and up directions, respectively. Due to the high variability of the tracked satellites when the u-blox GNSS receiver is used, the Ublox-Stim solution shows better positioning performance than the Ublox-Xsens counterpart. The Ublox-Xsens solution is improved by 27%, 18%, and 21% in the east, north, and up directions, respectively, when the Stim300 IMU is used with the u-blox F9P GNSS receiver. As concluded from the first two field trials, the Trimble-based solutions achieve better positioning accuracy than the u-blox-based counterparts. The   Figure 22 shows the attitude errors for the Xsens-based integrated solutions, while the attitude errors for the Stim-based solutions are presented in Figure 23.
The attitude errors of both of the Trimble-Xsens and Ublox-Xsens solutions are less than 2 • for the pitch and roll components and less than 6 • for the azimuth. It can be seen that for both of the Trimble-Stim and Ublox-Stim solutions, the attitude errors are less than 0.1 • for the pitch and roll components and less than 0.5 • for the azimuth. However, at some points in time, the azimuth angle errors for the Ublox-Stim solution exceed 0.5 • , which is mainly due to the degradation in the positioning accuracy. m for the Trimble-Xsens solution. Additionally, the RMS errors of the Ublox-Stim posi tioning solution are improved from 0.749 m, 0.562 m, and 0.565 m in the east, north, and up directions, respectively, to 0.212 m, 0.289 m, and 0.132 m for the Trimble-Stim.  Figure 22 shows the attitude errors for the Xsens-based integrated solutions, whil the attitude errors for the Stim-based solutions are presented in Figure 23.    The RMS and the mean of the attitude errors of Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions are presented in Table 8. The RMS values o the attitude errors of the Ublox-Xsens solution are 0.719°, 0.334°, and 2.980° for the pitch roll, and azimuth, respectively, compared to 0.718°, 0.335°, and 2.985°, respectively, fo the Trimble-Xsens solution. As concluded from the first two field trials, the IMU grade i The RMS and the mean of the attitude errors of Ublox-Xsens, Trimble-Xsens, Ublox-Stim, and Trimble-Stim integrated solutions are presented in Table 8. The RMS values of the attitude errors of the Ublox-Xsens solution are 0.719 • , 0.334 • , and 2.980 • for the pitch, roll, and azimuth, respectively, compared to 0.718 • , 0.335 • , and 2.985 • , respectively, for the Trimble-Xsens solution. As concluded from the first two field trials, the IMU grade is the dominant factor that affects the attitude solution accuracy, regardless of whether a geodetic or low-cost GNSS receiver is used. For example, the RMS values of the attitude errors of the Ublox-Xsens are improved from 0.719 • , 0.334 • , and 2.980 • for the pitch, roll, and azimuth angles, respectively, to 0.034 • , 0.047 • , and 0.466 • for the Ublox-Stim solution. Additionally, the attitude accuracy for the Trimble-Xsens solution is improved by 97%, 95%, and 96% for the pitch, roll, and azimuth angles, respectively, compared to the Trimble-Stim solution.

Conclusions
In this study, an ultra-low-cost TC triple-constellation GNSS PPP/INS integrated system was developed. The developed system tightly integrates data from the newly developed low-cost dual-frequency u-blox F9P GNSS receiver and the industrial-grade xsens MTi670 IMU. The performance of the developed system was assessed through three land vehicular field trials. The first field trial was carried out to assess the integrated solution in an open-sky environment, while the second and third field trials were carried out to assess the integrated solution under challenging environments. It was shown that the proposed Ublox-Xsens system achieved sub-meter-level PPP accuracy in an open-sky environment, while it achieved meter-level PPP accuracy in challenging environments. The relatively low PPP accuracy is due mainly to the tracking capability of the u-blox F9P GNSS receiver, which led to tracking a limited number of GPS satellites at the time of observations. The average RMS of the attitude errors of the integrated solution was 0.878 • , 0.799 • , and 2.901 • for the pitch, roll, and azimuth angles, respectively. The positioning accuracy of the integrated solution showed a significant improvement when the Trimble R9s geodetic receiver was used along with the xsens MTi670 IMU. The positioning solution of the Trimble-Xsens integrated system consistently achieved decimeter-level positioning accuracy under both open-sky and challenging environments. The Ublox-Stim integrated system, in which the u-blox F9P GNSS receiver and Stim300 tactical-grade IMU were used, achieved comparable positioning accuracy to the Ublox-Xsens solution. In contrast to the Ublox-Xsens solution, however, the Ublox-Stim integration smoothed the majority of the jumps that were present in the positioning solution during the epochs when a low number of satellites were tracked. Moreover, the accuracy of the obtained attitude solution was significantly improved when the Stim300 IMU was used. The average RMS of the attitude errors of the Ublox-Stim solution was 0.042 • , 0.048 • , and 0.394 • for the pitch, roll, and azimuth angles, respectively. The Trimble-Stim integrated system showed the best positioning and attitude solution accuracy. However, the relatively high cost of the Trimble-Stim integrated system might limit its usage for a wide range of applications.