## Abstract

Global navigation satellite system (GNSS) integrity monitoring (IM) has been introduced in aviation, but remains challenging for urban scenarios because of limited satellite visibility and strong multipath and non-line-of-sight effects. Consequently, factors such as limited measurement redundancy and inaccurate uncertainty modeling significantly compromise positioning and IM performance. To alleviate these issues, this paper proposes an integrity-constrained factor graph optimization model for GNSS positioning augmented by switch variables. In contrast to conventional IM methods, this method enhances redundancy through the factor graph structure. Instead of directly excluding measurements, the proposed method reweights the measurements by using switch variables to satisfy a chi-square test constraint within the optimization, ultimately yielding optimal positioning accuracy. Moreover, a proper protection level that conservatively bounds the positioning error can be derived by using the modified weighting matrix under a single-fault assumption. The effectiveness of the proposed method was verified based on data sets collected in open-sky and urban-canyon areas in Hong Kong.

- factor graph optimization
- fault detection and exclusion
- GNSS positioning
- integrity-constrained positioning
- integrity monitoring
- switch variables
- urban canyons

## 1 INTRODUCTION

Global navigation satellite systems (GNSSs) have been widely utilized in diverse applications such as aviation and ground vehicles because of their mature capability for all-day high-accuracy global positioning, navigation, and timing (Morton et al., 2021). While accuracy is the primary focus during the early phase of positioning applications, integrity concepts are vital in safety-critical applications (Zabalegui et al., 2020), notably for autonomous driving vehicles (Geiger et al., 2013). Specifically, GNSS integrity indicates the trustworthiness of the accuracy of a navigation solution (Morton et al., 2021). While integrity-monitoring (IM) techniques enable a navigation system to warn users against an unreliable GNSS status (Morton et al., 2021), incorrect IM could result in devastating consequences (Zabalegui et al., 2020). By exploiting the redundancy of available observations, IM can also be performed at the user level, which primarily includes both fault detection and exclusion (FDE) and protection level (PL) calculations (Zhu et al., 2018).

The concept of receiver autonomous IM (RAIM) (Brown, 1992; Brown, 1996; Parkinson & Axelrad, 1988; Walter & Enge, 1995) was initially investigated in aviation. In particular, RAIM provides the ability to monitor the system integrity and to timely alert users against faults based on the consistency check of an over-determined solution (Zhu et al., 2018). The main algorithms utilized in RAIM include the range-comparison-based method, the residual-based method, and the parity method based on snapshot pseudorange measurements (Brown, 1992). Whereas the majority of RAIM methods focus on detecting and isolating a single fault in the Global Positioning System (GPS) (Walter & Enge, 1995), recent techniques consider multiple faults due to applications of multiple constellations and more visible satellites (Hewitson & Wang, 2006; Lee, 2004; Wang & Xu, 2020). The advanced RAIM (ARAIM) technique (Blanch et al., 2012) extends the monitoring ability to multiple-fault assumptions and multi-constellation GNSS observables and enables an explicit computation of the integrity risk allocation based on the multiple-hypothesis solution-separation (MHSS) method (Pervan et al., 1998). Compared with conventional RAIM, which utilizes a single test statistic for detection, ARAIM establishes multiple test statistics in the position domain for each fault hypothesis (Joerger et al., 2014). Except for RAIM/ARAIM approaches that consider snapshot pseudorange measurements, sequential IM applies a Kalman filter (KF) or particle filter (PF) as a recursive estimator to make use of the measurement sequence and the system motion model to obtain optimal accuracy and integrity in the current state estimates (Gupta & Gao, 2019; Tanil et al., 2018; Wang et al., 2023). In general, KF-based sequential IM can be separated into techniques that utilize an MHSS-based detector, which maintains filter banks for the fault hypotheses (Wang et al., 2023), or an innovation-based detector with a single KF (Tanil et al., 2018). PF-based IM estimates integrity using direct measurement probabilistic modeling via particles (Gabela et al., 2021; Gupta & Gao, 2019; Gupta & Gao, 2021; Zhang & Gui, 2015). Pushing sequential IM beyond KF methods, one recent work (Hafez et al., 2020) also extended IM for fixed-lag smoothing to quantify robot localization safety by considering all measurements within a sliding window; both solution-separation and chi-squared IM methods were developed.

The above-mentioned IM techniques are suitable for open-sky environments with satellite redundancy and reliably modeled noise distributions of nominal measurements. However, these approaches encounter significant challenges in degraded urban scenarios because of unique circumstances such as limited satellite availability and multipath and non-line-of-sight (NLOS) effects, leading to difficulty in obtaining a reliable uncertainty model of nominal measurements (Zhu et al., 2018). Although previous works used robust estimation to ease the effect of faulty satellites via formulation of an equivalent weighting matrix (Crespillo et al., 2020; Crespillo et al., 2018; Wang et al., 2018), such implementations require tackling the parameters of specific robust cost functions and an additional inertial navigation system (INS). In addition, conventional GNSS IM directly excludes faulty GNSS measurements until the chi-square test is passed under a pre-defined probability of risk, which can significantly degrade the satellite geometry in urban canyons. Moreover, both RAIM and ARAIM tend to provide a loose bound of integrity risk, which is suitable in aviation but could be too conservative for vehicles in urban scenarios (Pullen et al., 2011; Reid et al., 2019). Recent advances in IM for ground vehicles in urban canyons also extend its application to multi-modal sensors (e.g., INS, radar, lidar, camera) (Meng & Hsu, 2020; Wang et al., 2023) to improve safety redundancy. However, these approaches still suffer from the high cost of lidar, high computation loads, and the lack of GNSS measurement redundancy brought by frequent and multiple multipath and NLOS receptions (Zhu et al., 2018).

