An innovation-based cycle-slip, multipath estimation, detection and mitigation method for tightly coupled GNSS/INS/Vision navigation in urban areas
Bo Xu Shoujian Zhang Jingrong Wang Jiancheng Li
















An innovation-based cycle-slip, multipath estimation, detection and mitigation method for tightly coupled GNSS/INS/Vision navigation in urban areas
    Bo Xu, Shoujian Zhang, Jingrong Wang, Jiancheng Li
Bo Xu and  Shoujian Zhang are with School of Geodesy and Geomatics, Wuhan University, Wuhan 430079, China; Corresponding author: Shoujian Zhang, Email: shjzhang@sgg.whu.edu.cn

Jingrong Wang is with the GNSS research center, Wuhan University, Wuhan 430079, China;
Jiancheng Li is with the School of Geodesy and Geomatics, Hubei Luojia Laboratory, Wuhan University, Wuhan 430079, China;
    September 9, 2024
========================================================================================================================================================================================================================================================================================================================================================================================================================================================



§ ABSTRACT 
Precise, consistent, and reliable positioning is crucial for a multitude of uses. In order to achieve high precision global positioning services, multi-sensor fusion techniques, such as the Global Navigation Satellite System (GNSS)/Inertial Navigation System (INS)/Vision integration system,
combine the strengths of various sensors. This technique is essential for localization in complex environments and has been widely used in the mass market. 
However, frequent signal deterioration and blocking 
in urban environments exacerbates the degradation of GNSS positioning and negatively impacts the performance of the multi-sensor integration system. For GNSS pseudorange and carrier phase observation data in the urban environment, we offer an innovation-based cycle slip/multipath estimation, detection, and mitigation (I-EDM) method to reduce the influence of multipath effects and cycle slips on location induced by obstruction in urban settings. The method obtains the innovations of GNSS observations with the cluster analysis method. Then the innovations are used to detect the cycle slips and multipath. Compared with the residual-based method, the innovation-based method avoids the residual overfitting caused by the least square method, resulting in better detection of outliers within the GNSS observations. The vehicle tests carried out in urban settings verify the proposed
approach. Experimental results indicate that the accuracy of 0.23m, 0.11m, and 0.31m in the east, north and up components can be achieved by the GNSS/INS/Vision tightly coupled system with the I-EDM method, which has a maximum of 21.6% improvement when compared with the residual-based EDM (R-EDM) method.  














multi-sensor fusion; tightly coupled integration; cycle slip detection; multipath mitigation; innovation-based method. 




§ INTRODUCTION

Autonomous driving has captured significant attention in both research and practical applications, where GNSS serves as the primary positioning technology to provide continuous, accurate, and stable positioning results <cit.>. Despite the worldwide prevalence and extensive utility of GNSS, its signals are susceptible to disruptions such as multipath effects and signal loss in complex urban settings, resulting in decreased positioning accuracy <cit.>. Consequently, the utilization of multi-sensors provides a more comprehensive and robust positioning capability, which becomes the optimal choice for pose estimation in urban scenarios <cit.>.


Integrating different sensors, such as INS and cameras, with GNSS can fully utilize the locally accurate characteristics of INS/Visual Inertial Odometry (VIO) and the global drift-free characteristics of GNSS. Fusing GNSS/INS has been widely studied. Significant improvements in accuracy, continuity, and reliability of localization are carefully analyzed and evaluated  <cit.>. However, in an integrated system where GNSS serves as the cornerstone for providing global positioning, the positioning accuracy is compromised or even worsens due to the rapid error accumulation of the microelectromechanical system inertial measurement unit (MEMS-IMU) when GNSS observations are affected by multipath effects or lost tracking  <cit.>.  


To mitigate the drifts of MEMS-IMU, the low-cost visual camera is applied in GNSS/INS integration system to offer redundant measurements. The multi-sensor fusion extended Kalman filter (MSF-EKF) framework proposed by Lynen et al. <cit.> allows for the seamless loose sensor-feed integration of GNSS, INS, and visual sensors. A semi-tightly coupled framework of multi-GNSS/INS/Vision based on graph optimization is proposed by Li et al. <cit.>. The method produces consistent and accurate global positioning outputs and operates well in GNSS-challenged environments. A tightly coupled GNSS/INS/Vision system with open-source code is proposed by Cao et al. <cit.>; However, the method only incorporates Doppler shift and code pseudorange measurements. To enhance the navigation accuracy during GNSS outages, Liao et al. <cit.> propose a tightly coupled framework to make full use of measurements from real-time kinematic (RTK)/INS/Vision. The method improves the accuracy and robustness of the system under conditions where GNSS is unavailable.  However, these methods regard GNSS observations as equal weight, without considering the impact of outliers on the system, thereby failing to make the most of the GNSS observations.

