The Design a TDCP-Smoothed GNSS / Odometer Integration Scheme with Vehicular-Motion Constraint and Robust Regression

: Global navigation satellite system (GNSS) is widely regarded as the primary positioning solution for intelligent transport system (ITS) applications. However, its performance could degrade, due to signal outages and faulty-signal contamination, including multipath and non-line-of-sight reception. Considering the limitation of the performance and computation loads in mass-produced automotive products, this research investigates the methods for enhancing GNSS-based solutions without signiﬁcantly increasing the cost for vehicular navigation system. In this study, the measurement technique of the odometer in modern vehicle designs is selected to integrate the GNSS information, without using an inertial navigation system. Three techniques are implemented to improve positioning accuracy; (a) Time-di ﬀ erenced carrier phase (TDCP) based ﬁlter: A state-augmented extended Kalman ﬁlter is designed to incorporate TDCP measurements for maximizing the e ﬀ ectiveness of phase-smoothing; (b) odometer-aided constraints: The aiding measurement from odometer utilizing forward speed with the lateral constraint enhances the state estimation; the information based on vehicular motion, comprising the zero-velocity constraint, fault detection and exclusion, and dead reckoning, maintains the stability of the positioning solution; (c) robust regression: A weighted-least-square based robust regression as a measurement-quality assessment is applied to adjust the weightings of the measurements adaptively. Experimental results in a GNSS-challenging environment indicate that, based on the single-point-positioning mode with an automotive-grade receiver, the combination of the proposed methods presented a root-mean-square error of 2.51 m, 3.63 m, 1.63 m, and 1.95 m for the horizontal, vertical, forward, and lateral directions, with improvements of 35.1%, 49.6%, 45.3%, and 21.1%, respectively. The statistical analysis exhibits 97.3% state estimation result in the horizontal direction for the percentage of epochs that had errors of less than 5 m, presenting that after the intervention of proposed methods, the positioning performance can fulﬁll the requirements for road level applications.


Introduction
To realize safe intelligent transportation systems (ITS), the demand for land vehicle navigation has been rapidly increasing. Global Navigation Satellite System (GNSS) is commonly used to provide global three-dimensional (3D) navigation with good long-term accuracy. It is currently the principal technique for providing globally referenced positioning without other aiding information. In kinematic reckoning (DR) aiding to enhance the performance without increasing the cost. Moreover, robust regression as an adaptive weighting method is combined with odometer-aided measurement, providing more reliable weighting indicators for state estimation in a scenario with increased multi-false measurements, without significantly burdening the navigation system with increased computation. This research investigates the feasibility of a cost-effective onboard solution for land-vehicle navigation applications without integrating INS and complex statistical models. Figure 1 illustrates the GNSS/odometer integration scheme with weighted least squares (WLS)-based robust regression, which includes the measurements implemented in this paper. In the proposed design, GNSS acts as a major sensor, comprising pseudorange, pseudorange rate, and TDCP measurements, for providing absolute positioning state (position and velocity). Whereas, the odometer with the forward speed acts as the constraint in the position and velocity domain. Those measurements will be processed by the WLS-based robust regression to reweight the scale indicators adaptively, and the state-augmented EKF integrates all of the inputs for the state estimation. In the following, the design of state-augmented EKF corresponding to delta measurement (TDCP) for estimating the position change is introduced in Section 2.1. Section 2.2 describes the measurement model, including conventional measurements and TDCP. Section 2.3 introduces the measurement model for the odometer, and the corresponding designs of vehicular-motion constraint and FDE mechanism in this study are also included. The robust regression based on the odometer aiding method proposed in this study is in Section 2.4. The process flow comprised of all of the proposed methods will be illustrated afterward.

Methods
Remote Sens. 2020, 12, x FOR PEER REVIEW 4 of 29 aided measurement, providing more reliable weighting indicators for state estimation in a scenario with increased multi-false measurements, without significantly burdening the navigation system with increased computation. This research investigates the feasibility of a cost-effective onboard solution for land-vehicle navigation applications without integrating INS and complex statistical models. Figure 1 illustrates the GNSS/odometer integration scheme with weighted least squares (WLS)based robust regression, which includes the measurements implemented in this paper. In the proposed design, GNSS acts as a major sensor, comprising pseudorange, pseudorange rate, and TDCP measurements, for providing absolute positioning state (position and velocity). Whereas, the odometer with the forward speed acts as the constraint in the position and velocity domain. Those measurements will be processed by the WLS-based robust regression to reweight the scale indicators adaptively, and the state-augmented EKF integrates all of the inputs for the state estimation. In the following, the design of state-augmented EKF corresponding to delta measurement (TDCP) for estimating the position change is introduced in Section 2.1. Section 2.2 describes the measurement model, including conventional measurements and TDCP. Section 2.3 introduces the measurement model for the odometer, and the corresponding designs of vehicular-motion constraint and FDE mechanism in this study are also included. The robust regression based on the odometer aiding method proposed in this study is in Section 2.4. The process flow comprised of all of the proposed methods will be illustrated afterward.

State-Augmented EKF Design
In this study, the position and velocity error states of GNSS positioning are estimated using EKF. Conventionally, the EKF design presents only the system error at a particular time. Whereas, a delta phase or delta position represents the integrated velocity over time [5]. To sufficiently use such delta measurement for estimating the change in the receiver position, the conventional system model should be modified and extended to include the estimated position and receiver clock bias from the last epoch as follows [5,37]:

State-Augmented EKF Design
In this study, the position and velocity error states of GNSS positioning are estimated using EKF. Conventionally, the EKF design presents only the system error at a particular time. Whereas, a delta phase or delta position represents the integrated velocity over time [5]. To sufficiently use such delta measurement for estimating the change in the receiver position, the conventional system model should be modified and extended to include the estimated position and receiver clock bias from the last epoch as follows [5,37]: δx = δr e u (t k ) δv e u (t k ) δr e u (t k−1 ) δb u (t k ) δb d (t k ) δb u (t k−1 ) T

1×12
(1) Remote Sens. 2020, 12, 2550 5 of 29 where δx is the error state vector including the position errors δr e u , velocity error δv e u , receiver clock bias δb u , and receiver clock drift δb d ; t k and t k−1 denote the corresponding epochs of the error states; each of δr e u and δv e u is a 3 × 1 vector, where the superscript e denotes a vector in the e-frame (earth-centered-earth-fixed coordinate system); Φ p 1 v and Q p 1 v denote the transition matrix (t k to t k+1 epoch) and the corresponding process noise matrix under assumed constant velocity conditions, for the current time; q p and q v are the spectral densities of the above process noise sequences; Φ t and Q t denote the transition matrix and process noise matrix for the receiver clock state; ∆t is the time interval between adjacent epochs; Φ p 1 vp 0 t and Q p 1 vp 0 t are the augmented models from Φ p 1 v and Q p 1 v to incorporate with the delta-measurement update; I is the identity matrix, corresponding to each component of the current position error state. The subscript p 1 , v, and p 0 denotes that the matrix contains the components corresponding to the current position, velocity, and previous position error state. This design both supports the previously defined dynamic system for the random walk model, and transfers the current position error vector to the previous spot in the error state vector during the propagation. In this design, the delta observables are directly related to the position increment as a constraint of the current position, relative to the previous position. Based on the phase-smoothing approach in the position domain, this EKF design can maximize the effectiveness of such measurements without assumptions about the dynamics of the vehicle [5] and control the growth of position errors. The correspondingly adopted measurement model in this study is described in the following subsection.

Measurement Model for Pseudorange, Doppler Shift, and TDCP
The models for GNSS measurements, including pseudorange, Doppler shift, and TDCP, which are related to the state of a vehicle are considered in the SPP mode in this study. The models of pseudorange and pseudorange rate measurement are well-documented in the literature, which can be found in [38]. The matrix form of the fundamental observations corresponding to the designed EKF can be expressed as, ρ i , R code/doppler , σ 2 0 , el i , and e i (t k ) denote the design matrices of the pseudorange and pseudorange rate, the observed pseudorange, geometric range between the receiver and the satellite, the observed pseudorange rate, geometric pseudorange rate, weighting matrix of the pseudorange/pseudorange rate measurement, a priori variance factor, the elevation angle of the i th satellite with respect to the rover, and LOS unit vector from the satellite to the receiver at t k , respectively. The design of the weighting matrix of pseudorange rate is similar to R code , but with a different σ 2 0 . The TDCP measurements, which are the delta-phase measurements from the same satellites between consecutive epochs, are conventionally used to estimate the average velocity [39][40][41][42][43]. In other words, the difference between the previous and current positions with mm-level noise can be observed through such phase differences. It could be partially observable if fewer than four satellites are continuously available in the absence of cycle slip. The advantage of this approach over conventional techniques is that carrier-phase measurements from the same satellites are merely available since the previous time epoch, rather than several seconds, to average out the code multipath error [5,37]. The original measurement equation of the carrier phase is expressed as follows: where λ, Φ i , and N i represent the wavelength, measured carrier phase in cycles, and integer ambiguity, respectively. The term η i includes a multipath/NLOS-contaminated source and receiver noise for the carrier phase. The difference between the carrier-phase measurements at consecutive epochs t k and t k−1 is described as follows: where ∆ indicates the differencing operation. Under the condition of no cycle slip, the integer ambiguity is eliminated by time differencing. Other common error sources, which include satellite clock bias, ephemeris error, tropospheric error, and ionospheric error, vary slowly within a sampling, and are, ∆δb i s , ∆δd i eph , ∆δd i iono , and ∆δd i trop are negligibly small. Notably, the transition between the adjacent ephemeris sets causes the discontinuity of ∆δd i eph , especially for ∆δb i s , when different sets are respectively selected between consecutive epochs for TDCP measurements [43]. Selecting the same ephemeris can prevent this problem during the transition. In transforming (10) into a relation of receiver position change, ∆ρ i has to be recombined as [43], where r e u and r e,i s are the receiver position vector and satellite position vector, respectively, including the x, y, and z components in the e-frame. The term ∆g reflects the orientation change of the LOS vector in the relative satellite-receiver geometry, and ∆D represents a change in the range and is proportional to the average Doppler shift caused by the relative satellite-receiver motion along the LOS vector [44]. The derivation of (11) is presented in detail in [39,43]. By substituting (11) into (10), the TDCP measurement equation becomes (12), and the corresponding measurement model is expressed by (13) and (14). The weighting matrix of the TDCP measurements, R TDCP , is also similar to the design of (17) with a different σ 2 0 .
The effectiveness of TDCP is affected by cycle slip. The cycle-slip detection method applied in this study is "Doppler integration" for a single-frequency receiver. This method uses an observed Doppler shift to predict the growth in a carrier phase observation as an indicator of a cycle slip, whose fundamental equation is presented in [38]. Nonetheless, due to the vehicular dynamics and the environment, the signal tracking quality of Doppler shift strongly affects the accuracy of detection. When a cycle slip occurs with false detection, new ambiguities are introduced in delta-phase measurements. On the other hand, the path delay of the carrier phase measurements due to NLOS reception could be similar to the pseudorange error, as mentioned previously in [8]. These errors could be introduced into carrier phase ambiguities and further degrade GNSS positioning if the phase-smoothing approach is applied. In this situation, standalone GNSS positioning cannot satisfy the requirement of accuracy and reliability in vehicular navigation. This study proposes the aiding from odometer velocity via OBD II to enhance the robustness of positioning in case of occurrence of such a situation or frequent cycle-slip. The corresponding aiding model based on the proposed EKF design is introduced in the following subsection.