Recent studies have utilized factor graph optimization (FGO) for GNSS positioning (Bai et al., 2022; Das et al., 2021; Watson & Gross, 2017; Wen & Hsu, 2021; Wen, Pfeifer, et al., 2021; Wen, Zhang, & Hsu, 2021), which improves measurement redundancy and obtains better positioning accuracy compared with filtering-based sequential IM by considering historical information and applying multiple iterations and linearization (Wen, Pfeifer, et al., 2021). By further applying FDE for historical measurements, Wen, Meng, and Hsu (2021) conducted a chi-square test over multiple epochs and showed that more redundant healthy measurements result in better FDE performance. However, performance improvement is limited when most of the measurements are faulty or when the given uncertainty modeling is unreliable. In contrast, conventional RAIM/ARAIM techniques operate separately alongside the positioning component, meaning that the positioning results obtained from all measurements provide an initial guess for IM and then faulty measurements are excluded via FDE. After these steps, PL and new positioning results are obtained correspondingly. This type of separation is not robust in two aspects: 1) When healthy measurements are excluded, they cannot be reconsidered in the next iteration of IM. 2) When the corresponding measurement variance mismatches with the actual variance, FDE may fail, and the calculated PL cannot adequately bound the position error. Although prior works on particle RAIM (Gabela et al., 2021; Gupta & Gao, 2019) have performed joint state estimation and IM using particle reweighting, they are hindered by the following issues: 1) high computation complexity due to numerous particles and resampling steps, 2) suboptimality of the PF owing to analytical intractability, and 3) neglect of previous measurement faults in consequence of the first-order Markov property from the PF-based scheme. Thus, the following question arises: *Can we automatically select/exclude faults and estimate measurement uncertainties simultaneously in a tractable mechanism with extended redundancy to alleviate the challenges listed above?* In fact, a recently developed switchable constraint (SC) algorithm (Sünderhauf & Protzel, 2012a) shed some light on this problem.

The SC algorithm, which was recently recognized as a robust estimation method (Yang et al., 2020), was initially introduced in simultaneous localization and mapping to detect and remove outlier constraints from factor graph problems (Sünderhauf & Protzel, 2012a) and was later exploited in GNSS positioning to identify and remove multipath- or NLOS-affected measurements in urban canyons, without any additional a priori knowledge or additional sensor information (Sünderhauf et al., 2013). The SC is implemented by associating a switch variable with each pseudorange measurement factor that could be an outlier. Specifically, the switch variable is designed to continuously range from 0 to 1, where 0 indicates exclusion of the corresponding measurement, 1 denotes measurement selection as an inlier, and a value between 0 and 1 accounts for measurement variance modification (Sünderhauf & Protzel, 2012b; Suzuki, 2021). A closed solution for switch variables has been obtained from an analysis of the partial derivatives for switch variables when the optimizer converges (Agarwal et al., 2013), indicating that switch variables can iteratively reweight measurements without changing the geometrical relationship between the state estimation and measurements. Therefore, it is feasible to apply the switch variable to incorporate FDE into GNSS positioning by FGO: *The test statistic is represented as the sum of normalized residuals of switch-variable-associated pseudorange factors within the factor graph, and its threshold depends on the number of measurements and a given false alarm probability*. Hence, the inequality constraint arising from the chi-square test in FDE can be formulated by ensuring that the chi-square test statistic is equal to or less than the threshold. In such a way, the additional switch variables can be iteratively estimated to satisfy the chi-square test inequality constraint, meaning that all involved measurements are iteratively identified during optimization, avoiding the exclusion of healthy measurements. The PL can be either calculated or estimated accordingly, following the idea of RAIM/ARAIM while considering the reweighting effect of switch variables.

In this study, switch variables will be applied to a conventional factor graph-based GNSS positioning framework as an attempt to demonstrate its potential to improve RAIM performance. In this paper, only pseudorange measurements are considered within the batch FGO. First, FGO will be implemented with switch variables, together with switch prior factors that aim to avoid local minimum. Then, the above-mentioned FDE inequality constraint will be analytically formulated. Next, the PL will be calculated via a geometry matrix. The entire FGO framework will be evaluated in terms of the performance of the FDE constraint and PL compared with conventional FGO-based positioning without switch variables (Wen, Meng, & Hsu, 2021).

The contributions of this paper are summarized as follows:

This paper proposes to incorporate FDE into FGO through switch variables to improve measurement redundancy by dynamically reweighting measurements through SCs. Instead of directly excluding the potentially faulty measurements, the reweighting inferred from switch variables helps to maintain the satellite geometry, consequently leading to improved positioning accuracy.

This paper presents an approach for PL calculation that considers the reweighting effect of switch variables to conservatively bound the position error, which is further validated by using data sets collected in both open and urbanized areas of Hong Kong.

The remainder of this paper is organized as follows: Section 2 includes an overview of the proposed method. The GNSS positioning method using FGO with switch variables is presented in Section 3. A performance evaluation of the proposed method based on experiments conducted on real data sets for Hong Kong is provided in Section 4. Finally, conclusions and future work are summarized in Section 5.

## 2 OVERVIEW OF THE PROPOSED METHOD