Effective GNSS measurement quality control is key to ensuring that the estimator is not affected by outliers in urban situations. A practical quality control technique is to assign weights according to the quality of the signals.  The common GNSS weighting strategies include: signal to noise ratio (SNR) model <cit.>, satellite elevation angle model <cit.>, SNR and satellite elevation angle hybrid model <cit.>  et al. However, computing the appropriate weights for GNSS observations in urban situations is challenging. Leveraging sky plots to divide the satellites into line-of-sight (LOS) and non-line-of-sight (NLOS), and decreasing the weights of NLOS can increase the positioning accuracy in urban settings <cit.>. Nevertheless, this approach is limited to adjusting weights at the satellite levels, which means it can not modify the weights corresponding to different frequencies from the same satellite. An additional technique for measurement quality control is fault detection and exclusion (FDE) <cit.>. A multiple-fault GNSS FDE technique is proposed by Sun et al.    <cit.>  for an integrated GNSS/IMU system. A parallel GNSS FDE technique for tightly coupled GNSS/INS/Vision integration through factor graph optimization is presented by Jiang et al.  <cit.>. However, they only perform quality control on the pseudorange observations and do not consider the carrier phase observations.

Even though significant progress has been made in the area of positioning in challenging situations, the GNSS quality control still requires extensive attention. In our earlier work, a unified cycle slip, multipath estimation, detection, and mitigation (EDM) method is proposed <cit.>, in which we utilize the clustering method to separate the cycle slips and multipath from the carrier phase observations aided by the predicted VIO positioning. However, this methodology entails repeatedly estimating the parameters using the least square method, followed by cycle slips detection and multipath modeling with the observation residuals. Whereas, the residuals are inevitably absorbed into the estimated parameters (e.g. positioning and receiver clock) with the least square method, leading to inaccurate outlier culling and cycle slip detection. Meanwhile, the method relies on the predicted VIO position to preprocess the GNSS observations. In the case of GNSS blocking frequently, setting the uncertainty of the predicted position becomes challenging. This contribution proposes an innovation-based EDM method. The innovations are applied to identify the cycle slips and multipath in the GNSS observations. To obtain the innovations, the inter-frequency bias (IFB) is estimated with the clustering method, and the receiver clock is obtained using the satellite with the highest elevation angle. Furthermore, considering the persistent impact of multipath on ambiguities, we mark the ambiguity as the cycle clip when the accumulated multipath exceeds a certain threshold to mitigate the influence of multipath effect on localization. Finally, we conduct a detailed analysis of the effectiveness of the innovation-based EDM method and residual-based EDM method. 


Our main contributions are summarized as follows: 


  * We proposed an innovation-based EDM method for the pseudorange and carrier phase observations in the tightly coupled GNSS/MEMS/Vision system.  

  * We compare and analyze the residual-based (R-EDM) and innovation-based (I-EDM) methods in terms of cycle slip detection, multipath estimation, positioning accuracy, and computational efficiency in the tightly coupled multi-sensor fusion system. 

  * Extensive road vehicular experimental evaluations are conducted in the urban area to evaluate the performance of our method. The experimental results demonstrate the superior performance of our GNSS observation preprocessing method in complex urban settings.


The introduction is followed by a description of the RTK/MEMS-IMU/Vision integration approach and our proposed innovation-based EDM method. Then, the detailed description of the vehicle-borne experimental setups and processing strategies are described in detail. The experimental outcomes of different GNSS observation weighting schemes in typical urban settings are examined, and the effectiveness of the innovation-based EDM and residual-based EDM approaches is contrasted. The conclusions are given in the end. 










§ METHODS

This section begins with an overview of the tightly coupled RTK/MEMS/Vision system. The error models of all the associated senors are next presented, followed by the time update and measurement update model in the tightly coupled filter. Lastly, we provide the specifics of our proposed innovation-based EDM algorithm. 





 §.§ System Overview

Fig.<ref> shows the system architecture of the tightly coupled RTK/MEMS/Vision system. Based on the Multi-State Constraint Kalman Filter (MSCKF), the pseudorange and carrier phase observations from GNSS, the raw data of MEMS-IMU, and the images from the stereo camera are fused. The INS mechanization is used to carry out the state propagation once the system has been initialized. The predicted state variables will also help with feature matching and tracking in the image processing, as well as the preprocessing of outlier culling and cycle slip detection in GNSS processing. The state variables and associated covariance in the estimator will be augmented upon receiving new visual or GNSS observations. Once the new observations are deemed to be available, the corresponding state variables and the covariance will be updated in the measurement update. After completing the filter, the state variables are ultimately fixed through ambiguity resolution. 




 §.§ INS Error Model