Measurement Model for Odometer Observation with Vehicular-Motion Constraint
The odometer-derived velocity indicates that unless a vehicle is off the ground or slides on the ground, the velocity of the vehicle in the plane perpendicular to the forward direction is nearly zero [44,45]. Because an odometer provides the forward speed of a vehicle, the measurement is expressed as follows: where v v denotes the velocity measurement in the v-frame, and the lateral and vertical speeds are presumed to be zero. Because (15) need not be obtained in the e-frame, it cannot be used directly in the state estimation. In this study, based on the fact that v is the observation as the resultant speed Remote Sens. 2020, 12, 2550 8 of 29 of x, y, and z components for the velocity state in (1), the corresponding measurement model can be formulated as, with where H odo is the design matrix for the odometer; η odo comprises the measurement noise and abnormal sensing such as sliding on the ground; v e x , v e y , v e z are the initial GNSS velocity at the current epoch in the three directions; v 3d is the initial resultant speed. The determination of σ 2 odo can be given by a pre-calibrated process (e.g., comparison between the resultant speed of the reference system and the OBD II odometer). It is worth mentioning that this model only constrains the combined velocity for an estimation that is non-directional, and could face the local minimum problem during state estimation. Due to the lateral velocity of a land vehicle is nearly zero under normal conditions, it can act as an additional aiding source to enhance such state estimation. The relationship between the vehicle frame and the local-level frame (l-frame) in the horizontal direction is depicted in Figure 2, and the corresponding equation given by, where v lateral is the lateral velocity derived from the velocity state estimation; η lateral is the combination of velocity estimation error and heading error; H lateral is the design matrix of the lateral constraint based on the transformation of the e-frame, the l-frame, and vehicle frame; θ, ϕ, and Λ denote the vehicular heading, latitude, and longitude angles at t k−1 epoch. The determination of σ 2 lateral is based on the uncertainty of the estimated heading, which is an empirical coefficient in this study. θ is derived from the velocity state transformed into the l-frame and is thus vulnerable to the quality of GNSS measurement and significant vehicular motion. To secure the performance of the lateral constraint, the moving average filter (MAF) and the detection of heading condition are designed before (19)- (21) update EKF, which is formulated as, where θ i−1 , θ i , dθ thre , γ, and β denote the averaged heading angle at i − 1 epoch, heading angle derived from the estimated velocity at the current epoch, threshold of heading change between the current heading and the previously averaged angle, window size, and the number of θ i used in MAF. The implementation of the proposed lateral velocity constraint in (19)- (22) further enhances the directivity of smoothing between adjacent epochs, because a priori information is used. In this design, when the heading change is less than dθ thre , and β exceeds γ, the lateral constraint will use the averaged heading based on the outcome of (22); otherwise, the lateral constraint will be closed until the counts of qualified θ i reach γ. The intervention of dθ thre is used to assure that the designed constraint will not be activated during significant vehicular motion, such as a turning maneuver.
Notably, the calculation of angular distance in (22) must consider the vibration of angular conversion from the estimated velocity between adjacent epochs, such as the case of θ i − θ i−1 > 180 • or < −180 • , averting a false output/determination in (22).
formulated as follows: where ∆ , , and indicate the displacement between adjacent epochs in the lateral direction of the vehicle, scale coefficient, and instantaneous velocity from an odometer at the last GNSS epoch, respectively. As shown in Figure 2, the concept of (23) is based on the fact that, theoretically, for a land vehicle, the displacement in the lateral direction between consecutive GNSS epochs will not surpass the resultant speed in the forward. While, is a variable coefficient depending on different vehicular conditions, such as low or high dynamics, coping with various scenarios in urban environments. If an unsatisfactory state estimation is detected and excluded by (23) or a GNSS outage occurs, the information, including odometer-based velocity and smoothedheading from the last epoch can provide the DR solution. Although this solution provides only a 2-D solution with a constant height state, which is vulnerable to variations in height and heading, such design can bridge the gap in positioning and maintain the estimation accuracy over a short time. Another additional aid from an odometer is the zero-velocity detection. When a vehicle is stationary, the velocity state estimation is nearly zero. However, in the densely packed urban areas, NLOS/multipath-contaminated reception becomes more severe for the duration that the vehicle is stationary. Such continuous faulty-signal reception not only degrades the velocity estimation but also misleads the position state due to the path delay. These measurements falsely update the EKF during the static motion, causing the EKF to diverge. To address this situation, the use of an odometer can provide precise constraints by assigning the velocity to be zero in any direction, thereby preventing error growth. This zero-velocity constraint is expressed as, where the identity matrix corresponds to each component of the velocity state in (1). Moreover, during the activation of zero-velocity constraint, GNSS measurements will not update EKF and fix the position state and heading output, thereby preventing an erroneous position change. This position-fix design can be formulated in (26), where denote the velocity threshold of determining the stationary motion. This mechanism can maintain stability for the state and its The advantage of the odometer-based methods mentioned above is that they provide additional measurements (constraints) and improve the state estimation result under the contamination of multipath/NLOS reception. However, if a limited number of tracked satellites with bad geometry exists and false GNSS measurements update EKF, unsatisfactory estimation could still emerge. Using the similar concept of (19)- (22), odometer-based FDE can be implemented. The velocity, based on an odometer, is a forward speed in the v-frame, as given in (15), following the forward direction of the vehicle. If the result of the position state estimation from GNSS presents an abnormal motion in the lateral, the estimation result is dominated by unhealthy signals. Such an FDE method can be formulated as follows: where ∆d lateral , s, and v k−1 indicate the displacement between adjacent epochs in the lateral direction of the vehicle, scale coefficient, and instantaneous velocity from an odometer at the last GNSS epoch, respectively. As shown in Figure 2, the concept of (23) is based on the fact that, theoretically, for a land vehicle, the displacement in the lateral direction between consecutive GNSS epochs will not surpass the resultant speed in the forward. While, s is a variable coefficient depending on different vehicular conditions, such as low or high dynamics, coping with various scenarios in urban environments.
If an unsatisfactory state estimation is detected and excluded by (23) or a GNSS outage occurs, the information, including odometer-based velocity and smoothed-heading from the last epoch can provide the DR solution. Although this solution provides only a 2-D solution with a constant height state, which is vulnerable to variations in height and heading, such design can bridge the gap in positioning and maintain the estimation accuracy over a short time.
Another additional aid from an odometer is the zero-velocity detection. When a vehicle is stationary, the velocity state estimation is nearly zero. However, in the densely packed urban areas, NLOS/multipath-contaminated reception becomes more severe for the duration that the vehicle is stationary. Such continuous faulty-signal reception not only degrades the velocity estimation but also misleads the position state due to the path delay. These measurements falsely update the EKF during the static motion, causing the EKF to diverge. To address this situation, the use of an odometer can provide precise constraints by assigning the velocity to be zero in any direction, thereby preventing error growth. This zero-velocity constraint is expressed as, where the identity matrix I corresponds to each component of the velocity state in (1). Moreover, during the activation of zero-velocity constraint, GNSS measurements will not update EKF and fix the position state and heading output, thereby preventing an erroneous position change. This position-fix design can be formulated in (26), where v thre denote the velocity threshold of determining the stationary motion. This mechanism can maintain stability for the state and its covariance matrix of EKF, and the accuracy of the solution. The flowchart depicting the utilization of zero-velocity constraints is presented next (Figure 3):