A flowchart of the proposed method is depicted in Figure 1. The system inputs comprise raw pseudorange measurements received from the GNSS receiver. These measurements, along with satellite position, clock bias, and a weighting matrix **W** based on the signal-to-noise ratio (SNR) and elevation angle (Herrera et al., 2016), are fetched from the raw GNSS data. Subsequently, the pseudorange measurements are associated with switch variables to construct switchable pseudorange factors. Moreover, based on the concept of the chi-square test in FDE, the chi-square test factor is formulated with the total number of pseudorange measurements *M*, the switchable pseudorange residual , and the given false alarm rate *P _{FA}*. Consequently, the chi-square test factor constrains the batch test statistic ψ

_{n}to not surpass the batch test threshold

*T*. In addition, the incorporation of switch variables motivates two types of factors: switch prior factors and switch reliable factors. Later, after the factor graph composed of the four aforementioned types of factors has been solved, the state estimation and the updated geometry matrix after reweighting are obtained. The behavior of FDE through the chi-square test factor is illustrated in Figure 1(c). In contrast to the skyplot for FDE, which directly excludes the faulty satellites (Figure 1(b)), the proposed method assigns diminished switch values (weightings) to faulty satellites, as denoted by the smaller radius of the circles, thereby preserving the geometry. Afterward, the modified weighting matrix is derived from the switch variable estimations and the weighting matrix

_{D}**W**. Finally, the horizontal PL (HPL) is derived based on the modified weighting matrix

**W**

_{s}, the updated geometry matrix

**G**, and the snapshot test threshold . This calculation is specifically conducted based on the single-fault assumption (Walter & Enge, 1995). The geometry matrix

**G**and take into consideration snapshot switchable pseudorange factors. Note that HPL calculations employ the snapshot test threshold , whereas the chi-square test factor utilizes the batch threshold

*T*for batch FDE.

_{D}## 3 METHODOLOGY

This section presents the methodology of integrity-constrained FGO for GNSS positioning. This paper exclusively utilizes GPS and BeiDou pseudorange measurements. The states of the GNSS receiver are represented as follows:

1

where **χ** denotes the set of states of the GNSS receiver *r* from the first epoch to the current epoch *n*. The state of the GNSS receiver *r* at a single epoch *t* can be expressed as follows:

2

where *t* ∈ [1, n] and the state **x**_{r,t} includes the state variable *x*_{r,t}, which comprises the Earth-centered Earth-fixed coordinates **p**_{r,t}(*p _{r,t,x}*,

*p*,

_{r,t,y}*p*) and the receiver clock biases of GPS and BeiDou, along with the switch variable associated with satellite

_{r,t,z}*j*. Each switch variable specifically ranges from 0 to 1, where 0 denotes measurement exclusion, 1 indicates measurement selection, and a value between 0 and 1 represents measurement reweighting. The initial switch value is set to 1. The superscript

*j*denotes the satellite index, and

*m*is the total number of received satellites at epoch

*t*.

The graph structure of the proposed framework is illustrated in Figure 2. Each blue circle denotes a satellite, with only two satellites drawn at each epoch for simplicity. The green rectangles represent the switchable pseudorange factors (e.g., ), and the purple rectangles denote the switch prior factors (e.g., ). The light blue rectangles indicate the switch reliable factors (e.g., ), and the red rectangle signifies the chi-square test factor (e.g., e^{cst}). The white circles represent the state of the GNSS receiver.

If the same satellite *j* appears in consecutive states, the switch variable associated with satellite *j* is connected via the switch reliable factor. The chi-square test factor is linked to all switch variables in the FGO. Specifically, our proposed method utilizes the chi-square test factor to conduct a batch FDE by considering all pseudorange measurements across epochs. Assuming that at most one fault with large weighting remains in each epoch after all possible faulty measurements have been reweighted via switch variables, the HPL calculation is simply based on the single-fault assumption (Walter & Enge, 1995) to indicate the reweighting effect of switch variables on the HPL.

### 3.1 Switchable Pseudorange Factor

Each raw pseudorange measurement from satellite *j* is obtained by the GNSS receiver and is denoted as follows:

3

where is the geometric range between satellite *j* and the GNSS receiver *r* at epoch *t*. *c* denotes the speed of light. For simplicity, δ_{r,t} represents the receiver clock bias from a single satellite constellation (GPS or BeiDou), and indicates the clock bias of satellite *j*. and represent the ionospheric and tropospheric delay errors, respectively, which are modeled according to the method in RTKLIB (Takasu & Yasuda, 2009). represents the errors introduced by multipath and NLOS receptions, receiver noise, and antenna phase-related noise.

The observation model for pseudorange measurement of satellite *j* is given as follows:

4

where **p**_{r,t} and denote the positions of the receiver and satellite at epoch *t*, respectively. The term represents Gaussian white noise with . Here, denotes the covariance of satellite *j*, which is calculated based on the SNR and satellite elevation angle (Herrera et al., 2016). Therefore, the error function of the switchable pseudorange factor for a given measurement is derived as follows:

5

where represents the switch variable associated with satellite *j* and denotes the corresponding error function of the pseudorange factor.

### 3.2 Switch Prior Factor

Initially, all satellites are accepted as healthy satellites. The switch prior factor anchors switch variables at an initial value of 1, preventing the optimization process from becoming stuck in a local minimum with all zero switch values (Sünderhauf et al., 2013; Suzuki, 2021). The corresponding error function is given as follows:

6

where is the covariance of the switch prior factor. Note that the choice of significantly affects performance; a smaller value indicates that the corresponding satellite is more likely to be healthy, leading the switch variable to be estimated closer to 1. The setting of can vary depending on whether the scenario is more or less urbanized. For simplicity, this value is set as 1 for all satellites, and its role is further discussed in the experimental validation for different settings.

### 3.3 Switch Reliable Factor

Because pseudorange measurements from the same satellite are time-correlated among epochs, switch variables from the same satellite in consecutive epochs are likely to be equal (Sünderhauf et al., 2012). The switch reliable factor acts as a constraint to prevent the switch variable from rapidly changing, leading to the following error function formulation:

7

where is the covariance of the switch reliable factor.

### 3.4 Chi-Square Test Factor

The chi-square test in classic RAIM (Walter & Enge, 1995) considers snapshot pseudorange measurements; the chi-square test factor, however, involves all pseudorange measurements in the graph. The detailed formulation is as follows:

First, the *m* switchable pseudorange residuals for epoch *t* are denoted as follows:

8

Then, the normalized residuals can be formulated as follows:

9

where is the standard deviation of pseudorange measurement *j*. Considering *n* epochs in the graph, the normalized residuals among *n* epochs are represented as follows:

10

Therefore, the test statistic ψ_{n} is formulated as the square root of the weighted sum of squared errors over all switchable pseudorange residuals among epochs:

11

Specifically, assuming *M* pesudorange measurements in the graph and a pre-given false alarm rate *P _{FA}*, the test threshold

*T*can be calculated through the following equation:

_{D}12

where *f _{χ}*(

*s*;

*v*) denotes the chi-squared distribution with the degree of freedom (DOF)

*v*=

*M*– 5

*n*. Then, a comparison with the threshold

*T*(

_{D}*v*,

*P*) is conducted. If ψ

_{FA}_{n}≤

*T*, then the position estimation is assumed to be safe under the false alarm rate

_{D}*P*; otherwise, a fault is declared. Therefore, the error function e

_{FA}^{cst}of the chi-square test factor can be formulated as follows:

13

where *α* is a step-size parameter for modifying the weighting and Σ^{cst} denotes the covariance of the corresponding residual.

### 3.5 GNSS Positioning via FGO

Based on the factors derived above, the objective function for GNSS positioning via FGO is formulated as follows:

14

The variable **χ*** denotes the optimal state set estimation, obtained by solving the above objective function. In this paper, Ceres (Agarwal & Mierle, 2012) is utilized as the nonlinear optimization solver.

The complete Jacobian matrix of the optimization problem is formed by stacking the partial Jacobians of the four types of factors: , , , and *e ^{cst}*, as shown in Figure 3 (left). Except for the Jacobian of the chi-square test factor, the Jacobians are sparse, given that most of the matrix entries are zero. Compared with the Jacobian of the single-point positioning (SPP) algorithm shown in Figure 3 (right), correlations among switch variables across epochs are facilitated by the switch reliable factors and chi-square test factors.

To prevent the optimization from rejecting all measurements with zero switch values, the following initialization strategy is leveraged when the new epoch *t* +1 comes:

The previous estimation results are set as initial guesses for the previous states

**x**_{r,1},…,**x**_{r,t}.If satellite

*j*is observed at both epochs*t*and*t*+1, then serves as the initial guess for the state ; otherwise, 1 is assigned as the initial guess.If the satellite number at epoch

*t*+1 exceeds 4, then the corresponding SPP result is used as the initial guess for**p**_{r,t+1}and**δ**_{r,t+1}; otherwise, the initial guess is provided by the last estimation results and .

### 3.6 PL Calculation

After the optimization problem in Equation (14) has been solved, the estimated switch variable is used to construct the modified weighting matrix **W**_{s} as follows (Sünderhauf & Protzel, 2012b):

15

where the weighting matrix **W** denotes the inverse of the covariance matrix of *m* obtained pseudorange measurements at epoch *t*. Next, we let denote the vector containing *m* pseudorange residuals in a single epoch and represent the five-dimensional offset vector (east, north, up, GPS clock, and BeiDou clock) near a fixed linearization point from the last iteration. The linearized pseudorange observation model in a single epoch is formulated as follows:

16

where **G** is the geometry matrix with dimensions *m* × 5 and **ε** is the ranging error with the modified covariance matrix . Then, the weighted least-square estimate of can be written as follows:

17

where:

18

The estimate of the ranging errors is given by the following:

19

From the error estimates, we can calculate the snapshot test statistic ψ_{r,t} for a single epoch (Angus, 2006; Walter & Enge, 1995):

20

with:

21

If the batch version of the chi-square test in Section 3.4 passes, indicating that ψ_{n} ≤ *T _{D}*(

*v*,

*P*), then the following inequality also holds:

_{FA}22

where is the test statistic threshold calculated from Equation (12) with the snapshot pseudorange measurement number *m* at epoch *t* and the DOF *v** = *m* −5. After passing the chi-square test described in Section 3.4, the HPL that protects for both bias and noise under the single-fault assumption is derived as follows:

23

Here, *k _{H}* denotes the number of standard deviations, which is empirically set to 3 based on the 3

*σ*rule.

*H slope*indicates the relationship between the horizontal position error and the pseudorange bias of the faulty satellite, and

_{i}*H slope*denotes the maximum of these parameters.

