Performance Analysis of MADOCA-Enhanced Tightly Coupled PPP/IMU

3.1 Doppler-Aided Dual-Frequency Observation Model

This study employed multi-constellation MADOCA-PPP to maximize the utilization of the augmentation message. To obtain precise estimates of the receiver’s position, both the GNSS pseudorange (P) and carrier phase (Φ) were used. These two characteristics can be combined to enhance the performance of the system. Generally, the undifferenced and uncombined PPP observation model can be easily scaled to integrate more frequencies (Vana & Bisnath, 2023). However, ionospheric constraints are required to eliminate negative impacts on the navigation system, which can cause positioning offsets on the order of meters (Kubo et al., 2024; Mayer et al., 2008). The observation model of MADOCA-PPP stated in the QZSS interface control document (ICD) (Cabinet Office, 2024) is followed. Therefore, in this study, the L1 and L2 bands for GPS, GLONASS, and QZSS were adopted for ionospheric-free observation measurements. Furthermore, Doppler measurements enable precise velocity estimation with less noise. A Doppler-aided model can smoothen position estimations during the initialization stage, thus reducing the positioning error (Chi et al., 2023). Hence, the multi-constellation ionospheric-free MADOCA-PPP can be expressed as follows:

PIF=‖rM−ra‖+c(dtr−dtMq)+Tslant+ϵPIFΦIF=‖rM−ra‖+c(dtr−dtMq)+Tslant+NIFλIF+ϵΦIF−λ1D1=‖r˙M−r˙a‖+c(dt˙r−dt˙Mq)+T˙slant+I˙+ϵD7

with:

‖r˙M−r˙a‖=(ers)T[vMs(t)−va(t)]+ωec(vysxa+ysvx,r−vxsya−xsvy,r)8

where the abbreviation IF denotes the ionospheric-free combination, the superscript q denotes the specific constellations, and G, R, and J indicate GPS, GLONASS, and QZSS, respectively. r_M and dtMq are the precise satellite position and clock bias, respectively, and the subscript M indicates that the parameters are corrected by MADOCA. Note that r_M = [x^s, y^s, z^s]^T in Equation (8). r_a and dt_r are the antenna position and receiver clock bias, respectively, with r_a = [x_a, y_a, z_a]^T. T_slant is the troposphere total delay, and I is the ionospheric delay or advance. N_IF and λ_IF are the float ambiguity and wavelength with the ionospheric-free combination. ϵ for each measurement is the sum of the noise and unexpected effect, for example, NLOS and multipath.

A superscript dot indicates the rate of change in Equation (7) and Equation (8). Doppler measurements in only the L1 band (denoted D₁ with wavelength λ₁) are used. The rates of change of the atmospheric terms, Ṫ_slant and İ, are negligible. ers is the line-of-sight vector from the receiver to the satellite, and ω_e is the Earth’s rotation rate in scalar form. vMs(t) is the satellite velocity vector with [vxs,vys,vzs]T, and v_a(t) is the receiver velocity vector with [v_a,x, v_a,y, v_a,z]^T, both in the ECEF frame.

In addition, the slant troposphere delay comprises the zenith hydrostatic (dry) delay (T_d) and zenith wet delay (T_w) of air and can be calculated via a corresponding mapping function (MacMillan & Ma, 1997). The magnitude of the dry component may be several times larger than that of the wet component (Gunning, 2021). However, the dry delay can generally be precisely predicted via existing models, whereas the wet delay might have a larger uncertainty. Consequently, the delta wet component (δT^w) can be modeled as a random walk (Gao et al., 2017) or the entire wet delay (T^w) can be directly estimated by the EKF. The total slant delay can be rewritten as follows using the zenith component:

Tslant=md(E)Td+mw(E)(Tw+δT^w)=md(E)Td+mw(E)T^w9

where m_d(E) and m_w(E) are the satellite elevation-angle-dependent mapping function for the dry and wet delay components, respectively.

Based on the GPS clock bias, intersystem bias (ISB) might be caused by differences between the time systems in each constellation (Dalla Torre & Caporali, 2014). This parameter should also be considered when processing multi-GNSS observations. The QZSS time system is intended to be synchronized with that of GPS. Hence, only the ISB between GLONASS and GPS must be calculated, as follows:

δISBG,R=dtrR−dtrG10

where δISB_G,R is the time difference between GLONASS and GPS and the superscripts R and G represent GLONASS and GPS, respectively.

3.2 MADOCA-PPP EKF Formulation