The vehicle platform's biased noisy angular velocity and linear acceleration make up the measurements of the IMU. Because of the low-cost IMU measurement noise, the Coriolis and centrifugal forces caused by Earth's rotation are disregarded in IMU's formulation. Therefore, the kinematic model of the IMU's error state can be written as  <cit.>:

    δ𝐩̇^n    = δ𝐯^n 
    δ𝐯̇^n    = -𝐑^n_b(𝐚-𝐛_a)^∧δθ-𝐑^n_bδ𝐛_a-𝐑_b^n 𝐧_a 
    δθ̇   = -(ω - 𝐛_w)^∧δθ-δ𝐛_w-𝐧_w 
    δ𝐛̇_w    = 𝐧_b_w
    δ𝐛̇_a    = 𝐧_b_a

where b and n are the IMU body frame and navigation frame (East-North-Up frame), respectively. δ𝐩̇^n, δ𝐯̇^n, δθ̇, δ𝐛̇_w  and δ𝐛̇_a represent the derivative of the position, velocity, attitude, gyroscope bias and accelerometer bias error in the navigation frame, respectively. 𝐑^n_b indicates the rotation from the body frame to the navigation frame. 𝐚 and ω are the acceleration and angular velocity measurement, respectively. The vector  𝐧_a and 𝐧_w represent the Gaussian noise of accelerometer and gyroscope measurement, while 𝐧_b_a and 𝐧_b_w are the random walk rate of the accelerometer and gyroscope measurement biases. 𝐚^∧ is the skew symmetric matrix of 𝐚. Therefore, the INS error state vector has the following writing:


    δ𝐱_ins = [δθ  δ𝐯^n   δ𝐩^n   δ𝐛_w   δ𝐛_a]^⊤




 §.§ Visual Observation Model

Taking into account that the stereo camera observes a visual point feature f^j, the related visual observation measurement 𝐳_i^j can be written as follows:


    𝐳_i^j= (
        u_i,1^j 
     v_i,1^j 
    u_i,2^j 
     v_i,2^j 
        ) = (1/z_j^C_i, 1𝐈_2 × 2   0_2 × 2
    0_2 × 2   1/z_j^C_i, 2) (
        x_j^C_i, 1
     y_j^C_i, 1
    x_j^C_i, 2
     y_j^C_i, 2) + 𝐧^j_i

where [u_i, k^j, v_i, k^j]^⊤, k ∈{1,2} indicate the feature observations of the left and right cameras on their normalized projective plane. 𝐧^j_i indicates the visual measurement noise. [x_j^C_i, n,  y_j^C_i, n, z_j^C_i, n]^⊤, n ∈{1,2} are the visual landmarks in the camera frame, which can be calculated as follows:


    (
       x_j^C_i, 1
    y_j^C_i, 1
    z_j^C_i, 1)    = (𝐑^n_C_i,1)^⊤(𝐩^n_j - 𝐩^n_C_i, 1) 
    (
       x_j^C_i, 2
    y_j^C_i, 2
    z_j^C_i, 2)    = 𝐑^C_i,2_C_i,1(𝐩^C_i,1_j - 𝐩^C_i,1_C_i, 2)

where 𝐑^C_i,2_C_i,1 and 𝐑^n_C_i,1    denote the rotation matrix that goes from the left camera frame to the right camera frame and navigation frame, respectively. The location of the left camera with regard to the navigation frame is 𝐩^n_C_i, 1. 𝐩^C_i,2_C_i, 1 is the position of the left camera frame with respect to the right camera frame. 𝐩_j^n and 𝐩_j^C_i,1 represent the visual landmarks' position in the navigation frame and left camera frame, respectively. 


In order to construct the visual reprojection residuals between relative camera poses, we use the algorithm that is proposed by <cit.>. The following is the description of the visual state vector: 


    δ𝐱_vis = [δθ^n_C_1 δp^n_C_1 δθ^n_C_2 δp^n_C_2 ⋯ δθ^n_C_k δp^n_C_k]^⊤

where δθ^n_C_i and δp^n_C_i represent the error state of the left camera's rotation and position at various time stamps, respectively. k indicates how many camera poses there are in the sliding window overall. 

The following is the expression for the projection residual of the visual measurement:


    𝐫_vis = 𝐳_vis - 𝐳̂_vis = 𝐇_vis𝐱_vis + 𝐧_vis