_{max}*H slope*can be calculated as follows:

_{i}24

Note that the HPL calculation is based on the single-fault assumption (Walter & Enge, 1995), except that the weighting matrix **W**_{s} considers the reweighting effect of the switch variables. While there exists a method considering multiple faults (Angus, 2006) that may be more appropriate for urban settings, we assume only a single fault to quantify the impact of the switch-variable-modified weighting matrix **W**_{s} on the HPL.

## 4 EXPERIMENTAL RESULTS AND DISCUSSION

To evaluate the performance of the proposed method, experiments were conducted with real pseudorange observables in an open-sky scenario and our recently released open-source UrbanNav data sets (Hsu et al., 2021). The receiver is static in the open-sky scenario, and *P _{FA}* is set to 10

^{−4}for both scenarios. The performance assessment focuses on the east, north, and up frames. Because GNSS positioning targets ground vehicles, the evaluation specifically considers horizontal positioning and integrity performance.

### 4.1 Open-Sky Experiment

#### 4.1.1 Experimental Setup

During a static experiment with the open-sky data set, raw GPS/BeiDou measurements were collected by a u-blox M8T GNSS receiver at a frequency of 1 Hz, observing a total of 16 satellites. Gaussian random errors were manually injected into observables from satellites of pseudorandom noise (PRN) 2, 6, and 15 with a mean *μ* = 15 m and a standard deviation *σ* = 5m to introduce faulty measurements. For simplicity, measurements with injected biases are considered faulty, and potential faults in other satellites are ignored. The ground truth location represented in latitude, longitude, and altitude is [114.17°, 22.31°, 63 m]. We aim to use the open-sky data set with a controllable bias to evaluate the effectiveness of the proposed method in reweighting outlier measurements and calculating the HPL. The experiments were conducted on a desktop computer with an Intel Core i9-12900K central processing unit. To evaluate the contribution of this paper, we compare the following methods:

**FGO**: Batch integration of pseudorange factors via FGO (Sünderhauf et al., 2012).**RAIM-FGO**: Batch integration of pseudorange factors via FGO. Faulty measurements are manually excluded at each epoch to simulate snapshot residual-based FDE (Blanch et al., 2015).**SW-FGO**: Batch integration of switchable pseudorange factors , switch prior factors , and switch reliable factors via FGO (Sünderhauf et al., 2012).**SWFDE-FGO (proposed method)**: All four types of factors illustrated in Figure 2 are used within the FGO for GNSS positioning. In contrast to RAIM-FGO, the proposed method conducts batch FDE with the assistance of the chi-square test factor*e*.^{cst}

#### 4.1.2 Evaluation of Positioning Accuracy and HPL

The mean horizontal positioning error (HPE) and HPL of GNSS positioning obtained via the four above-mentioned methods are shown in Table 1. A mean HPE of 3.23 m is achieved using only pseudorange factors in FGO, with a minimum mean HPL of 11.70 m. By applying residual-based FDE (Blanch et al., 2015) in FGO, RAIM-FGO manually excludes faulty measurements, leading to a lower mean HPE of 2.72 m, with a slightly increased mean HPL of 13.11 m. Interestingly, both SW-FGO and SWFDE-FGO achieve better positioning accuracy but a larger HPL than the first two methods, indicating that switch variables can reweight the measurements and influence the positioning performance. Additionally, the first two methods both maintain the satellite geometry instead of directly excluding outlier measurements. Specifically, SWFDE-FGO obtains the lowest mean HPE of 1.31 m and a less conservative HPL than SW-FGO. This finding reveals that the proposed method could lead to more conservative HPL values compared with the conventional methods (FGO and RAIM-FGO) because of the increased GNSS measurement uncertainties arising from the switchable variables. The bound rate of the HPL is defined as the percentage of epochs for which HPL ≥ HPE. GNSS IM is satisfied if the calculated HPL can bound the actual HPE. The bound rates are 100% for all four methods, indicating that the applied weighting method (Herrera et al., 2016) for this open-sky scenario is conservative.

Table 1 shows that SWFDE-FGO achieves a less conservative HPL than SW-FGO, primarily because of the differing switch value estimation of these two methods.

Figure 4 shows the snapshot test statistic ψ_{r,t} throughout the experiment, with the initial test threshold denoted by a black dashed curve, where *v** = 16 −5 and *P _{FA}* = 10

^{−4}. FGO has the largest statistic because of faulty measurements. RAIM-FGO (cyan curve) stays close above because it excludes the faulty measurements but keeps the original conservative weighting matrix

**W**. Because of potential faults beyond PRN 2, 6, and 15, the statistic of RAIM-FGO slightly exceeds the threshold. It is anticipated that the fault-free statistic would remain below the threshold if the experiment were conducted using simulated open-sky data. Both SW-FGO and SWFDE-FGO achieve statistics lower than , as they reweight all measurements within the geometry. The chi-square test factor helps keep the statistic closer to , resulting in relatively larger switch variable estimations.

Figure 5 illustrates the relationship between the absolute normalized pseudorange residuals and the switch variables . It is shown that for a large residual, SWFDE-FGO estimates a larger switch value than SW-FGO. Here, we let **W**_{s,SWFDE}, HPL_{SWFDE} and **W**_{s,SW}, HPL_{SW} denote the modified weighting matrix and HPL of SWFDE-FGO and SW-FGO, respectively, with the following:

25