In this study, a kinematic, time-uncorrelated EKF was selected. According to the observation model described in Section 3.1, the estimated system state vector (x_GNSS) for standalone MADOCA-PPP can be expressed as follows:

xGNSS=[r3×1v3×1dtrGδISBG,Rδt˙rT^wNIFλIF](10+n)×1T11

Here, r and v are the receiver’s position and velocity, respectively; both are 3 × 1 vectors in the ECEF frame. δṫ_r represents the receiver clock drift, which can be estimated by Doppler observations using multi-GNSS data (Gao et al., 2016). The state vector has n +10 dimensions given n observable satellites, as stated in Equation (11), where the subscript (10 + n) × 1 indicates the total number of states.

The measurement noise for pseudorange, carrier-phase, and Doppler measurements are elevation-angle-dependent (el) and are assigned as follows (Takasu & Yasuda, 2009):

σP2=σtrop2+σeph2+αp(a2+b2sin2(el))σΦ2=σtrop2+σeph2+αΦ(a2+b2sin2(el))σD2=σeph2+αD(a2+b2sin2(el))12

Here, small values can be assigned to the standard deviation of the tropospheric delay σtrop2. σeph2 is the SSR user range accuracy of ephemeris provided by the augmentation product. α is the error ratio between each measurement, which can be assigned as 300 for α_P, 3 for α_D, and 1 for α_Φ. a and b are the standard deviation of the carrier phase and its elevation-angle-dependent term, respectively. Fine parameter tuning should be executed if the measurement quality is not consistent between constellations; this issue will be addressed in Section 5.3.2.

A conventional EKF was applied to estimate the system state vector in Equation (11). For brevity, the propagation model in the EKF will be described in Section 4.

4.1 Inertial Propagation

INSs provide redundant information, such as precise attitude and velocity estimation. However, if a standalone INS is used, the error in the IMU drift over time severely degrades the navigation performance. Therefore, the error model and its mechanization should be clearly specified before integration. The major error sources are the bias and scale factor, which both exist in accelerometers and gyroscopes. The measurement model (Groves, 2013) can be represented as follows:

[f^bω^b]=[(I+δSa)fb+ba+wa(I+δSg)ωb+bg+wg]13

where f^b and ω^b are the raw measurements from the accelerometer and gyroscope in the body frame (denoted as the b-frame), respectively, f^b and ω^b are the unbiased acceleration and angular rate, I is a unit matrix, δS is the scale factor error, b is the bias, w is the white noise in the sensors, and the subscripts a and g represent the accelerometer and gyroscope, respectively.

The clean measurement without bias can be obtained from Equation (13), which can derive the propagation of attitude, velocity, and position. In this study, the IMU dynamic was propagated in the ECEF frame and derived from small perturbation analysis (Groves, 2013; Shin, 2005). Therefore, the error model can be described as follows:

[δψ˙eδv˙eδr˙e]=[Cbeδωb−ωiee×δψeCbeδfb−(Cbefb)×δψe−2ωiee×δve+δgeδve]14

where the superscript e denotes the ECEF frame and δψ˙e,δv˙e, and δṙ^e are the time derivatives of the attitude, velocity, and position errors, respectively. Cbe is the transformation matrix from the b-frame to the ECEF frame, and δω^b and δf^b are the systematic errors of the gyroscope and accelerometer in the b-frame, which can be derived from Equation (13). ωiee is the Earth rotation rate vector, and δg^e is the gravity error vector.

This study implements the error-state EKF (ES-EKF) to propagate the state vector. The major difference between the EKF and ES-EKF is that the ES-EKF only calculates the error terms in each state relative to the previous epoch. The time update is always conducted after mechanization to update the covariance and the IMU-corresponding state vector. In this study, the IMU state vector (δx_IMU) was formulated as follows:

δxIMU=[δψeδveδrebabgSaSg]21×1T15

The time update in discrete-time form (Groves, 2013) can be written as follows:

δxIMU,k=Φk/k−1δxIMU,k−1+wk−116

Pk−=Φk/k−1Pk−1+Φk/k−1T+Qk−117

where Φ_k/k–1 is the state transition matrix from epoch k −1 to epoch k, w_k−1 is the system noise matrix, P_k is the error-state covariance matrix, and Q_k−1 is the noise covariance matrix.

In addition, the sensor bias and scale factor are both varied gradually; thus, the variation in the discrete form can be modeled as a first-order Gauss–Markov process:

[ba,kbg,kSa,kSg,k]=e−Δt/Tcorr[ba,k−1bg,k−1Sa,k−1Sg,k−1]+[wbwbwsws]18