where 𝐳_vis and 𝐳̂_vis are the observation and reprojection visual measurements, respectively, and 𝐇_vis is the Jacobian of the relevant camera states, which are provided in <cit.>.


 §.§ Between-station Single-difference Multi-GNSS Observation Model

From the undifferenced pseudorange and carrier phase models, we will derive the between-station single-difference GNSS model.  The undifferenced pseudorange P^s_r,i and carrier phase L_r,i^s model is first given as follows:


    P^s_r,i   = ρ+c(t_r-t^s)+T_r^s+I_r,i^s+e_r,i^s,σ^2_P^s_r,i
      
         L_r,i^s   =ρ+c(t_r-t^s)+T_r^s-I_r,i^s+ λ_i · N_r,i^s+ϵ_r,i^s, σ^2_L_r,i^s

where the satellite and receiver are denoted by s and r, respectively. The carrier frequency is i=1,2,3. The speed of light is c. The pseudorange and carrier phase observations are denoted by P_r,i^s and L^s_r,i. The distance between the receiver and the satellite is  ρ. The receiver clock and satellite clock are t_r and t^s. The tropospheric and ionospheric delays on the i frequency are T^s_r and  I^s_r,i.  The wavelength of the carrier phase is  λ_i. The float ambiguity at frequency i is denoted by N^s_r,i. The noise of the pseudorange and carrier phase measurement is e^s_r,i and ϵ_r,i^s, with the variance of σ^2_P^s_r,i and σ^2_L_r,i^s, respectively.


The GNSS observations from the base station help to significantly improve the positioning accuracy of the rover station. Satellite orbit and clock biases, ionospheric and tropospheric delays in (<ref>) are eliminated by performing single-difference between the base station b and rover station r when the baseline is less than 10km.  The following is how we get the single-difference measurement model:


    Δ P^s_r,i   = Δρ+c Δ t_r + Δ e_r,i^s,Δσ^2_P^s_r,i
    Δ L_r,i^s   =Δρ+c Δ t_r +  λ_i ·Δ N_r,i^s+ Δϵ_r,i^s, Δσ^2_L_r,i^s

where Δ(·) represents the single-difference operator.

We also introduce the IFB <cit.> for the RTK model with multi-system and multi-frequency observations. Subsequently, the state vector of single-difference RTK can be expressed as follows:

    δ𝐱_rtk = [δΔ𝐩^n_r δΔt^G_r δΔ𝐈𝐅𝐁_r,i δΔ𝐍_r,i^s]^⊤
 
where δΔ𝐩^n_r indicates the error state of the baseline from a base station to a rover station. δΔt^G_r is the error state of the single-difference GPS receiver clock, i.e., the datum receiver clock.  δΔ𝐈𝐅𝐁_r,i is the error state of the single-difference IFB. And δΔ𝐍_r,i^s represents the error state of the single-difference ambiguities. 

It is worth noting that we estimate the single-difference float ambiguity in (<ref>). After the measurement update at each epoch, we will fix the between-satellite between-station ambiguities to integers, if possible, then the fixed solutions will be obtained. In this study, the integer least-square estimation of the float ambiguities is searched using the least-square ambiguity decorrelation adjustment (LAMBDA) method <cit.>. Although fixing ambiguities can enhance the positioning accuracy of RTK, if the float ambiguities are excessively erroneous, they tend to result in incorrect fixes, leading to fixed positioning solutions worse than the float solutions. Therefore, reliable initial float ambiguities are crucial for high-precision RTK positioning. 



 §.§ Tightly Coupled Filter Model for RTK/MEMS/Vision


The MSCKF is employed to estimate the real-time state variables with single-difference observations of GNSS, raw observations of MEMS IMU, and the stereo camera. To realize the tightly coupled estimation, all the corresponding state variables must be held in one estimator, which is different from loosely coupled filter <cit.>  or semi-tightly coupled filter <cit.>. The state variables in the tightly coupled filter are defined as:

    δ𝐱= [δ𝐱_ins δ𝐱_vis δΔ t_r^G  δΔ𝐈𝐅𝐁_r,i δΔ𝐍_r,i^s_δ𝐱_rtk]^⊤
 

The  Kalman filter's time updates and measurement updates are used to estimate the optimal state variables. In the time updates, the state prediction and error state covariance are propagated. The INS states are propagated forward by INS mechanization. The RTK-related state variables and the camera poses in the sliding window are regarded as constant without process noises. Regarding the error state covariance, the error state variables are propagated using a continuous system model, which is provided by: 


    [δ𝐱̇_ins
    δ𝐱̇_vis
    δ𝐱̇_rtk] = 
        [𝐅_ins   0   0
    0   0   0
    0   0   0] [δ𝐱_ins
    δ𝐱_vis
    δ𝐱_rtk] + [𝐧_ins
    0
    𝐧_rtk]
 