Robust Regression with Odometer Aiding
WLS estimation is more susceptible to outliers compared with KF-based estimation. Robust regression, based on WLS residuals analysis, can adjust the weight measurements probably affected by multipath/NLOS reception whose corresponding residuals are large [18,19]. Such a process is a generalization of the maximum likelihood estimator (M-estimator), and is applied in this study. Considering that the state estimation for WLS is an iterative process (until the convergence is reached), the weightings depending on the residuals are adaptively reweighted in each iteration. The weighting function of this robust regression in each iteration can be expressed as follows [18,19]: where w i j denotes the diagonal element of the weighting matrix corresponding to the i th satellite measurement at the j th iteration;σ 0 is the a posteriori scale factor given by median absolute deviation due to the spread of the residuals; α is the tuning constant and is set to 1.345 based on [46]; α/ r i j /σ 0 is the weighting scale. Using (27), the residuals are considered as an approximation of the measurement errors, and the results until WLS convergence are implemented to adaptively adjust the measurement covariance matrix used in EKF estimation [19]. A potential limitation of this method emerges if abnormal measurements are in majority, as mentioned previously. To improve the performance of robust regression under the emergence of multiple faulty measurements, this study proposes the use of accumulated distance based on the odometer information as a positional-displacement constraint to the WLS filter during the robust regression operation. When the vehicle does not change its heading orientation significantly, the moving distance from the wheel is similar to the difference between the corresponding position states of the adjacent epochs. Considering the error accumulation, due to the non-linear acceleration and the driving path challenges this assumption, such distance information is used only as an aided constraint for obtaining a weighting scale in (27), rather than as a direct measurement update in EKF estimation. The wheel information from the OBD II interface is transformed into instantaneous velocity; therefore, the accumulated distance can be calculated as follows, where ∆d 3d , ξ, and v κ indicate the averaged displacement between consecutive GNSS epochs, odometer rating between GNSS consecutive epochs, and instantaneous velocity measurement through the OBD II interface. Similar to (16)- (18), this measurement model of position constraint for the WLS filter can be formulated as, with r 3d = (∆r e x ) 2 + (∆r e y ) 2 + (∆r e z ) 2 where H WLS , δx WLS , η ∆d , and r 3d are the design matrix of averaged displacement, error state, including position and receiver clock bias for WLS filter, combined errors comprising sliding, heading change and noise, and initial distance derived from the change of initial rover position between consecutive epochs; ∆r e x , ∆r e y , and ∆r e z denote the difference between the initial position vectors in the three directions of the e-frame. The measurement model in the WLS filter is similar to the design in EKF, as mentioned previously, and is not differentiated here. The benefits of odometer aiding for robust regression, include the improvement in observability of geometry, which combines the ground-based (odometer) and space-based (GNSS) sources, and the enhanced connection of relative-position change in addition to the TDCP measurements. Those contributions enhance WLS filter to find more accurate results; therefore, a priori weighting indicators in (27) will be more reliable to implement in EKF estimation.
The flowchart depicting the combination of the proposed methods in this study is shown in Figure 3. The data flow from process 1 to process 6 aims to determine the condition of combining TDCP measurements to update EKF (process 6) according to the result of receiving carrier-phase measurements (process 3) and cycle slip detection (process 4). After odometer measurements are introduced, the priority is the detection of stationary motion (process 8), since the strong intervention from zero velocity constraint updates EKF without any GNSS measurements (process 9 and 13) and fix the position output (process 14), ending the whole process flow for the current epoch. While the vehicle is not stationary, odometer-aided robust regression process based on WLS will be conducted to obtain weighting scales (process 11), and the corresponding GNSS measurements, odometer-aided velocity constraint (based on the resultant speed, process 10), and lateral velocity constraint (process 12) will be used to update EKF accordingly (process 13). After the position and velocity states have been estimated, and the heading (derived from the velocity states) does not change significantly (process 15), odometer-based FDE in (23) will be activated (process 16). Once an unsatisfactory estimation is detected, the odometer-based DR solution will bridge the gap during the estimation (process 17). Otherwise, MAF will be activated (process 18) and produce the averaged heading for other processes (process 12, 15, and 17) for the next epoch. 12) will be used to update EKF accordingly (process 13). After the position and velocity states have been estimated, and the heading (derived from the velocity states) does not change significantly (process 15), odometer-based FDE in (23) will be activated (process 16). Once an unsatisfactory estimation is detected, the odometer-based DR solution will bridge the gap during the estimation (process 17). Otherwise, MAF will be activated (process 18) and produce the averaged heading for other processes (process 12, 15, and 17) for the next epoch.