where Δt is the time interval, T_corr is the correlation time, and w_b and w_s are the white noise power spectral densities corresponding to the bias and scale factor, respectively. Therefore, the state transition matrix Φ_k/k−1 is derived from Equation (14) and Equation (18).

4.2 Tightly Coupled PPP/IMU

The measurement update is implemented after the GNSS and IMU measurements become available. The tightly coupled scheme simultaneously estimates the PPP and IMU state during the integration process. According to the above-described model, the combined state vector can be constructed as follows:

δxTC=[δxIMUδtrGδISBG,Rδt˙rT^wNIFλIF](25+n)×1T19

where δx_TC is the error state in the tightly coupled scheme and δx_IMU is as defined in Equation (15). The vector has (25 + n) × 1 dimensions given n observable satellites. Additionally, the Doppler measurement is integrated to calibrate the velocity of the vehicle. Therefore, the measurement model in the EKF can be expressed as follows:

δZk=HTCδxTC+wk20

where δZ_k is the residual of the measurement vector at epoch k, H_TC is the geometry matrix in the tightly coupled scheme, and w_k is the measurement noise vector, which is obtained from the measurement covariance and assumed to have a zero mean. The lever-arm effect should be carefully removed while considering the integration model; then, the measurement vector can be described as follows:

δZk=[PGNSS−P^IMUΦGNSS−Φ^IMUDGNSS−D^IMU]3n×1=[PGNSSΦGNSSDGNSS]−[PIMU−|δpri|ΦIMU−|δpri|DIMU−|δvri|]21

with:

|δpri|=Cbelb|δvri|=Cbe(ωb)×lb−(ωiee)×Cbelb22

where P_GNSS, Φ_GNSS, and D_IMU are as shown in Equation (7) with MADOCA augmentation. P^IMU,Φ^IMU, and D^IMU are the predicted pseudorange, carrierphase, and Doppler measurements using inertial propagation as introduced in Section 4.1. In Equation (22), |δpri| and |δvri| are the lever-arm compensation vector corresponding to range and velocity, respectively. l^b is the lever arm in the b-frame, from the IMU to the GNSS antenna. The superscript “×” denotes a skew-symmetric matrix.

According to Equation (7) and Equation (21), the geometry matrix can be derived from the partial derivatives of each element in the observation model, which can be defined as follows:

HTC=[HPHΦHD]3n×(25+n)=[−(ers)T(Cbelb)×0(ers)T01×12Hδtrmw0−(ers)T(Cbelb)×0(ers)T01×12Hδtrmw1HδψHδvHδrHimuHδt˙r00]3n×(25+n)23

with:

Hδψ=−[(ers)T−ωec(rM)TA][Cbe(ωb)×lb)×−(ωiee)×(Cbelb)×]Hδv=(ers)T−ωec(rM)TAHδr=−ωec(rM)TAHδbg=−((ers)T−ωec(rM)TA)(Cbe(lb)×)Himu=[01×3Hδbg01×301×3]1×12A=[0−10100000]24

where r_M is the satellite position vector in the ECEF frame, with the same notation as in Equation (8). Hδtr corresponds to the receiver clock bias term, which can be [1, 0, 0] or [0, 1, 0]. Hδt˙r corresponds to the clock bias drift term, which is [0, 0, 1].

4.3 Motion Constraint Models

The error accumulation in an IMU increases rapidly with time. The error drift can be reduced by using motion constraint models based on prior knowledge from the detection of certain motions. In this section, three motion constraint models are introduced with regard to lateral, static, and angular rate motion. In land applications, a ground vehicle is assumed to be moving forward, and no sideslip is assumed to occur with respect to the b-frame. The non-holonomic constraint (NHC) model was constructed based on these assumptions. In addition, the ZUPT and ZARU models are common strategies for a static vehicle. These strategies are necessary to bound the error drift, particularly for vehicle or pedestrian dead-reckoning applications. In this paper, the generalized likelihood ratio test detector was employed to determine whether a vehicle is moving. The fundamental theorem for this test can be found in the study by Skog et al. (2010).

4.4 Proposed System Architecture