where the continuous-time state transition matrix of INS is represented by 𝐅_ins. The process noise of INS is 𝐧_ins, including the Gaussian noise of the accelerometer and gyroscope. And 𝐧_rtk is the process noise of RTK, containing the Gaussian noise of the receiver clock.



The discrete form of the covariance propagation based on (<ref>) can be written as: 


    𝐏_k = Φ_k, k-1𝐏_k-1Φ_k, k-1^⊤ + 𝐐_k-1

where Φ_k,k-1 is the discrete-time state transition matrix. 𝐐_k-1 is the discrete-time noise covariance matrix. The covariance matrix and state variables are augmented each time the new visual or RTK observation is captured. RTK augmentation is simple. The optimal value and covariance are carried over from the previous epoch if the state variables from that epoch are present in the current epoch, as seen in Fig. <ref>. For the newly added state variables in the current epoch, the optimal error state value is set to 0, and the variance is set as Gaussian noise. The INS mechanization would be used to initialize the camera pose for the new add image, and the augmented covariance can be written as follows:


    𝐏'_k = (𝐈_15+γ + 6k
    𝐇)  𝐏_k(𝐈_15+γ + 6k
    𝐇)^⊤

with 𝐇 = [𝐑^b_c^⊤   0   0   0   0   0_3 × (γ + 6 k)
 -𝐑_b^n (𝐩^b_c)^∧   0   𝐈   0   0   0_3 × (γ + 6 k)]. γ and k are the number of the RTK-related parameters and camera poses.  𝐑_c^b and 𝐩_c^b are the offline-calibrated extrinsic parameters that connect the IMU and camera <cit.>. 

The tightly coupled measurement update can 
be represented as follows when the GNSS or visual measurements are available: 


    [ 𝐫_vis
    𝐏^s_r,i
    𝐋^s_r,i]    =  [ 𝐳_vis - ẑ_vis
    Δ𝐏^s_r,i -  Δ𝐏̂^s_ins,i
    Δ𝐋^s_r,i -  Δ𝐋̂^s_ins,i]  
       =  [ 𝐇_vis
    𝐇_Δ𝐏^s_r,i
    𝐇_Δ𝐋^s_r,i]  [ 𝐱_vis
    δ𝐱_vis
    δΔ t^G_r
    δΔ𝐈𝐅𝐁_r,i
    δΔ𝐍_r,i^s
        ] +  [ 𝐧_vis
    Δ𝐞^s_r,i
    Δϵ^s_r,i]

where  Δ𝐋̂^s_ins,i and Δ𝐏̂^s_ins,i represent the INS-predict single-difference GNSS carrier phase measurement and the pseudorange measurement, respectively. 𝐇_vis, 𝐇_Δ𝐏^s_r,i and 𝐇_Δ𝐋^s_r,i represent the Jacobian matrix of vision representation error, pseudorange error and carrier phase error, respectively. Since the INS central position 𝐩_ins^n does not overlap with the GNSS receiver antenna reference point 𝐩_rtk^n, we also consider the lever-arm correction 𝐥^b in our study.






 §.§ Innovation-based EDM Method

Based on the hybrid model of SNR and satellite elevation angle <cit.>, we propose an innovation-based EDM method for GNSS pseudorange and carrier phase observations, which avoids residual overfitting with the least square method. As illustrated in Fig.<ref> and Fig.<ref>, the pseudorange and carrier phase data are processed independently by the innovation-based EDM approach. For pseudorange processing, we assume all the pseudorange observations share the same receiver clock Δt_r, and the observations from different frequencies and systems are modeled with IFB. There are three steps in the pseudorange processing algorithm: 



  
  1) 
 IFB initialization: at the initial epoch, IFB needs to be computed first. We choose the code observations on B1I frequency of the BDS system as the reference observations, which are the most among all the observations. The DBSCAN cluster analysis method <cit.> is utilized to obtain the stable mean value m_P_B1I. Then we clustered the observations from different frequencies after calibrating with m_P_B1I. The mean IFBs of different frequencies are computed from the most groups in the clustering results: 
    

    ΔIFB_i = 1/n∑^n_i=1(Δ P^s_r,i - m_P_B1I)


 
  
  2)
   Outlier culling: we assume the observations of satellites with the highest elevation angle are almost not affected by the multipath effect. Then the receiver clock Δt_r is computed with these observations after calibrating the IFB.  After obtaining  ΔIFB_i and Δt_r, innovations of different observations, i.e. the prefit-residuals, can be computed by:  

    ℰ_Δ P^s_r,i = Δ P^s_r,i - ΔIFB_i - Δ t_r