Based on Equations (23) and (24), we can conclude that HPL_{SWFDE} ≤ HPL_{SW}. This phenomenon primarily arises from the aggressive behavior of SW-FGO in assigning switch values to measurements with relatively larger residuals, although these measurements may still contribute to positioning accuracy. By utilizing reweighting-based FDE via the chi-square test factor, SWFDE-FGO can provide more appropriate switch values that promise a lower HPE and less conservative HPL.

#### 4.1.3 Evaluation of the Chi-Square Test Factor

Figure 6 illustrates the HPEs (solid curves) and HPLs (dashed curves) during the experiment. The shaded gray area denotes the bias injection period. Throughout this injection period, RAIM-FGO (cyan), SW-FGO (green), and SWFDE-FGO (blue) consistently exhibit a lower HPE than FGO (magenta). This trend indicates that both switch variables alone and the combination of switch variables with a chi-square test factor act as an FDE mechanism to mitigate the impact of faulty measurements.

To further demonstrate the effectiveness of FDE with the proposed method, skyplots of the four methods at epoch A are presented in Figure 7. Note that the circle’s radius represents the switch value, which is assigned as 1 for FGO and RAIM-FGO. It can be clearly observed that RAIM-FGO conducts FDE by simply excluding the faulty measurements (PRN 2, 6, and 15). In contrast, both SW-FGO and SWFDE-FGO conduct FDE by assigning small switch values to the faulty measurements, thereby maintaining the geometry and obtaining optimal weightings during positioning. Note that even a healthy satellite can obtain a small switch value (e.g., PRN 52), primarily because of its poor geometry. Additionally, SWFDE-FGO tends to assign larger switch values than SW-FGO (e.g., PRN 15, 19, and 62) as larger values keep the statistic closer to its threshold, thereby reducing the error of the chi-square test factor *e ^{cst}*.

### 4.2 Urban Scenario Experiment

In this subsection, we evaluate the performance of the proposed method on the challenging UrbanNav data set (Hsu et al., 2021). Raw pseudorange measurements from GPS and BeiDou were collected at a frequency of 1 Hz by a u-blox M8T GNSS receiver, while ground truth information was provided by the NovAtel SPAN-CPT system. We chose to utilize the UrbanNav-Medium and UrbanNav-Harsh parts of the data set. The UrbanNav-Medium data set consists of environmental elements such as wide streets, medium-height buildings, and one-sided buildings. In contrast, the UrbanNav-Hash data set features narrow streets, bridges, and tall buildings. Specifically, we calculate the mean mask elevation angle *μ _{MEA}* induced by building blockages to indicate their urbanization rates (Wen et al., 2020) as follows:

26

where *θ _{a}* represents the mask elevation angle at the azimuth angle

*a*and

*N*indicates the number of azimuth angle increments in the sky mask. Here,

*N*= 360 because the resolution of the azimuth angle is 1°. For instance,

*μ*= 0° indicates a purely open-sky scenario, whereas

_{MEA}*μ*= 90° depicts an indoor scenario. Details of the selected data sets in terms of environmental elements and urbanization rates are given in Table 2.

_{MEA}The experiments were implemented on a desktop computer equipped with an Intel Core i9-12900K central processing unit. As there is no prior knowledge of faulty measurements, we skip the evaluation of RAIM-FGO and compare only the following five methods:

**SPP:**Single-point positioning (Enge, 1994).**RAIM-SPP:**SPP assisted with the snapshot residual-based FDE (Blanch et al., 2015).**FGO**: Batch integration of pseudorange factors via FGO (Sünderhauf et al., 2012).**SW-FGO**: Batch integration of switchable pseudorange factors , switch prior factors , and switch reliable factors via FGO (Sünderhauf et al., 2012).**SWFDE-FGO (proposed method)**: All four types of factors illustrated in Figure 2 are used within the FGO for GNSS positioning. In contrast to RAIM-SPP, the proposed method conducts batch FDE with the help of the chi-square test factor.

Figure 8 shows the trajectory obtained by the five methods and the ground truth for both data sets. The black curve denotes the ground truth trajectory. Because the UrbanNav-Medium data set (Figure 8(a)) has a lower urbanization rate, SWFDE-FGO achieves a smoother trajectory than the other methods overall, with a significant performance improvement at approximately epoch 95655 s. In contrast, the performances of all methods are severely degraded on the UrbanNav-Harsh data set (Figure 8(b)) because of the more challenging scenario and fewer healthy satellites. Detailed evaluations are presented in the following two subsections.

#### 4.2.1 Performance Evaluation for UrbanNav-Medium

The positioning results evaluated for the UrbanNav-Medium data set are shown in Table 3. FGO and SPP both attain comparable mean HPEs of 11.80 and 11.90 m. By incorporating switch variables, SW-FGO and SWFDE-FGO reduce the mean HPE to 10.35 and 9.01 m, respectively. The proportion of errors smaller than 10 m increased to 72.73% and 77.69%, respectively, from 55.37%. Notably, RAIM-SPP achieves a mean HPE of 9.72 m and a proportion of 72.31%. This result demonstrates that incorporating the chi-square test factor helps FGO in effectively performing FDE alongside conventional residual-based RAIM and benefits GNSS positioning in light urban canyons via redundant healthy measurements.

