A factor set-based GNSS fault detection and exclusion for vehicle navigation in urban environments

With the rapid development of safety–critical applications of Intelligent Transportation Systems, Global Navigation Satellite System (GNSS) fault detection and exclusion (FDE) methods have made navigation systems increasingly reliable. However, in multi-fault cases of urban environments, FDE methods generally demand massive calculations and have a high risk of missed detection and false alarm. To deal with this issue, we proposed a factor set-based FDE algorithm for integrating GNSS and Inertial Measurement Units (IMU). The FDE is first performed efficiently via consistency checking over far fewer subsets of the pseudorange. Afterward, the FDE results are validated by missed-detection and false-alarm checks. The missed-detection check factor is designed by predicting the maximum horizontal GNSS positioning error, while the false-alarm check factor is designed with the aid of IMU mechanization. Following FDE, a loosely coupled GNSS/IMU integration is carried out to output the final estimation of the vehicle's position, velocity and attitude. The proposed algorithm improves both horizontal and 3D positioning accuracy by more than 50% in the field test, compared to the traditional GNSS/IMU loosely coupled scheme. Additionally, with the proposed algorithm, the accuracy improvements of the velocity and the heading are over 20% and 50%, respectively.


Introduction
For decades, the fusion of Global Navigation Satellite Systems (GNSS) and Inertial Measurement Units (IMU) has been of essential importance in its applications for vehicular navigation (Sun et al. 2010;Chen et al. 2020). The necessity for reliable navigation systems has become more stringent with the development of safety-critical applications of Intelligent Transportation Systems (ITS) (Feng and Ochieng 2007;Wang et al. 2020). However, it should be noted that GNSS pseudoranges may contain gross errors due to multipath interferences and non-line-of-sight (NLOS) receptions in urban environments (MacGougan et al. 2002;Sun et al. 2022). Faulty GNSS measurements significantly reduce the reliability of vehicular navigation systems in urban environments, with the safety of life being an issue in the worst cases. Therefore, it is of urgent necessity and importance to develop an efficient algorithm to detect and exclude faulty GNSS measurements to enhance navigation safety.
Fault detection and exclusion (FDE) methods typically work by checking the consistency of GNSS measurements (Sabatini et al. 2017;Li et al. 2020). Generally, FDE methods can be classified into two categories: snapshot FDE and recursive FDE (Zabalegui et al. 2020). Snapshot FDE only checks the consistency of current measurements, while recursive FDE utilizes both current and previously recorded historical measurements. Originally, FDE was applied in civil aviation as a major part of Receiver Autonomous Integrity Monitoring (RAIM) (Feng et al. 2006;Wang and Ober 2009). In the 1980s, classical FDE algorithms, including the pseudorange comparison method (Lee 1986), the least squares residual (LSR) method (Parkinson and Axelrad 1988) and the parity vector method (Sturza 1988), were proposed successively. These three classical FDE algorithms are all based on single-fault assumptions. With the development of multi-frequency and multi-constellation GNSS, more satellites and signals are available: This does, however, mean that the risk of multiple faults gets higher and cannot be ignored at the same time. Consequently, advanced RAIM (ARAIM) was proposed based on multiple-hypothesis solution separation (Blanch et al. 2012(Blanch et al. , 2013. In theory, ARAIM can detect multiple faults, but as many subsets are involved in the consistency check, the consumption of computational resources is very high. Notably, the FDE schemes of classic RAIM and ARAIM are all snapshots. Alternatively, recursive Kalman filter (KF)-based FDE was also developed (Bhattacharyya and Mute 2020); however, it may fail due to undetected faults in historical epochs.
To improve the reliability of vehicular navigation systems, applying FDE to land transportation, on top of civil aviation, has also been proposed. There are, however, limitations when doing so, as FDE designed specifically for civil aviation cannot be implemented directly into urban environments, where measurement redundancy is low and the possibility of simultaneous multiple faults is much higher (Zhu et al. 2018). Information space projection-based method (Kaddour et al. 2015) and recursive consistency check-based FDE (Blanch et al. 2015) have been proposed to deal with multi-fault cases. These methods, however, have a high demand for computing resources and have only been verified by simulated data. FDE based on the innovation of GNSS/IMU fusion has also been applied to improve the performance of integrated navigation systems (Hwang et al. 2005;Zhu et al. 2017;Sun et al. 2021). Though these innovation-based FDE methods are valid for multiple faults and avoid numerous subsets-based consistency checks, it should be noted that errors in IMU measurements could result in potentially false FDE. In addition, the difficulty of detecting and excluding simultaneous multiple faults also increased the risk of false alarm and missed detection greatly. However, for GNSS FDE applied in integrated navigation, most of the research directly inputs the FDE results to the filter of integration, which significantly reduces the robustness of the navigation system.
As demonstrated above, current FDE methods are defective in dealing with the performance degradation of vehicle navigation in urban areas due to these methods demanding large amounts of calculations and having high risk of false alarm and missed detection. To overcome this problem, we proposed a new factor set-based FDE algorithm. In particular, simultaneous multiple faults can be detected and iteratively excluded by consistency checking over the universal set and single-fault hypothesis subsets of the pseudorange with much fewer computational resources required. Also, the improvement of correctness and robustness of FDE is achieved by reducing the possibility of missed detection and false alarm. The missed-detection check is performed by predicting the maximum of the horizontal GNSS positioning error, while the false-alarm check is based on the vehicular position reckoned by IMU mechanization.

Algorithm framework
The framework of the proposed algorithm is shown in Fig. 1. In this figure, N B and N G are the current number of BDS and GPS pseudoranges, respectively, which could be reduced due to fault exclusion. In addition, F LS and F update are two variables used in the judgment segments of the algorithm. On the whole, the proposed FDE scheme has two steps, preliminary FDE and FDE validation. Specifically, in step one, when N B ≥ 4 and N G ≥ 4 , the fault-detection factor S LS k,0 is calculated and compared with the predetermined threshold > T LS , the fault-exclusion iteration needs to be continued until all faults have been excluded or the condition N B + N G > 6 , N B ≥ 1 and N G ≥ 1 is not satisfied.
In step two of the proposed FDE scheme, a missed-detection check or false-alarm check is performed. When F LS = 1 , the missed-detection check factor S MD k is calculated and compared with the threshold T MD . If S MD k ≤ T MD , the variances of pseudorange residuals are updated and F update is set to 1. Otherwise, F update is set to 0. When F LS = 0 , the false-alarm check factor S FA k is calculated and compared with the threshold T FA to determine the value of F update .
Above is the basic procedure of the proposed FDE. If F update = 1 , the navigation state reckoned by IMU mechanization is updated by the GNSS with the proposed FDE. The final estimation of vehicular position, velocity and attitude is then calculated by the loosely coupled GNSS/IMU fusion.

Preliminary FDE
In step one of the proposed FDE, preliminary FDE is performed and the flag F LS used in step two is determined. If N B ≥ 4 and N G ≥ 4 , the least squares positioning solution is calculated with all the pseudorange whose ionosphere delays, troposphere delays, and satellite clock errors have been corrected with corresponding models. Then, the faultdetection factor is defined as: where the subscript k represents the index of the epoch; r i is the i th element of the pseudorange residual vector r ; i is the standard deviation of the pseudorange residual of satellite i . The value of i is determined with experience initially and updated by the formula (15) in the following epochs. It is assumed that pseudoranges are only affected by nominal errors, which obey zero-mean Gaussian distributions in fault-free cases, while one or more pseudoranges contain bias in faulty cases. Based on this assumption, it can be obtained that: (2) In fault − free cases, (S LS k,0 ) 2 ∼ 2 (m, 0) where 2 (m, 0) represents central Chi-square distribution with m degrees of freedom; m equals N B +N G − 5 ; 2 (m, ) represents non-central Chi-square distribution with m degrees of freedom and non-centrality parameter . The threshold of S LS k,0 is obtained by where F −1 2 (m,0) is the inverse of the cumulative probability distribution function of 2 (m, 0)and P FA is the predetermined false-alarm probability.
(3) In faulty cases, (S LS If S LS k,0 < T LS , all pseudoranges are marked as fault-free measurements, and F LS is set to 1. Otherwise, faulty measurements are believed to exist, and fault exclusion should be performed. First, all single-fault-hypothesis subsets which exclude one pseudorange and include at least one BDS pseudorange and one GPS pseudorange are constructed: where A k,i (i = 1, 2, ⋯ , N 1 A ) represents the i th single-faulthypothesis subset; subscripts k and i mean indexes of the epoch and the subset, respectively; N 1 A is the count of subsets in the first iteration.
The test statistics corresponding to all these subsets can be obtained according to formula (1): Their thresholds can be calculated according to formula (4), but it should be pointed out that the degrees of freedom of Chi-square distribution should be reduced by one.
If the minimum test statistics, (S LS k,j ) min , is lower than the threshold T LS , the pseudorange excluded by the subset corresponding to (S LS k,j ) min is regarded as a faulty measurement, and F LS is set to 1. Otherwise, the next fault-exclusion iteration is performed on the subset corresponding to (S LS k,j ) min . And the iteration is continued until all faults have been excluded or N B + N G > 6 , N B ≥ 1 and N G ≥ 1 is not satisfied (single-fault-hypothesis subsets with redundancy cannot be constructed under this condition).

Missed-detection check
If F LS is equal to 1 after the preliminary FDE, all the remaining pseudoranges are marked as fault-free measurements by the preliminary FDE. However, it is possible that faulty measurements still exist. Hence, the theoretical maximum of undetected horizontal GNSS positioning errors is estimated to reduce the risk of missed detection. The theoretical maximum of vertical or 3D positioning errors is not used in the missed-detection check after attempts in the field test. With the predetermined possibility of missed detection, P MD , the maximum undetectable non-centrality parameter, max , can be obtained by solving: where T LS f inal and m f inal represent the threshold and degrees of freedom of Chi-square distribution corresponding to the final pseudorange set, respectively.
The observation matrix H can be obtained by are direction cosine vectors of the BDS satellite and the GPS satellite in the local-level navigation frame (north, east and down), respectively.
Let matrix H + be the pseudo inverse of the observation matrix H: and matrix can be obtained by: where I is an identity matrix whose size is the same as that of matrix HH + .
The k slope,i , which can project the pseudorange error of satellite i onto the domain of horizontal positioning error, can be calculated with the following equation (Brown and Chin 1998): where H + 1,i represents the element in the 1st row and i th column of the matrix H + ; H + 2,i represents the element in the 2nd row and i th column of the matrix H + ; S i,i is the i th diagonal element of the matrix S.
Then, the estimated maximum of the horizontal GNSS positioning error was defined as the missed-detection check factor: The threshold of S MD k is obtained by: where MD is the mean of S MD k ; MD is the standard deviation of S MD k ; is an empirical coefficient taking a value from 3 to 5.
The variances of pseudorange residuals are initialized based on experience at the first epoch. In the later epochs, they get updated by the data in the sliding window: where C B and C G are the number of satellites in BDS and GPS constellations, respectively; is the pseudorange residual of satellite i at epoch j ; L is the empirical length of the sliding window generally taking the value of 1000; null means corresponding pseudorange residual is empty because the satellite is invisible or the pseudorange is faulty.
The variance of the pseudorange residual of satellite i is updated by: where N i represents the count of non-empty pseudorange residuals of satellite i in the sliding window; is the empirical coefficient taking the value between 0.7 and 0.9; min i is the predetermined minimum of i ; pre i is the standard deviation of the pseudorange residual of satellite i in the last epoch.

False-alarm check
In the case, when F LS is equal to 0 after the preliminary FDE, either the visible satellites are too few for the proposed fault detection scheme, or the fault-exclusion iteration cannot be performed. If the GNSS data are discarded directly in this case, a measurement update would not be performed at this epoch. It increases the risk of error divergence in the integrated navigation system. Therefore, we proposed a false-alarm check scheme as follows.
If the condition N B + N G > 4 , N B ≥ 1 and N G ≥ 1 is not satisfied, F update is set to 0. Otherwise, based on the GNSS position solution and the lever arm between the GNSS antenna and the IMU, the position of the IMU can be estimated: are, respectively, estimated X-axis, Y-axis and Z-axis coordinates of the IMU in the (14) T earth-centered earth-fixed frame ( e-frame) based on the least squares position solution of GNSS.
At this epoch, the estimated position of the IMU can also be obtained by IMU mechanization: are, respectively, the estimated X-axis, Y-axis and Z-axis coordinates of the IMU in the e -frame based on IMU mechanization.
Then, the innovation of extended Kalman filter (EKF) in the position domain can be obtained as follows: The false-alarm check factor is defined as: The threshold can be obtained by: where FA is the mean of S FA k ; is an empirical coefficient that takes a value from 3 to 6; FA is the standard deviation of S FA k . If S FA k > T FA , F update is set to 0. Otherwise, the remaining GNSS measurements are remarked as fault-free measurements and F update is set to 1.

GNSS/IMU fusion scheme
As the FDE is performed in the position domain, EKF-based loosely coupled (LC) scheme is utilized for the fusion of GNSS and IMU. The state vector is defined as: where r e IMU , v e IMU and e IMU are the position error vector, the velocity error vector and the attitude error vector of IMU mechanization in the e-frame, respectively; b g and b a are the vectors of gyroscope and accelerometer three-axis biases, respectively; s g and s a are the vectors of gyroscope and accelerometer three-axis scale factors, respectively.
After the missed-detection or false-alarm check, if F update = 0 , the measurement update is not performed, and the navigation state based on IMU mechanization would be directly output.
If F update = 1 after the proposed FDE, the state vector is updated with the measurement vector: is the velocity solution in the e-frame based on GNSS Doppler shifts and the lever arm effect has been corrected; is the estimated velocity by IMU mechanization in the e-frame. Next, the estimated position, velocity and attitude are corrected by the state vector, and the final estimation of the navigation state is output.

Field test and result analysis
A field test was conducted in order to validate the proposed algorithm on November 10, 2021, in Nanjing, China. The experimental vehicle and equipment are shown in Fig. 2. The raw GNSS data were collected with a BDStar Navigation C520-AT receiver, which uses one NovAtel OEM7500 multi-frequency GNSS receiver module at a sampling rate of 10 Hz, while the raw IMU data were collected with a MEMS IMU, STIM300, at a sampling rate of 125 Hz. The reference trajectory was determined by the post-processing kinematic mode of NovAtel Inertial Explorer software with raw data of the high grade inertial/GNSS navigator, Honeywell HGuide N580. Antennas 1 and 2 are both ZYACF-S806 antennas of Zhejiang ZhongYu Communication Technology Co., Ltd.    The N580 navigator was connected with antenna 1, and the BDStar Navigation receiver was connected with antenna 2. Key IMU parameters of STIM300 and N580 are given in Tables 1 and 2, respectively. The vehicle was driven in the area near the Nanjing South Railway Station. Figure 3 shows the driving trajectory of the experimental vehicle in this test. The GNSS observation condition along the test route is exhibited in Figs. 4, 5 and 6. In Fig. 4, the number of visible GNSS (BDS and GPS) satellites is reduced greatly in a large part of the test trajectory. According to the number of satellites, three typical test cases were selected for more detailed analysis in the following. The selected test cases are open sky, GNSS challenging and GNSS denied, respectively. The environment of each test case is shown in Fig. 5. Figure 6 indicates that the position dilution of precision (PDOP) value is lower than 3 in most parts of the test route. Table 3 shows the setting of relevant parameters of the proposed algorithm. Both P FA and P MD were set to 1 × 10 −5 since we referred to the classical RAIM in civil aviation and adjusted them according to the results of the experiment. The means and standard deviations of the missed-detection check factor and the false-alarm check factor were calculated by analyzing historical data collected in the same area on November 9, 2021. Then, the thresholds T MD and T FA were set to 9 and 30, respectively.
As illustrated before, the preliminary FDE is performed firstly. Figure 7 shows the comparison of the fault-detection factor and its threshold at the start of step one of the proposed FDE. Since the scale of the vertical axis is too large for the threshold in the upper subfigure, the part in the yellow dashed rectangle is zoomed in, as seen in the lower subfigure. This figure shows that the fault-detection factor exceeds the threshold for a large part of the time, about 75% specifically. Additionally, the value of the fault-detection factor exceeds 100 in many epochs and exceeds 1000 in a few epochs, while the threshold remains at approximately 7. These characteristics of Fig. 7 are due to the complex environment of the experimental area, where faulty GNSS measurements arise easily.
After the preliminary FDE, the fault-detection factor sharply drops as illustrated in Fig. 8. It is clear that the fault-detection factor is lower than the threshold most of the time, excepting little epochs. It should be noted that the  The possibility of false alarm, P FA 1 × 10 −5 The possibility of missed detection, P MD 1 × 10 −5 The threshold of missed-detection check factor, T MD 9 The threshold of false-alarm check factor, T FA threshold is also reduced due to the exclusion of the faulty pseudoranges.
The change in the number of GNSS satellites after the preliminary FDE is depicted in Fig. 9. It is common for more than one measurement to be identified as faults by the preliminary FDE. In addition, the change of PDOP after step one is shown in Fig. 10. It is reasonable for the PDOP to get higher because some measurements are excluded. However, excluding faulty measurements generally benefits the navigation algorithm.
The preliminary fault exclusion is based on the singlefault hypothesis subsets of the pseudorange. Figure 11 shows the number of such subsets during the experiment. The maximum number of subsets at one epoch is 175, and fewer than 100 subsets are needed in most epochs.   However, far more subsets are needed with ARAIM; for example, when 10 BDS satellites and 6 GPS satellites are visible, just considering the possibility of fewer than 4 simultaneous faulty measurements, C 1 16 + C 2 16 + C 3 16 = 696 subsets should be included in calculations with ARAIM. It should be noted that the number of subsets rapidly increases when there are more visible satellites or more possible faulty measurements for ARAIM.
After the preliminary FDE, the missed-detection check or false-alarm check is performed according to the value of F update . Figure 12 depicts the missed-detection check factor and the threshold during the experiment. It is clear that the value of the missed-detection check factor is smaller than the predetermined threshold in a great number of epochs.
There are also some epochs when the missed-detection check factor exceeds the threshold. Figure 13 shows the value of false-alarm check factor and the threshold. The count of epochs when the falsealarm check was performed is less than that when the missed-detection check was performed. The false-alarm check factor is lower than the threshold in the major part of those epochs.
In order to evaluate the performance of the proposed algorithm, (1) traditional EKF-based fusion algorithm and (2) EKF-based fusion with LSR FDE are two candidate algorithms. The settings of the EKF parameters for the three algorithms are the same. The detailed steps of the algorithms are shown in Table 4.
Positioning errors of the candidate algorithms are depicted in Figs. 14 and 15. The positioning accuracy improvement is significant for the proposed algorithm, especially in the GNSS challenging area, where the maximum improvement of 3D positioning accuracy is over 40 m. The positioning root mean square error (RMSE) results and corresponding accuracy improvements compared to the traditional EKF-based fusion algorithm in the whole trajectory are shown in Table 5. It is clear that the proposed algorithm outperforms the other algorithms in the overall positioning performances. In particular, the horizontal and 3D position RMSEs of the proposed algorithm are 3.296 m and 4.562 m, respectively, corresponding to improvements of 52.2% and 56.9% over the traditional EKF-based fusion, while they are 16.6% and 13.4% for the EKF with LSR FDE.
The RMSEs and corresponding accuracy improvements over the traditional EKF-based fusion in velocity and attitude domains are shown in Table 6. Notably, more accurate velocities and attitudes are obtained with the proposed algorithm. In particular, the velocity and heading accuracy are improved by over 20% and 50% for the proposed algorithm, while the improvements are about 10% and 18.0% for the EKF with LSR FDE. Although roll and pitch accuracy improvements are 2.7% and -0.3% for the proposed algorithm, the magnitudes of roll and pitch errors are far lower than that of the heading errors in the test.
The performance comparison in the open sky area of test case 1 is shown in Tables 7 and 8. These algorithms perform similarly in the open sky area. Specifically, the 3D position RMSEs of the proposed algorithm and the EKF with LSR FDE are 1.141 m and 1.103 m, corresponding to improvements of − 0.6% and 2.7%, respectively. The improvements in velocity and attitude domains are lower than 4% for both the proposed algorithm and the EKF with LSR FDE.
As for the GNSS challenging area of test case 2, Tables 9 and 10 indicate that the proposed algorithm outperforms the other algorithms. In particular, the proposed algorithm improves both horizontal and 3D position accuracy by over 60%, while the EKF with LSR FDE improves horizontal and 3D position accuracy by 20.5% and 13.9%. In addition, the proposed algorithm also provides a heading accuracy improvement of 25.1%. However, for the EKF with LSR FDE, the heading accuracy is reduced by 0.8% .  Tables 11 and 12 show the performance comparison of these algorithms in the GNSS-denied area of test case 3. The 3D position RMSE of the proposed algorithm is 10.356 m, corresponding to an accuracy improvement of 6.6%, while they are 10.489 m and 5.4%, respectively, for the EKF with LSR FDE. Besides, the accuracy improvements Step 2: GNSS/IMU fusion is performed with EKF. The measurement vector is based on the position and velocity solutions of fault-free GNSS measurements determined by LSR FDE Proposed algorithm Step 1: Proposed FDE is implemented Step 2: GNSS/IMU fusion is performed with EKF. The measurement vector is based on the position and velocity solutions of fault-free GNSS measurements determined by the proposed FDE  of the heading for the proposed algorithm and the EKF with LSR FDE are 2.5% and 3.0%, respectively. The performance improvement of the proposed algorithm in the GNSS-denied area is much smaller than the one in the GNSS challenging area. This is because the vehicle was driven through a short tunnel during this time, causing GNSS positioning was unsolvable in most epochs of this case, and the FDE of the proposed algorithm is invalid if GNSS positioning is unavailable. Besides RMSE, other error statistics like maximum are also necessary to verify the proposed algorithm fully. For brevity of illustration, just horizontal positioning error (HPE) and heading error, which are relatively more important for the navigation of land vehicles, are further analyzed. The mean of errors, the maximum absolute value of errors (Max) and the standard deviation of errors (STD) in the horizontal position domain are shown in Table 13, respectively. As the table indicates, the Max of HPE is 2.16 m for the proposed algorithm in the open sky area, corresponding to an improvement of 13.3%. In the GNSS challenging area, the Mean of HPE is 5.02 m for the proposed algorithms corresponding to an improvement of 64.5%. As for the GNSSdenied area, the proposed algorithm provides improvements of over 6% in Mean, Max and STD of HPE. Compared with the EKF with LSR FDE, there is a much greater improvement for the proposed algorithm in the GNSS challenging area. This is because simultaneous multiple faults are excluded correctly. In addition, false alarm and missed detection are avoided at the same time. The improvement of the proposed algorithm in the open sky is smaller than that in the GNSS challenging area due to the lower possibility of simultaneous multiple faults. The improvement is also smaller in the GNSS-denied area for the proposed algorithm. This is because insufficient satellites disable the proposed algorithm at these epochs.
As Table 14 shows, in the GNSS challenging area, the Mean and Max of heading errors of the proposed algorithm are 1.67 degrees and 3.68 degrees, respectively, corresponding to improvements of over 25%. As for the GNSS-denied  area, the proposed algorithm provides an improvement of 15.1% in the STD of heading errors.

Conclusion
We have developed a novel factor set-based FDE scheme for integrated navigation of vehicles in urban environments. Simultaneous multiple faults can be detected and excluded efficiently with the proposed algorithm since far fewer subsets are included in the consistency check. Significantly, the missed-detection check factor and the false-alarm check factor are also designed to enhance the correctness and robustness of the FDE. The performance of the proposed algorithm is validated by the real-world field test. The horizontal and 3D positioning accuracy of the proposed algorithm is 3.296 m and 4.562 m, respectively, in the urban field test.
The results correspond to improvements of over 50% compared to the traditional EKF-based GNSS/IMU loose fusion algorithm. Furthermore, the proposed algorithm provides an over 20% improvement in velocity accuracy and a 52.1% improvement in heading accuracy.