Figure 4 shows the proposed framework for the MADOCA-enhanced tightly coupled PPP/IMU. First, the standalone MADOCA-PPP was applied for 5 s to perform absolute positioning and initial IMU alignment. In addition, several parameters, such as the heading angle and velocity of the vehicle, were initialized using the GNSS. The IMU was activated after system initialization. The raw accelerometer and gyroscope measurements were compensated for using the error model described in Section 4.1. Subsequently, the IMU states in Equation (15) were predicted in the ECEF frame at a frequency higher than the GNSS observation frequency. After the inertial propagation, the ES-EKF time update is always performed to compute the corresponding covariance. Thereafter, the bias and scale factor are provided as feedback to correct the IMU errors in the next epoch. The GNSS receiver obtains raw measurements of the dual-frequency pseudorange, carrier phase, and single-frequency Doppler shift from the multi-constellation. A cycle slip should be detected before ionospheric-free correction is performed; this detection is commonly performed using geometry-free and Melbourne–Wubbena combination (Teunissen & Montenbruck, 2017). Subsequently, the ionospheric-free measurement vectors were constructed to eliminate ionospheric error. The orbit and clock bias correction in the MADOCA message were applied in the PPP model for precise estimation of the range from the receiver to the satellite with the broadcast ephemeris. Thereafter, the residual vectors described in Section 4.2 were computed to implement the EKF measurement update when the GNSS observations are available. In this stage, the motion constraint model was used to improve the navigation performance. Consequently, through the aforementioned steps, a MADOCA-enhanced tightly integrated position, velocity, and attitude estimation system can be obtained. Section 5 describes the experimental validation of this system and defines the related parameters for the GNSS and IMU.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 4

Tightly coupled MADOCA-PPP/IMU integration framework

5.1 Experimental Setup

In this study, experiments were performed in two stages: static and dynamic field tests. The first experiment was aimed at verifying the ability of the augmentation message combined with the broadcast ephemeris to reduce convergence time and increase navigation performance. The static experiment was conducted in an open-sky environment with the GNSS receiver NovAtel PwrPak7 and a GNSS-850 antenna, which are capable of processing multi-frequency and multi-constellation observations.

The dynamic field experiment was performed by installing the navigation system on a sport utility vehicle (SUV, as shown in Figure 5(a)) to collect raw measurements for experimental analysis and to generate a ground truth trajectory using NovAtel’s synchronous position, attitude, and navigation (SPAN) technology in dual-antenna mode. In addition, the NovAtel PwrPak7 receiver and a tactical-grade IMU, namely the LN200C, were integrated into the system platform (Figure 5(b)), and the experiment was performed at an IMU sampling rate of 100 Hz.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 5

(a) Experimental setup on an SUV; (b) experimental platform

The processing strategies for standalone MADOCA-PPP and the tightly coupled system are presented in Table 2. To evaluate the performance of MADOCA-PPP, data provided by the augmentation message for all three constellations were used. Note that our analysis did not account for Galileo satellite corrections, because the analysis was conducted before August 4, 2023, when Galileo corrections were included in the RTCM-formatted MADOCA product. Therefore, only three constellations were considered.

First, the satellite orbit and clock bias were computed from the broadcast ephemeris and then corrected based on the MADOCA product. The elevation mask angle was set at 10° to mitigate signal degradation, the sampling rate was set to 2 Hz, and all data could be acquired in real time if required. Furthermore, the instability and random walk parameters for the IMU in Table 2 were estimated using the Allan variance (El-Sheimy et al., 2008) for 6-h raw accelerometer and gyroscope measurements. Additionally, the parameters for motion constraint models were negligibly low because the vehicle remained stationary in the static motion and moved straight in the lateral motion constraint model. However, these parameters must be nonzero, as the IMU measurement noise should be accounted for in the model. These parameters were also determined in part through empirical observation or manual tuning. The dynamic field experiments were conducted in Tainan, Taiwan, and lasted for at least 30 min each. The trajectories included various obstacles, such as trees, footbridges, utility poles, and tall buildings. For each dynamic scenario, the performance of both kinematic MADOCA-PPP and its tight integration with an IMU was evaluated. The second experiment was conducted in a suburban area, and the third experiment was conducted in a GNSS-challenging urban area.

5.2 Static Scenario