Experiment
In the experiment, the reference system uses a high tactical-grade INS (NovAtel SPAN ® LCI, Canada) with a multi-band reception antenna (NovAtel GPS-GGG-703-HV), and the evaluated system is an automotive-grade Allystar EVK-2024 module with a single-band active patch antenna (Allystar AGR 6301, China); these sensors are mounted on the same vehicle platform (Figure 4). Table  1 lists the specifications of the reference system. Table 2 shows the performance of the OBD II odometer applied in this study, which is evaluated from the other individual tests in an urban area.

Experiment
In the experiment, the reference system uses a high tactical-grade INS (NovAtel SPAN ® LCI, Canada) with a multi-band reception antenna (NovAtel GPS-GGG-703-HV), and the evaluated system is an automotive-grade Allystar EVK-2024 module with a single-band active patch antenna (Allystar AGR 6301, China); these sensors are mounted on the same vehicle platform (Figure 4). Table 1 lists the specifications of the reference system. Table 2 shows the performance of the OBD II odometer applied in this study, which is evaluated from the other individual tests in an urban area. To evaluate the performance with the proposed methods, an experimental route (red trajectory) is designed under a GNSS-challenging environment (Tainan city), as shown in Figure 5. The referenced position solution from tactical-grade INS is based on post-smoothing processing with full constellations and frequencies receiving, whereas the evaluated EVK-2024 module uses two satellite constellations (GPS, Beidou) with single-frequency (L1 band) reception in the SPP mode. All of the results are synchronized with GPS time, and the relationship between the reference system and the EVK-2024 antenna (as shown in Figure 4) are measured by the surveying tapes. The performance analysis can be implemented after the referenced results are projected to the evaluated antenna.     The total experiment duration is approximately 40 min. The tested results are based on forwardfilter estimation with the proposed EKF design and without any post-smoothing processing. The scenarios A-E are selected to discuss the contribution of the proposed methods, where are multipath/NLOS-contaminated environments, and their locations are represented by white blocks in Figure 5. The corresponding street-view images of these areas are shown in Figure 6a-f. Moreover, as shown in Figure 7a,b, respectively, the position dilution of precision (PDOP) increases, and the number of visible satellites significantly decreases after the vehicle approaches urbanized areas, although two constellations are used. The observability of healthy satellite signals is frequently unsatisfactory during the operation. For analysis, the implemented methods are classified into proposed odometer-aided constraints and robust regression, i.e., method A, and method B, respectively. The tested results include code-based measurements only (conventional GNSS), code/TDCP measurements (TDCP-aided GNSS), code/TDCP measurements with odometer-aided constraints (TDCP-aided GNSS with method A), code/TDCP measurements with robust regression (TDCP-aided GNSS with method B), and code/TDCP measurements with both methods (TDCP-aided GNSS with both methods). The evaluated solutions are colored in cyan, orange, blue, green, and black, in order for the following analysis for the trajectory and the line chart of errors. Moreover, to evaluate the performance of these solutions during a GNSS outage (area I of scenario E, Figure 6e), the solutions without method A will use assumptions, such as constant velocity from the last epoch to bridge the gap during GNSS positioning. Whereas, the others with method A perform odometeraided DR in this experiment. For the error analysis, the results from the reference system at the first epoch in this experiment is selected as the origin of the l-frame, and all of the estimated results are transformed from the e-frame into the l-frame accordingly in the following. The total experiment duration is approximately 40 min. The tested results are based on forward-filter estimation with the proposed EKF design and without any post-smoothing processing. The scenarios A-E are selected to discuss the contribution of the proposed methods, where are multipath/NLOS-contaminated environments, and their locations are represented by white blocks in Figure 5. The corresponding street-view images of these areas are shown in Figure 6a-f. Moreover, as shown in Figure 7a,b, respectively, the position dilution of precision (PDOP) increases, and the number of visible satellites significantly decreases after the vehicle approaches urbanized areas, although two constellations are used. The observability of healthy satellite signals is frequently unsatisfactory during the operation. For analysis, the implemented methods are classified into proposed odometer-aided constraints and robust regression, i.e., method A, and method B, respectively. The tested results include code-based measurements only (conventional GNSS), code/TDCP measurements (TDCP-aided GNSS), code/TDCP measurements with odometer-aided constraints (TDCP-aided GNSS with method A), code/TDCP measurements with robust regression (TDCP-aided GNSS with method B), and code/TDCP measurements with both methods (TDCP-aided GNSS with both methods). The evaluated solutions are colored in cyan, orange, blue, green, and black, in order for the following analysis for the trajectory and the line chart of errors. Moreover, to evaluate the performance of these solutions during a GNSS outage (area I of scenario E, Figure 6e), the solutions without method A will use assumptions, such as constant velocity from the last epoch to bridge the gap during GNSS positioning. Whereas, the others with method A perform odometer-aided DR in this experiment. For the error analysis, the results from the reference system at the first epoch in this experiment is selected as the origin of the l-frame, and all of the estimated results are transformed from the e-frame into the l-frame accordingly in the following.