Finally, the weights of innovations greater than 1m are decreased from w to 1e^-3w. The reason for choosing pseudorange values greater than 1m as outliers is that the standard deviation of pseudorange is set to 0.3m in our research. Therefore, the innovation which is greater than 3 times standard deviation needs to be culled.  
  
  3)
  IFB update: we compute the m_P_B1I and update the mean IFBs of different frequencies as (<ref>) in each epoch to obtain the robust estimation of IFB. 


For carrier phase observations, the algorithm is simplified due to the adoption of between-station between-epoch double difference carrier phase observations Δ∇ L^s_r,i.
In order to determine the innovation of the observations, we further assume that all the carrier phase observations share the same receiver clock, with medium value serving as the datum. Next, using the following criteria <cit.>, the observations that belong to the good, multipath, and cycle slips are identified.




    f(ℰ_Δ∇ L^s_r,i)={
     good,       |ℰ_Δ∇ L^s_r,i| < δ_th
    
    multipath,      δ_th < |ℰ_Δ∇ L^s_r,i| <  3 ·δ_th
    
    cycle slip,       |ℰ_Δ∇ L^s_r,i| >  3 ·δ_th.

δ_th is set to 0.05m in our experiments. It is worth noting that the biased ambiguities will still have an impact on the system's positioning performance if the multipath is absorbed by the ambiguities during the Kalman filter. The multipath will be gathered and designated as the cycle slip if the aggregated value is above 0.2m to mitigate the multipath effect on ambiguity. 




§ EXPERIMENTS

Road vehicular experiments in urban regions are carried out in Wuhan, China, to assess the performance of the tightly coupled RTK/IMU/Vision system with the proposed I-EDM method. As seen in Fig. <ref>, the experimental vehicle was outfitted with two FLIR BFS-U3-31S4C-C cameras,  a low-cost ADIS-16470 MEMS IMU, a tactical-grade NovAtel SPAN-ISA-100C IMU, a time synchronization board and a Septentrio mosaic-X5 mini GNSS receiver with a NovAtel GNSS-850 antenna.  The vehicle-borne mobile system aboard the vehicle was used to collect the experimental data. Tab.<ref> contains the specification information of the consumer-grade and tactical-grade IMUs. 





Furthermore, a single  Septentrio PolaRx5 GNSS receiver was installed a mere 5km away from the rover station, which functions as the base station with exact coordinates. The single-difference GNSS measurements between the base and rover stations are produced. We used the commercial Inertial Explorer (IE) 8.9 software <cit.> to acquire the accurate smoothed solutions of the tightly coupled multi-GNSS post-processing kinematic (PPK) and tactical-grade IMU integration, with which the experimental results are verified. When
compared to MEMS-IMU, the tactical-grade IMU can sustain a particular level of pose output
accuracy for a comparatively extended period of time without the need for external corrective data, ensuring a trustworthy comparison in the assessment that follows.






The vehicular dataset, including both open-sky and urban settings, is collected between 18:04:06 and 18:41:19 on September 3, 2023, is used for a complete evaluation. Fig. <ref> displays the number of available satellites and the position dilution of precision (PDOP). The average number of LOS, NLOS, and total satellites are 14.57, 7.77, and 22.34. The average number of available satellites for GPS, BDS, Galileo, and QZSS satellite system is 4.25, 11.86, 3.64, and 2.59, respectively. The PDOP has an average value of 5.84. There is a sharp decline in the number of LOS satellites while the car is traveling between towering structures. The experimental scenarios are severely affected by GNSS NLOS, multipath effect, and cycle slip issues, negatively impacting the RTK/INS/Vision integration system's positioning performance. The trajectory's top view is shown in Fig. <ref>. Additionally, the vehicle trajectory's typical experimental situation are depicted in Fig. <ref>. The trajectory's overall length is roughly 8587m. Tab. <ref> displays the specific parameters for the GNSS, INS, and Vision modules. The processing frequency for GNSS, INS, and Vision is 1HZ, 100HZ, and 10HZ, respectively. The experiments are carried out on a PC with an Intel Core i7-9750 @ 2.6 GHz and 16 GB of RAM.



 §.§ Performance of Different Weighting Strategies

In urban canyon scenarios, the GNSS observations are susceptible to multipath and outliers, consequently decreasing the positioning performance of the GNSS/INS/Vision integration system. Therefore, choosing appropriate weights is essential for GNSS observations. We compared the positioning performance of the current mainstream weighting methods of GNSS pseudorange observations, including signal-to-noise ratio (SNR) model <cit.>, 1/sinθ^2 satellite elevation angle model <cit.>, SNR and satellite elevation angle hybrid model <cit.> and the proposed I-EDM method in RTK/INS/Vision integration system. The trajectory errors in east, north, up for SNR, satellite elevation angle, hybrid method, and innovation-based method are displayed in Fig. <ref>. The root mean square error (RMSE) of the trajectory error is computed and displayed in Tab.<ref>. We can observe that the proposed I-EDM method achieves the highest positioning accuracy of 0.31m. This is because our method effectively models the multipath effects in the GNSS signals by computing the innovation of the pesudorange observations and carrier phase observations. Therefore the accuracy of the positioning is significantly improved by eliminating the impact of the multipath effects. On the contrary, the other three methods do not account for the significant multipath effect caused by the obstructions in urban situations, resulting in the positioning accuracy only reaching the meter level. Then, by comparing the hybrid method with the SNR and satellite elevation angle method, we can see that the positioning performance of the hybrid method is better than that of the SNR and satellite elevation angle method, achieving on average a 28.1% and 38.4% improvement in 3D direction. And the performance of the SNR method is better than that of the satellite elevation angle method. The SNR method allows for assigning reasonable weights to the observations of different frequencies and types from the same satellite. However, the satellite elevation angle method only assigns the weights according to the satellite type. Consequently, the SNR method provides more accurate modelling of observations at the signal level, resulting in higher precision in positioning compared to the satellite elevation angle method. Furthermore, the hybrid method combines the advantages of both methods, achieving better positioning accuracy. 




 §.§ Comparison of R-EDM and I-EDM Method

We conduct experiments to evaluate the performance of R-EDM and I-EDM in terms of cycle slip detection, multipath estimation, positioning accuracy, and running time consumption. The aim is to provide a detailed analysis of the differences between the residual-based EDM method and the innovation-based  EDM  method in urban environments.




  §.§.§ Cycle Slip Detection Results

We compared the performance of the cycle slip detection of R-EDM and I-EDM methods. The reference trajectory generated from the tactical-grade IMU is employed to obtain the real cycle slips and multipath <cit.>. For fairness, the same criteria described in equation (<ref>) are utilized to access the cycle slips and multipath of R-EDM, I-EDM, and real values. The false detection rate of cycle slip is displayed in Fig. <ref>. The false detection rate means that the cycle slips exist within the real data, yet the method fails to detect them. This omission leads to incorrect ambiguity calculation, which is harmful to navigation performance, especially in urban experiments. The periods when cycle slip detection errors occur essentially coincide with intervals of high PDOP values, as shown in Fig.<ref>. This correlation demonstrates that multipath effects and signal loss in GNSS observations often result in cycle slips. More specifically, R-EDM and I-EDM have average error rates of 0.07% and 0.03%, respectively. This demonstrates that in the residual-based method, the cycle slips are partially absorbed into the estimated parameters such as IFB and clock biases when solving least squares, making it challenging to detect the minor cycle slips within the GNSS observations. Conversely, the innovation-based method only leverages the raw GNSS observations to detect the cycle slips, which are more sensitive to minor cycle slips.













  §.§.§ Multipath Estimation Results

We also compared the performance of the I-EDM and R-EDM methods for the multipath estimation of pseudorange and carrier phase observations. In Fig.<ref>, we plot the probability distribution histogram of the pseudorange multipath obtained by I-EDM and R-EDM methods. The multipath exceeding ±1m  detected by I-EDM is greater than R-EDM method, which is  15.70% and 11.87%, respectively. This indicates that the innovation-based method can mitigate the impact of the multipath effect in pseudorange observations on the positioning system by
modeling the multipath effect effectively. The cumulative distribution of the multipath in carrier phase observations acquired by I-EDM, R-EDM, and the actual value is displayed in Fig.<ref>. It can be observed that the multipath estimated by I-EDM is smaller than R-EDM. However, due to the between-station between-epoch double difference operator on carrier phase observations, the RMSE of the multipath estimated by I-EDM and R-EDM is similar, which is 0.011m and 0.009 m, respectively.





  §.§.§ Positioning Results

The positioning accuracy of different preprocessing methods is estimated in the tightly coupled RTK/MEMS/Vision mode. Due to the R-EDM method relying on the prior position constraints, we set the standard deviations of 0.1, 0.2, 0.3, and 0.4m to the predicted position of VIO, and analyzed the impact on the R-EDM method. Tab. <ref> displays the RMSE of the approximated trajectory errors. The position errors calculated with different uncertainty of the prior position show a maximum difference of 6cm in the vertical direction. This indicates that the performance of the R-EDM method is influenced by the uncertainty of prior position constraints. When GNSS loses tracking frequently, setting the uncertainty of the prior position estimated by IMU and vision observations becomes challenging. The I-EDM method, on the other hand, avoids the problem of setting uncertainty of prior position constraints and the deficiency estimation of the outliers. The error trajectory of the I-EDM method decreases by 21.57% compared with the R-EDM-0.1 method. Fig.<ref> is the error trajectory for epochs from 37300-37800 with severe multipath effect, which presents the trajectory errors of different methods more clearly. The benefits of the innovation-based method are further validated. 




  §.§.§ Running Time Results

We counted the time consumption for the pseudorange process and carrier phase process module in the R-EDM and I-EDM algorithms. The computation cost (in milliseconds) of different modules is shown in Table <ref>.  In the R-EDM method, employing the clustering analysis avoids searching for the outliers in the observations with the iterative method, thereby decreasing the algorithm’s time consumption. However, the algorithm still introduces the least square method to compute the residuals of pseudorange and carrier phase observations, which consumes most of the time. In the I-EDM method, we only computed the innovation of the pseudorange and carrier phase observations. Compared with the R-EDM method, the pseudorange processing, carrier phase processing, and total processing time decreased by 46.5%, 51.5%, and 49.9%, respectively.





§ CONCLUSION

The utilization of the GNSS/INS/Vision integration system provides a more comprehensive and robust positioning capability. However, multipath and cycle slips brought on by obstructions deteriorate GNSS location performance in urban areas, which further impairs multi-sensor integration positioning. For GNSS pseudorange and carrier phase measurements, this work develops an innovation-based cycle slip, multipath estimation, detection, and mitigation (I-EDM) technique. The cluster analysis method is utilized to extract the innovations from GNSS measurements, which are subsequently employed to identify cycle slips and multipath. 


The experimental results show that the proposed strategy effectively reduces the impact of outliers within observations in urban test scenarios. Our proposed method significantly enhances the positioning accuracy when compared to the signal-to-noise ratio (SNR) model, 1/sinθ^2 satellite elevation angle model, and the hybrid model.  Moreover, a detailed comparison between the performance of residual-based and innovation-based EDM methods is conducted. Compared with the residual-based EDM method, the innovation-based EDM method obtains a decrease in average error rates of cycle detection by 57.1%, a 21.6% improvement in positioning accuracy, and a 49.9% reduction in running time consumption. These results show the significant improvements achieved by the proposed innovation-based EDM method over the residual-based approach.





§ ACKNOWLEDGMENTS

Upon reasonable request to the corresponding author, the experimental data used in this research is available. 
This work is sponsored in part by The National Key Research and Development Program of China (2021YFB2501100), the National Natural Science Foundation of China (Major Program, Grant No. 42192533), and the Fellowship of China National Postdoctoral Program for Innovative
Talents (Grant No. BX20200251).




























 
 











IEEEtran











 



[
    < g r a p h i c s >
]Bo Xu received the B.E. degree in the School of Land Science and Technology from China University of Geosciences, Beijing, China, in 2014 and 2018, and the M.E. degree in the School of Geodesy and Geomatics from Wuhan University, Wuhan, China, in 2018 and 2021, respectively, where, he is currently pursuing the Ph.D degree with the School of Geodesy and Geomatics. His research interests include GNSS precise positioning, visual SLAM, visual inertial odometry (VIO) and multi-sensor fusion algorithm.


[
    < g r a p h i c s >
]Shoujian Zhang received the B.Eng. and Ph.D.
degrees in the School of Geodesy and Geomatics from Wuhan University, Wuhan, China, in 2004 and 2009, respectively. He is an associate professor at the school of geodesy and geomatics from Wuhan University. His research interests are GNSS data processing,  multiple sensor  data fusion algorithms and applications.




[
    < g r a p h i c s >
]Jingrong Wang  received the B.E. degree in the School of Land Science and Technology from China University of Geosciences, Beijing, China, in 2014 and 2018, and the M.E. degree in  the State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing (LIESMARS) from Wuhan University, Wuhan, China, in 2018 and 2020, respectively, where, he is currently pursuing the Ph.D degree with the GNSS Research Center. His research interests include GNSS precise positioning, visual inertial odometry (VIO) and multi-sensor fusion algorithm. 


[
    < g r a p h i c s >
]Jiancheng Li obtained his PhD degree from Wuhan University in 1993 and is a professor in geodesy at Wuhan University, Wuhan, China. He was honored as an academic member of the Chinese Academy of Engineering in 2011. His main research interests are the earth gravity field and its engineering applications.