For a preliminary performance assessment, the static MADOCA-PPP was validated in an open-sky environment for 6 h. The validation process lasted 43,200 epochs. The process noise of position terms in the state vector was defined to be zero because the position was previously known to be stationary. In accordance with the criterion proposed in the QZSS ICD for a dual-frequency receiver in an open-sky environment, the convergence time was defined as the time when the positioning errors in the horizontal and vertical directions reached 30 and 50 cm, respectively. Figure 6 presents a comparison of the positioning error with and without the augmentation message correction. Specifically, the upper plot shows the horizontal positioning error, and the lower plot shows the vertical direction error. This figure reveals the enhanced positioning result obtained after augmentation message correction compared with that achieved through broadcast ephemeris PPP without MADOCA. For clarity, only two combinations are shown in Figure 6. Furthermore, the inner plots present the performance of MADOCA-PPP in a 2-h window, demonstrating that MADOCA-PPP can achieve a much shorter convergence time. This result reveals a considerable difference in the precision of the satellite orbit and clock bias estimation of the adopted methods. Obtaining an accurate orbit and clock bias can remarkably improve the positioning accuracy and convergence time. Table 3 lists the RMSE and convergence time obtained by for a single-constellation (GPS, denoted by G), dual-constellation (GPS + GLONASS, denoted by G + R), and multi-constellation (GPS + GLONASS + QZSS, denotes by G + R + J) configuration. The preliminary results reveal that the augmentation message enhanced the convergence speed and positioning accuracy. Particularly, only approximately 5 min was required to reduce the error to the specified level, which is 95% lower than that required by the broadcast ephemeris-only method. The most significant enhancement in the positioning accuracy was observed in the vertical direction, as the vertical error decreased from the meter range to approximately 0.3 m. For the horizontal RMSE, which is the essential metric for terrestrial applications, the MADOCA-PPP can achieve decimeter-level accuracy, and the performance was comparable to that of conventional PPP. However, in the MADOCA-PPP configuration, the single-constellation test slightly outperformed the dual-constellation test, which might be attributed to the inconsistent measurement quality of the constellations. This aspect is discussed further in Section 5.3.2.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 6

Position error for the static PPP with and without augmentation

The inner plots present the MADOCA-PPP results for a 2-h window, demonstrating a shorter convergence time with augmentation.

Brdc: Using broadcast ephemeris without MADOCA

The convergence time and 95th percentile solution for each configuration are presented in Table 3. The MADOCA-PPP approach reduced the 95th percentile of the horizontal and vertical solution from meters to 0.2 m, and the performance was similar to that of the MADOCA configurations. The most significant enhancement was observed for the vertical direction. Compared with the several hours required to converge the navigation filter toward decimeter-level accuracy without augmentation, the MADOCA-PPP approach could converge the positioning error within 12 min, achieving the shortest convergence time of 267 s for the multi-constellation configuration.

5.3 Suburban Scenario

5.3.1 Standalone MADOCA-PPP

The first dynamic field experiment for a vehicle was conducted in a suburban environment for a distance of approximately 18.5 km; the test lasted 30 min, and the trajectory of the vehicle is depicted in Figure 7(a). The red points denote the start and end points of the trajectory, and the arrows indicate the vehicle orientation. Some images recorded from a dashcam are shown in the insets of Figure 7(a); these images indicate the effects of environmental factors on the positioning performance. This field experiment was conducted in a suburban environment, and low-lying environmental structures, such as buildings, trees, and poles, may have occluded the signal around the receiver, leading to slightly polluted measurements.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 7

Details of the first field experiment: (a) vehicle trajectory in the suburban scenario; (b) number of visible satellites in view and PDOP in the suburban scenario

Figure 7(b) presents the number of satellites in view and the position dilution of precision (PDOP) for different constellation combinations. The average number of available GPS satellites and the average PDOP were approximately 8.64 and 2.45, respectively. The visibility for the dual constellation with GLONASS was higher than that for the single constellation: the average number of visible satellites in the dual constellation was 13.71, and the PDOP for this constellation decreased to 1.65. In the multi-constellation test with QZSS, the average number of visible satellites increased by only 3 (to 16.62) compared with that for the dual constellation; this small increase was attributed to the limited number of QZSS satellites. Moreover, the average PDOP in the multi-constellation test was 1.42, which is marginally lower than that in the dual-constellation test.

Figures 8(a) and 8(b) show the horizontal and vertical positioning error and cumulative distribution function (CDF) of the solutions for each configuration; the values of the 95th percentile of solutions are denoted in both figures. All configurations achieved submeter-level and lane-level accuracy, with the 67th and 95th percentiles of cross-direction errors being less than 1 and 1.5 m, respectively. Moreover, the horizontal and vertical RMSE values, presented as (horizontal RMSE, vertical RMSE), were (0.47, 0.88 m) for GPS-only, (0.38, 1.31 m) for GPS + GLONASS, and (0.33, 1.31 m) for GPS + GLONASS + QZSS. The multi-constellation configuration achieved the smallest horizontal RMSE and 95th percentile of solution. However, its increased observation does not imply that it exhibited the best overall performance.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 8