Results and Discussion
In this section, the analysis of the proposed methods will be categorized into separate subsections. Figures 8-12 show the trajectory and the corresponding positioning errors in the E, N, and U directions of the l-frame for the discussed scenarios. In addition, the statistical analyses are shown in Tables 3-6 accordingly.

Performance Analysis of TDCP-Aided GNSS
The contribution of TDCP aiding could be mainly found in a multipath/NLOS-contaminated environment. It can be shown in Figure 8 (scenario A, time 520-625 s), the TDCP-aided GNSS solution produced more accurate results than the conventional GNSS solution. Pseudorange measurements, suffering from such contamination, presented unsatisfactory estimation with approximately 20 m of

Results and Discussion
In this section, the analysis of the proposed methods will be categorized into separate subsections. Figures 8-12 show the trajectory and the corresponding positioning errors in the E, N, and U directions of the l-frame for the discussed scenarios. In addition, the statistical analyses are shown in Tables 3-6 accordingly.

Performance Analysis of TDCP-Aided GNSS
The contribution of TDCP aiding could be mainly found in a multipath/NLOS-contaminated environment. It can be shown in Figure 8  Most of the time, the use of TDCP measurements with state-augmented EKF lessened the multipath effect; however, it is similar to other phase-smoothing approaches that lead to unsatisfactory results. For example, in Figure 9  Taking scenario B as an example, despite using the TDCP measurements, the maximum error was still up to 15 m in height, and the maximum error is 13 m in the E-direction, which is similar to the conventional GNSS. The reason is that the multipath/NLOS reception and the limited satellite geometry result in a deviated estimation in these challenging environments.
In addition to the influence of multipath/NLOS signals, the potential problem of using the phase-smoothing approach would be the occurrence of cycle-slip detection. When unstable tracking qualities frequently occur for Doppler-shift measurements, the cycle-slip detection generates false outputs, and range errors due to the path delay are intervened, further degrading the state estimation. This situation can be found in Figures 9 and 10 (scenarios B and C, during time 872-875 s and 975 s), where TDCP-based GNSS presented a degraded position in the horizontal direction. Even worse, this degradation of state estimation could become severe in a frequent GNSS-outage environment. In Figure 12

Performance Analysis of TDCP-Aided GNSS with Method A
Under an open-sky condition (e.g., PDOP value was less than four, such as Figure 11 (scenario D)), the results of TDCP-aided GNSS with/without method A were acceptable. It indicates that phasesmoothing is a critical factor in obtaining a reliable result. However, in constrained environments, due to enhanced positional connection between the adjacent epochs, TDCP-aided GNSS plus method A exhibited a more reliable outcome (blue and black lines) compared with the others without method A. For instance, in Figure 9 (scenario B, time 815 s), when the phase-smoothing approach was unavailable, the use of method A was able to improve the position state estimation, reducing the Edirection error by around 3.60 m when comparing with TDCP-based GNSS. Moreover, for the situation of deteriorated results due to the degraded cycle-slip detection (during time 872-875 s and at time 975 s), in these signal-degraded environments, the positioning errors produced with method A mitigated from 7.82 to 3.58 m (at time 873 s) and 5.53 to 2.38 m (at time 975 s) in the horizontal direction, being able to reach the SPP-nominal level. A similar case could be found in Figure 10 (scenario C, during time 963-969 s).