It is crucial to ensure that the HPL accurately bounds the HPE for ground vehicle positioning integrity. The mean HPLs for SPP, RAIM-SPP, and FGO are 44.47, 49.37, and 49.23 m, respectively. While these are relatively small values, they fail to bound the HPE all the time (bound rate: ~97%). In contrast, both SW-FGO and SWFDE-FGO show higher mean HPLs with 100% bound rates due to the reweighting effect of switch variables. The proposed method (SWFDE-FGO) achieves a mean HPL of 74.45 m, which is less conservative compared with SW-FGO (165.61 m), and more tightly bounds the HPE, verifying the effectiveness of the chi-square test factor.

Figure 9 depicts the snapshot test statistic (a), HPE (b) and HPL (c) of the five methods in the experiment. The black dashed line in Figure 9(a) represents the snapshot statistic threshold . The statistics of SPP and FGO overlap during the experiment, while those of SPP and RAIM-SPP coincide when the statistic stays below the threshold because FDE is not activated. The statistics for FGO and SPP frequently surpass the threshold because of the presence of faulty measurements in urban scenarios. In contrast, SW-FGO and SWFDE-FGO maintain a lower statistic below the threshold compared with RAIM-SPP, as the switch variables downweigh most of the possible faulty measurements. Notably, SWFDE-FGO approaches closer to the threshold to minimize the cost of the chi-square test factor, ensuring optimal positioning accuracy at most epochs and a less conservative HPL overall. The evaluations above reveal the effectiveness of the proposed method in maintaining the geometry during FDE while providing a less conservative HPL to appropriately bound the HPE.

#### 4.2.2 Performance Evaluation for UrbanNav-Harsh

This subsection analyzes the effectiveness of our method on the more challenging UrbanNav-Harsh data set. Table 4 illustrates the evaluation metrics throughout the campaign. Owing to the lack of healthy measurements, both RAIM-SPP and SWFDE-FGO exhibit limited improvement in accuracy, with mean HPEs of 15.07 and 15.96 m, respectively. The percentage of HPEs smaller than 10 m rises to 47.62% for RAIM-SPP and 47.77% for SWFDE-FGO, yet no substantial correction is evident for epochs in which HPEs exceed 20 m. Compared with the other methods, SW-FGO and SWFDE-FGO yield larger mean HPLs of 166.04 and 95.71 m, respectively. As switch variables modify the original measurement weightings, the corresponding bound rate is consequently higher than that of others. The presence of the chi-square test factor prevents the switch variables from aggressively down-weighing the measurements, resulting in less conservative HPLs.

Figure 10 depicts the total satellite number (blue) and NLOS satellite number (red) labeled by the three-dimensional (3D) building model during the experiment in the top panel, the HPE of the five methods in the middle panel, and the HPL in the bottom panel. The effectiveness of FDE for both RAIM-SPP and SWFDE-FGO is closely related to the proportion of NLOS satellites among the total satellites. During periods in which the majority of satellites are healthy (e.g., shaded areas A and C with 2–3 NLOS satellites), both methods can mitigate the impact of faulty measurements and improve positioning accuracy. In contrast, during periods in which most satellites are contaminated (e.g., shaded areas B and D), both methods fail to correct the faulty measurements. SW-FGO and SWFDE-FGO exhibit larger HPLs than the other methods. SWFDE-FGO helps to ease the increase in HPL with the help of the chi-square test factor.

To further explore the behavior of the chi-square test factor faced with different numbers of NLOS satellites, skyplots of SPP, RAIM-SPP, SW-FGO, and SWFDE-FGO are presented in Figure 11 for epoch 185174 s in area A and epoch 185241 s in area B. In these plots, green circles indicate line-of-sight satellites, while red circles represent NLOS satellites, with the corresponding PRN numbers shown inside the circles. The circle’s radius represents the switch value, and the shaded area represents the sky mask due to building blockages.

In Figure 11(a), with three NLOS satellites, RAIM-SPP and SWFDE-FGO show similar behavior by excluding or de-weighting the most inconsistent satellite PRN 30, whereas SW-FGO de-weights both satellites PRN 30 and 17. Compared with SPP, all three methods obtain optimal positioning accuracy with redundant healthy satellites. Conversely, in Figure 11(b), where the majority of satellites are contaminated, RAIM-SPP excludes the faulty satellites PRN 17, 59, and 60, but fails to exclude all faulty measurements. SW-FGO de-weights additional faulty satellites PRN 8 and 39 and healthy satellite PRN 7, and SWFDE-FGO de-weights additional faulty satellite PRN 21. SW-FGO and SWFDE-FGO evaluate the switch value of each satellite to maintain the most consistent but erroneous group of satellites, resulting in a degraded positioning result.

The snapshot test statistic is shown in Figure 12. Similar to the results obtained for experiments conducted in open sky and UrbanNav-Medium, SW-FGO and SWFDE-FGO reduce the statistic below the threshold, with the statistic of SWFDE-FGO approaching closer to the threshold. However, as illustrated in Figure 10 and Figure 11, statistics below the threshold cannot ensure optimal positioning because of NLOS/multipath effects and a lack of measurement redundancy. To showcase the correction effect of our proposed method, the switchable pseudorange errors of SW-FGO and SWFDE-FGO are shown in Figure 13. Interestingly, SW-FGO and SWFDE-FGO confine the (switchable) pseudorange error to a narrower range than FGO. While SW-FGO restricts most pseudorange errors near 0 m, SWFDE-FGO exhibits a greater distribution of large pseudorange errors to keep the statistic closer to the threshold. This result indicates that satellites with large pseudorange errors still provide geometric benefits for positioning.