Standalone MADOCA-PPP in a suburban area with different configurations: (a) horizontal error and solution distribution; (b) vertical error and solution distribution

5.3.2 Residual Analysis

As indicated by the results for the standalone MADOCA-PPP, the multi- and dual-constellation configurations do not perform better than the GPS-only configuration, especially in the vertical direction. The positioning error is positively correlated with the dilution of precision (DOP). In general, a reduced DOP was observed with increasing observations owing to the values in the geometry matrix (Figure 7(b)). However, the result contradicts the inference. This may be attributed to the measurement noise for GLONASS, which is greater than that for other constellations. Additionally, the measurement quality might be affected by satellite orbit and clock errors, even after MADOCA correction. To analyze the measurement quality for each constellation, measurement residuals are commonly used (Richardson et al., 2016). The ionospheric-free measurement residuals of the code and carrier phases for GPS, GLONASS, and QZSS calculated using Equation (7) are shown in Figure 9. The figure indicates that GLONASS exhibited the largest measurement residuals, with an RMSE of 12.09 m for the code phase and 0.06 m for the carrier phase. Additionally, owing to its high elevation angle, QZSS exhibited stable measurement quality, with an RMSE of 2.78 and 0.01 m for the code and carrier phases, respectively. Lastly, the GPS RMS residual was 4.39 m for the code phase and 0.02 m for the carrier phase, which is slightly higher than those of the QZSS satellites. Some larger residual values for the code phase occur because of rising satellites. The red squares in Figure 9(d) indicate the residual RMSE; the error bars denote the standard deviation for each constellation. This figure further supports the aforementioned conclusions. The standard deviations of the code phase for GPS, GLONASS, and QZSS were 1.29, 1.63, and 1.06 m, respectively, whereas those of the carrier phase for GPS, GLONASS, and QZSS were 0.006, 0.024, and 0.007 m, respectively.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 9

Measurement residuals and statistics: (a) GPS; (b) GLONASS; (c) QZSS; (d) error bars for the three target constellations

STD: standard deviation

The residual analysis indicates that the measurement quality for GLONASS was lower than that of GPS or QZSS. The measurement noise of GLONASS should be specified as higher, whereas it exhibited larger RMSE and standard deviation values for the measurement residual. Therefore, it can be fine-tuned in the EKF as follows:

σP,GLO2=KGLOσP2σΦ,GLO2=KGLOσΦ229

where K_GLO is a constant that amplifies the measurement noise of the pseudorange, σP,GLO2, and carrier phase, σΦ,GLO2, in GLONASS. In this study, K_GLO was set to 5. This value was selected as a middle ground on the basis of empirical and manual tuning. The benefit of adding the GLONASS constellation was limited if a larger value of K_GLO was assigned. In addition, QZSS had small residual biases. The measurement noise of QZSS was slightly fine-tuned with a scaling factor of 1.5 according to the same method for amplifying the measurement noise of GLONASS. This value of 1.5 is not universal, and users can adopt a different value according to their use case.

The positioning errors obtained after the aforementioned modification are shown in Figure 10. Although an increase in the GLONASS measurement noise reduced the positioning errors, the convergence speed decreased because of an increase in the value of the 95th percentile of solution. Thus, the navigation filter tended to distrust GLONASS measurements. However, this modification improved the RMSE in the horizontal plane and vertical direction, to (0.50, 0.80 m) for GPS + GLONASS and (0.35, 0.68 m) for GPS + GLONASS + QZSS. This improvement in the RMSE with an increase in the number of constellations is more consistent and demonstrates that the parameter selection is reasonable. Table 4 presents the overall statistics for the standalone PPP, along with the results for the MADOCA-unaided broadcast ephemeris, denoted as Broadcast-PPP. The MADOCA product was the key factor for reducing the positioning error from several meters to less than 1 m, achieving lane-level accuracy.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 10

Standalone MADOCA-PPP after parameter tuning: (a) horizontal error and solution distribution; (b) vertical error and solution distribution

View this table:

• View inline
• View popup

TABLE 4

Comparison of the Position RMSE and Solution Distribution for Multiple Configurations

5.3.3 MADOCA-Enhanced PPP/IMU