Performance Analysis of TDCP-Aided GNSS with Method A
Under an open-sky condition (e.g., PDOP value was less than four, such as Figure 11 (scenario D)), the results of TDCP-aided GNSS with/without method A were acceptable. It indicates that phase-smoothing is a critical factor in obtaining a reliable result. However, in constrained environments, due to enhanced positional connection between the adjacent epochs, TDCP-aided GNSS plus method A exhibited a more reliable outcome (blue and black lines) compared with the others without method A. For instance, in Figure 9 (scenario B, time 815 s), when the phase-smoothing approach was unavailable, the use of method A was able to improve the position state estimation, reducing the E-direction error by around 3.60 m when comparing with TDCP-based GNSS. Moreover, for the situation of deteriorated results due to the degraded cycle-slip detection (during time 872-875 s and at time 975 s), in these signal-degraded environments, the positioning errors produced with method A mitigated from 7.82 to 3.58 m (at time 873 s) and 5.53 to 2.38 m (at time 975 s) in the horizontal direction, being able to reach the SPP-nominal level. A similar case could be found in Figure 10   With the help of using method A, the benefit of the proposed zero-velocity constraint with position-fix design in (24)-(26) was significantly presented when the vehicle was stationary. Due to continuously faulty signal reception, the major contribution of method A could include the precise detection for vehicle speed, assuring the timing of stoppage, the duration time, and preventing the estimated position from presenting a deteriorated change. As shown in  With the help of using method A, the benefit of the proposed zero-velocity constraint with position-fix design in (24)-(26) was significantly presented when the vehicle was stationary. Due to continuously faulty signal reception, the major contribution of method A could include the precise detection for vehicle speed, assuring the timing of stoppage, the duration time, and preventing the estimated position from presenting a deteriorated change. As shown in Figures 8-11 (scenarios A-D, during time 550-610 s, 835-863 s, 907-954 s, and 1636-1663 s, indicated by yellow circles), the solution with method A yielded a more accurate outcome. The aid of method A maintains the position output that is the same as the position state before transiting to the static, averting the use of unhealthy measurements. Table 3 shows the related analysis for the stationary condition in the experiment. The position accuracy in this situation improved because of method A, and the position accuracy in 3D attained 4.87 m with an improvement of 41.7%, demonstrating that the zero-velocity constraint with position-fix design successfully enhanced the performance in the challenging circumstance.
For method A, the odometer provides only a relative constraint in positional displacement, and this study uses the constraint design in (16)- (22) in the velocity domain to improve this constraint, as mentioned earlier. However, this constraint is along the forward direction of the vehicle, meaning that its aid is less influential in the lateral direction, especially for the case where the estimated heading is degraded. Table 4 shows that method A improves the positional accuracy in the forward (RMSE 1.87 m), but has limited effectiveness in the lateral (RMSE 2.10 m, similar to the TDCP-based solution without method A). A representative situation was found in Figure 8 (scenario A, during time 519-547 s), where the solution with method A still presented a deteriorated trajectory misled to the oncoming road (indicated by the white dotted line). In addition, there is still a potential risk on the position-fix design for method A. In Figure 8 (scenario A, during time 550-610 s), the height error based on the scheme with method A presented an error of 12.12 m during the static motion. In this representative case, because the result of position estimation had insufficient time to revert to the stable state before transiting to the stationary, this degraded result was selected to be the position-fix output. Although, the velocity error was corrected back to the stationary state, and averted to influence the state estimation further as the vehicle started to move, this erroneous position result was still unable to adjust until the end of stationary. Moreover, even if method A was implemented, an unsatisfactory result was still observed in the previously discussed GNSS-outage environment (area I of scenario E, during time 1945-1954 s). The proposed odometer-aided FDE did not detect the false result, because the lateral errors during time 1947-1954 s did not surpass the threshold. Specifically, the effectiveness of the detection is vulnerable to the quality of the initial position, such as an inaccuracy result at time 1946 s. Although the result was relatively stable than the GNSS solutions without method A in this area (cyan and orange result), the maximum error presented an unsatisfactory outcome (9.94 m and 5.36 m in the horizontal and vertical directions) due to the limited effectiveness of method A.

Performance Analysis of TDCP-Aided GNSS with Method B
The contribution of method B is to adjust the measurement weights adaptively. In this case, the influence of inconsistent measurements can be reduced, and the estimated position becomes closer to the designated path in comparison with TDCP-based GNSS. For instance, In Figures 8 and 11 (scenario A and D, during time 525-540 s and 877-884 s), the solution with method B (green and black lines) showed superior robustness during the operation in challenging environments, whereas the solutions without method B misdirected the vehicular path towards an oncoming-traffic road. In Figure 9 (scenario B, time 815 s), the positional distortion due to the multipath/NLOS contamination was improved after the intervention of method B, even when the TDCP-aided approach was inactive.
However, unsatisfactory estimation was still found in static situations. Although the performance of the TDCP-based solution with method B in the stationary condition was superior to the typical solutions without it (cyan and orange lines), the present finding exhibits that the remaining error was still up to 13.47 m in the horizontal direction (scenario B, time 853 s) and 10.62 m in the vertical direction (scenario C, time 932 s) for the harsh cases. This situation was attributed to the increase in the number of abnormal signals for the stationary condition. Another discussed case was in the GNSS-outage environment, as shown in Figure 12 (area I of scenario E, during time 1945-1954 s). In this area, the aid of method B was not obvious, where the positional error reached a maximum of up to 8.95 m and 5.43 m in the horizontal and the vertical directions. That was because method B failed to find appropriate scales for reweighting. Taking the situation of time 1948 for example (the corresponding information of measurements as shown in Table 5), it is seen that the range errors of five observations were over 10 m, even reaching 37 m. In this situation, for method B, residuals could not adequately reflect the approximation of measurement errors due to the influence of the contaminated signals, and some false measurements such as G31, G32, and C10 listed in Table 5 still significantly influenced the estimation. 1 "G" and "C" indicate GPS and Beidou satellite system, respectively. 2 Range measurements were corrected by the process of broadcast ephemeris. Moreover, considering that WLS residuals are derived from the estimated states (including the receiver clock bias), range measurements in this table are also adjusted by the estimated receiver bias for consistency. 3 Range error means the difference between range measurements and true ranges. 4 The definitions correspond to (27).