#### 4.2.3 Discussion of the Chi-Square Test Factor

In addition to the formulation of the chi-square test factor given in Equation (13), an alternative implementation involves inequality constraints through the hinge loss function (Xie et al., 2020), as depicted below:

27

Intuitively, we aim to penalize the objective function only if the statistic ψ_{n} exceeds the threshold *T _{D}*. Let SWFDEH-FGO denote the method utilizing Equation (27). Figure 14 illustrates the result of SW-FGO (green), SWFDE-FGO (blue), and SWFDEH-FGO (yellow) in terms of the statistic (top panel), HPE (middle panel), and HPL (bottom panel). It is evident that the curves of SWFDEH-FGO and SW-FGO almost overlap, indicating that the chi-square test factor, as formulated in Equation (27), introduces no cost for optimization, causing SWFDEH-FGO to degenerate into SW-FGO. The main reason for this behavior is that SW-FGO behaves so aggressively that it can ensure ψ

_{n}<

*T*, making

_{D}*e*always zero. For the chi-square test to work in switch-variable-assisted FGO, it is necessary to ensure that the statistic does not deviate too far from the threshold, even when the statistic is already below the threshold. This precisely represents the effect of SWFDE-FGO.

^{cst}#### 4.2.4 Discussion of Switch Prior Factors

As mentioned in Section 3.1, various configurations of can significantly impact the performance of the proposed method. This subsection examines these effects by setting as 0.5^{2}, 1^{2}, or 1.5^{2} in SWFDE-FGO.

Figure 15 illustrates the snapshot test statistic, HPE, and HPL of SWFDE-FGO for various values of . A smaller covariance (0.5^{2}) in switch prior factors improves the confidence in switch values close to 1, thereby aligning the statistic more closely with the threshold. This setting allows for a further reduction in HPL due to larger switch values. Conversely, when a weak constraint is applied with , the method encounters local minima with all zero switch value estimates (e.g., shaded areas A and B). Consequently, the corresponding statistics rapidly decrease to zero, causing the HPL and HPE to surge to nonsensically large values.

To further elucidate these results, we introduce a large horizontal alert limit of 350 m to assess the availability of the three settings. Table 5 provides a comprehensive evaluation encompassing availability along with the HPE, HPL, and bound rate when the system is available. Notably, when the covariance increases from 0.5^{2} to 1.5^{2}, the availability decreases from 96.37% to 93.36%, as a larger covariance suffers from the aforementioned local minimum issue. The setting of 0.5^{2} attains optimal positioning performance, with a mean HPE of 12.67 m and a mean HPL of 58.46 m. This result suggests that a small covariance proves advantageous for positioning under challenging scenarios with limited redundancy and frequently occurring faults.

## 5 CONCLUSIONS AND FUTURE WORK

In this study, a novel framework of RAIM for factor-graph-based GNSS positioning, inclusive of FDE and PL calculation, has been proposed. First, the chi-square test is integrated into the optimization process, facilitated by switch variables. Subsequently, the HPL is computed by considering the reweighting effect of switch variables following conventional residual-based RAIM.

The proposed method was evaluated on open-sky and urban data sets with a comparative analysis against the snapshot residual-based RAIM. The results demonstrate the effectiveness of the proposed method in executing reweighting-based FDE via switch variables. This approach avoids the exclusion of satellites, thereby preserving the satellite geometry. Furthermore, the calculated HPL adopts a conservative value, effectively bounding positioning errors, as switch variables modify the pre-defined measurement variance.

It is noteworthy that in more challenging scenarios, the performance improvement of the proposed method is constrained because of limited satellite visibility and a high number of faulty measurements. More importantly, the performance of our method relies heavily on the switch prior factors. In future work, we will investigate the tuning strategy of the switch prior factors. In addition, we plan to utilize a differential GNSS technique to mitigate system errors. We will also consider incorporating PL estimation into the factor graph while considering constraints from integrity risk requirements in ARAIM. This incorporation aims to provide a positioning solution and tight PL in the presence of changeable measurement uncertainties. Furthermore, we will investigate additional sensors, such as inertial measurement units, to establish correlation among states, which will contribute to FDE by leveraging historical measurement redundancy within the factor graph.

## HOW TO CITE THIS ARTICLE

Xia, X., Wen, W., & Hsu, L.-T. (2024). Integrity-Constrained factor graph optimization for GNSS positioning in urban canyons. *NAVIGATION, 71*(3). https://doi.org/10.33012/navi.660

## ACKNOWLEDGMENTS

This work was supported in part by the Guangdong Basic and Applied Basic Research Foundation under Grant No. 2021A1515110771, the University Grants Committee of Hong Kong through the Research Impact Fund under Grant No. R5009-21, and the Faculty of Engineering of the Hong Kong Polytechnic University through the Project Perception-based GNSS PPP-RTK/LVINS integrated navigation system for unmanned autonomous systems operating in urban canyons.

Xiao Xia and Weisong Wen are with the Hong Kong Polytechnic University Shenzhen Research Institute, China, and Hong Kong Polytechnic University, Hong Kong (e-mail: xxiao.xia{at}connect.polyu.hk; welson.wen{at}polyu.edu.hk).

Li-Ta Hsu is with Hong Kong Polytechnic University, Hong Kong (e-mail: lt.hsu{at}polyu.edu.hk).

The authors would like to thank Yihan Zhong and Liyuan Zhang for their support with data set preparation and Hoi-Fung Ng for processing the ground truth data.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.