The navigation performance can be further improved by tightly integrating the IMU. Figure 11 shows the positioning error and CDF of the solution in the horizontal plane and vertical direction, respectively. The results reveal that the sensor integration notably improved the navigation accuracy. The multi-constellation configuration achieved the best performance in the horizontal plane, with a 95th percentile of solution distribution of 0.56 m, outperforming the standalone MADOCA-PPP; however, in the tight integration framework, all configurations exhibited similar performance if a sufficient number of satellites was observed. The RMSE values in the horizontal plane and vertical direction were (0.42, 0.26 m) for GPS-only, (0.42, 0.32 m) for GPS + GLONASS, and (0.40, 0.35 m) for the multi-constellation configuration. The corresponding overall 3D positioning errors of these three configurations were 50.0%, 43.8%, and 30.3%, respectively, lower than that of the standalone MADOCA-PPP.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 11

MADOCA-enhanced PPP/IMU: (a) horizontal error and solution distribution; (b) vertical error and solution distribution

The redundant information provided by the IMU enabled a more reliable and accurate computation of the velocity and attitude of the vehicle. The proposed system could precisely estimate the dynamics owing to the tactical-grade IMU. In particular, the velocity RMSEs in the east, north, and upward directions were 0.013, 0.016, and 0.018 m/s, respectively. Moreover, the attitude RMSEs of the roll, pitch, and heading components were 0.17°, 0.23°, and 0.24°, respectively. Limited improvement can be observed for different combinations of constellations. Table 4 shows the overall statistics for PPP/IMU in multiple tests. A 3D positioning error of 0.50 m can be achieved with the IMU-aided configurations. However, the advantages and improvements of the integration scheme might not be evident in a suburban scenario with sufficient observations and low susceptibility to environmental interference. Therefore, an experiment in a GNSS-challenging scenario was performed (Section 5.4) to further assess the performance of the MADOCA-PPP framework.

5.4 GNSS-Challenging Scenario

5.4.1 Standalone MADOCA-PPP

The second dynamic field experiment was conducted in a GNSS-challenging urban environment. In this experiment, the test duration and distance were approximately 40 min and 8.1 km, respectively. Figure 12 shows the vehicle trajectory in this scenario. The complex environment shown in the attached dashcam record can severely degrade the navigation performance because of NLOS or multipath signals. During the test, the vehicle was pulled over beside a tall building or under a shaded area for approximately 8 min. When standalone MADOCA-PPP was used, the number of visible satellites was frequently insufficient. Figure 12(b) shows the corresponding information. The average number of visible satellites and PDOP, presented as (average number of satellites, PDOP), were (5.93, 4.17) for GPS-only, (8.80, 2.44) for GPS + GLONASS, and (11.20, 1.99) for GPS + GLONASS + QZSS, respectively. As expected, the highest number of satellites was observed for the multi-constellation configuration. The lowest PDOP was also observed for this configuration, even in the highly challenging environment.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 12

Information for the second field experiment: (a) vehicle trajectory in the urban scenario; (b) number of visible satellites in view and PDOP for the urban scenario

The positioning error and CDF of the solution for the standalone MADOCA-PPP are shown in Figure 13. The errors were compared only for the epochs with sufficient satellites. However, some outliers affected the navigation capability. For example, GNSS signals were frequently interrupted because of environmental factors. The observations became unavailable when there were cycle slips in the ambiguities. A navigation filter was required to achieve convergence because the related covariance was reset. The dual- and multi-constellation configurations notably mitigated the degradation in the challenging environment. The 95th percentiles of solutions are presented in Figure 13. The CDFs of the solutions in the horizontal plane for the GPS-only, GPS + GLONASS, and GPS + GLONASS + QZSS configurations were 4.61, 5.11, and 1.78 m, respectively. The CDF of the multi-constellation configuration was improved by 61.4% and 65.2% compared with that of the GPS-only and GPS + GLONASS configurations, respectively, indicating that the multi-constellation configuration considerably improved the positioning capability in an urban environment. Specifically, the positioning RMSE values in the horizontal plane and vertical direction were (4.47, 6.20 m), (2.75, 4.82 m), and (1.03, 2.24 m) for the GPS-only, GPS + GLONASS, and GPS + GLONASS + QZSS configurations, respectively.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 13

Standalone MADOCA-PPP in a GNSS-challenging scenario: (a) horizontal error and solution distribution; (b) vertical error and solution distribution

The velocity could not be precisely determined with the standalone MADOCA-PPP regardless of the configuration. The velocity RMSEs in the east, north, and upward directions were 3.70, 4.33, and 0.41 m/s, respectively, for the GPS-only configuration; 2.87, 3.48, and 0.42 m/s, respectively, for the GPS + GLONASS configuration; and 1.76, 1.80, and 0.27 m/s, respectively, for the multi-constellation configuration. Although the error decreased as the number of constellations increased, the estimations were not sufficiently accurate in driving situations. Therefore, additional information was provided by the integrated scheme to accurately estimate the velocity.