Performance Analysis of TDCP-Aided GNSS with Both Methods
The proposed combination of method A and method B presents complementary attributes of both the methods and outperforms the standalone schemes of each aiding approach. For a multipath/NLOS-contaminated environment, the improvement of the absolute position estimation is mainly contributed by method B. This situation can be obviously found in Figures 8 and 11 (scenario A and D), where the combined result showed that the trend of the black line was similar to that of the green line. Another contribution can also be observed under the stationary condition when method A and method B are adopted at the same time. The previous situations can be used as examples. In Figure 8  1636-1663 s), this combined scheme mitigated the limitation of method B and prevented the related error growth, and the maximum of the remaining error for these areas was reduced by approximately 8.93 m, 5.97 m, and 2.91 m in the horizontal direction, and 2.71 m, 8.04 m, and 3.99 m in the vertical direction, respectively. According to the previous results, it can be found that: (1) method B provides improved position before the dynamical motion transits to stationary, and therefore, method A can provide a superior position-fix result. (2) method A effectively averts the degradation of continuous reception of abnormal signals (where method B is limited) and maintains the performance under stationary motion. According to the statistical analysis (Table 3), the TDCP-based GNSS with both methods showed an RMSE of 4.06 m in 3D with an improvement of 51.4% under the static condition. The maximum error in 3D improved from 37.24 m (conventional solution) to 8.88 m, which is significantly better than the performance based on method A only (12.45 m).
The other significant advantage of the combined method is in the frequent GNSS-outage environment. In Figure 12 (area I of scenario E, during time 1945-1954 s), the unsatisfactory performance was improved by 63.2% and 82.3% at most in the horizontal and vertical directions in contrast with method B. The main reasons are the following mechanisms from both methods in this area: (1) odometer-aided robust regression, providing a more reliable scale for reweighting and more accurate position state (compared with method B only) before driving into the GNSS-degraded environment (epochs before time 1948 s); (2) odometer-based FDE, averting unsatisfactory position state estimation due to the majority of faulty signals; (3) odometer-aided DR, bridging the duration of GNSS outage. The combined performance significantly increased the availability of satisfactory results in such a harsh situation. Although the degradation of state estimation still emerged during time 1963-1967 s, where the vehicle was performing a turning motion, and the auxiliary constraint and FDE in (19)-(23) were deactivated due to the continuous heading change, the remaining aiding methods still provided better performance compared with others, as shown in the trajectory (Figure 12).
Based on this proposed solution, the overall analysis in Table 6 presents an RMSE of 2.51 m and 3.63 m in the horizontal and vertical directions with improvements of 35.1%, and 49.6%, respectively. The indicators of the percentage of epochs, including errors of less than 3 m (narrow-road accuracy) and 5 m (wide-road accuracy), were 86.1% and 97.3% in the horizontal direction and 63.9% and 84.4% in the vertical direction for these two levels, respectively. The corresponding analysis of the empirical cumulative distribution function (CDF) is shown in Figure 13. An analysis of vehicular tracking in Table 4 shows that the proposed combination scheme reduced the forward error and lateral error with an RMSE of 1.63 m and 1.95 m, corresponding to improvements of 45.3% and 21.1%. Statistically, this finding shows that, by using a low-cost automotive-grade receiver, the proposed GNSS/odometer scheme and the combined methods without any GNSS correction services can satisfy the accuracy requirement for locating a vehicle's position on a designated road in GNSS-challenging environments.

Conclusions
This study investigated GNSS/odometer integration with TDCP-designed EKF by utilizing vehicular motion constraint and robust regression for an automotive-grade receiver in an urban environment. The findings show that the phase-smoothing technique improved the state estimation and smoothed out some multipath effects after the TDCP updated the EKF state, but was subject to cycle-slip and NLOS reception. The degradation was apparent during the stationary motion. After the odometer-based aiding (method A) intervened, the aided constraints based on the resultant velocity and zero-assumption velocity of the lateral significantly enhanced the performance, when the state estimation reverted to the unsmoothed level due to a cycle slip. Moreover, the additional aiding from method A, including the zero-velocity constraint and the position-fix design, effectively

Conclusions
This study investigated GNSS/odometer integration with TDCP-designed EKF by utilizing vehicular motion constraint and robust regression for an automotive-grade receiver in an urban environment. The findings show that the phase-smoothing technique improved the state estimation and smoothed out some multipath effects after the TDCP updated the EKF state, but was subject to cycle-slip and NLOS reception. The degradation was apparent during the stationary motion. After the odometer-based aiding (method A) intervened, the aided constraints based on the resultant velocity and zero-assumption velocity of the lateral significantly enhanced the performance, when the state estimation reverted to the unsmoothed level due to a cycle slip. Moreover, the additional aiding from method A, including the zero-velocity constraint and the position-fix design, effectively averted the contamination of abnormal signals when the vehicle was stationary. However, the results show that method A had limited effectiveness in the lateral direction and presented a risk of fixing the position at an unsatisfactory level. When robust regression (method B) was added, the combination of method A and method B not only improved the accuracy in the lateral direction and the stationary condition, but also apparently enhanced the contribution of the odometer-aided FDE and DR, which effectively avoided erroneous outputs from the degradation of state estimation and performed seamless positioning, with adequate performance over a short period during a GNSS outage.
The experimental results show that the GNSS/odometer integration with the proposed methods can satisfy the demand of guiding a vehicle on designated roads with a single-frequency SPP mode in an urban area, without INS aiding. However, the performance of the proposed methods significantly depends on the empirical coefficients by the tuning process, including (21)- (23). An adequate setting for general cases is a time-consuming process. The auxiliary mechanism can include machine learning approaches, such as the random forest algorithm, for improving this situation in the future. Moreover, to mitigate the limitations of the proposed odometer-aided FDE and DR during a heading-change motion, such as a turning maneuver, the aid of z-channel gyro is considerable, and must be investigated in the future.