5.4.2 MADOCA-Enhanced PPP/IMU

The results in the previous subsection revealed that an integration scheme is essential in a GNSS-challenging environment if the satellites visible for positioning are insufficient. A tightly coupled system can continuously propagate the navigation solution with fewer than four visible satellites, a scenario that frequently occurs in urban environments. Figure 14 presents the positioning errors of the constellations and the CDFs of their solutions in the horizontal plane and vertical direction, respectively. The IMU-aided system had considerably better results than the standalone MADOCA-PPP. Additionally, the solutions were smoother for all configurations with the IMU-aided system, whereas several spikes were observed in the standalone MADOCA-PPP solution. The performance degradation for an elapsed time of 1250–1750 s was attributed to pulling over in the extremely challenging region, as shown in Figure 12(a). The ZUPT was executed to reduce the error drift caused by the IMU. The effectiveness of this approach will be discussed in the next subsection.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 14

MADOCA-enhanced PPP/IMU in a GNSS-challenging scenario: (a) horizontal error and solution distribution; (b) vertical error and solution distribution

The 95th percentile of the horizontal solution in the multi-constellation configuration was 1.62 m. Greater satellite availability improved the navigation performance in the challenging environment. Particularly, the positioning RMSEs in the horizontal plane and vertical direction were (3.03, 3.92 m) for GPS-only, (1.60, 3.29 m) for GPS + GLONASS, and (0.92, 1.47 m) for GPS + GLONASS + QZSS, respectively. The 3D positioning RMSEs for the multi-constellation configuration were 65% and 52% lower than those for the GPS-only and GPS + GLONASS configurations, respectively.

The poor velocity estimation results obtained with standalone MADOCA-PPP can be mitigated by using the IMU. The IMU considerably increased the accuracy of the velocity determination, especially in GNSS-challenging environments. Specifically, when the IMU was incorporated, the velocity RMSEs in the east, north, and upward directions were 0.21, 0.19, and 0.48 m/s, respectively, for the GPS-only configuration; 0.07, 0.09, and 0.17 m/s, respectively, for the GPS + GLONASS configuration; and 0.02, 0.02, and 0.03 m/s, respectively, for the GPS + GLONASS + QZSS configuration. The precision achieved with these configurations was considerably higher than that achieved with standalone MADOCA-PPP; the IMU reduced the RMSEs from the order of meters to centimeters in each direction. These results indicate the importance of using an IMU and a multi-constellation configuration.

The overall statistical results for several evaluation metrics for each constellation configuration in the GNSS-challenging environment of the field experiment are presented in Table 5. The integration scheme clearly reduced the RMSEs for position. Although the standalone MADOCA-PPP performed well in conditions with high signal quality, it should be supplemented with an IMU and should employ a multi-constellation configuration for high performance in various environments.

View this table:

• View inline
• View popup

TABLE 5

Comparing the Position RMSE and Solution Distribution for Multiple Configurations in a GNSS-Challenging Scenario

5.4.3 Effectiveness of Motion Constraints

The field experiment was conducted to evaluate the proposed system under challenging conditions. Specifically, the experimental vehicle was parked under a shadow and by a tall building for approximately 500 s, as indicated in Figure 14. For most of this time, the vehicle was linked to no more than four satellites in the GPS-only configuration. In particular, the receiver on the car was connected to an insufficient number of satellites for 200 s of the 500-s duration. Furthermore, NLOS and multipath disruptions corrupted the GNSS signals. Performance is compromised if a measurement update cannot be performed because of corrupted or absent code and carrier-phase measurements; the IMU can only remain accurate for a short period in the absence of these updates. The resulting drift can eventually result in a large positioning error. Therefore, ZUPT is essential in these conditions. Figure 15 presents a comparison of the results obtained with and without ZUPT. The positioning RMSE values were (43.11, 36.31 m) without ZUPT. The position readings drifted radically owing to error accumulation in the IMU (Figure 15(a)). ZUPT greatly reduced these RMSEs to (3.03, 3.92 m), as indicated in Figure 15(b). In general, the motion constraint model significantly improved the results, particularly under the most demanding conditions.

• Download figure
• Open in new tab
• Download powerpoint

FIGURE 15

Comparison of the trajectory: (a) without ZUPT; (b) with ZUPT