Performance Analyses of a RAIM Algorithm for Kalman Filter with GPS and NavIC Constellations

Abstract

This paper evaluates the performance of an integrity monitoring algorithm of global navigation satellite systems (GNSS) for the Kalman filter (KF), termed KF receiver autonomous integrity monitoring (RAIM). The algorithm checks measurement inconsistencies in the range domain and requires Schmidt KF (SKF) as the navigation processor. First, realistic carrier-smoothed pseudorange measurement error models of GNSS are integrated into KF RAIM, overcoming an important limitation of prior work. More precisely, the error covariance matrix for fault detection is modified to capture the temporal variations of individual errors with different time constants. Uncertainties of the model parameters are also taken into account. Performance of the modified KF RAIM is then analyzed with the simulated signals of the global positioning system and navigation with Indian constellation for different phases of aircraft flight. Weighted least squares (WLS) RAIM used for comparison purposes is shown to have lower protection levels. This work, however, is important because KF-based integrity monitors are required to ensure the reliability of advanced navigation methods, such as multi-sensor integration and vector receivers. A key finding of the performance analyses is as follows. Innovation-based tests with an extended KF navigation processor confuse slow ramp faults with residual measurement errors that the filter estimates, leading to missed detection. RAIM with SKF, on the other hand, can successfully detect such faults. Thus, it offers a promising solution to developing KF integrity monitoring algorithms in the range domain. The modified KF RAIM completes processing in time on a low-end computer. Some salient features are also studied to gain insights into its working principles.

1. Introduction

Since its inception, satellite navigation has become a critical infrastructure that supports a significant part of everyday life today. Its diverse array of applications includes power grid synchronization, precision farming, and emergency services, to name a few [1]. Although the system keeps the modern world running, it is vulnerable to numerous sources of signal degradation, anomalies, and intentional or unintentional interference. It is estimated that a failure of the global navigation satellite systems (GNSS)—a collective term for different satellite navigation systems—could cost modern economies billions of dollars [2]. Furthermore, growing reliance on GNSS for safety- and liability-critical applications, such as civil aviation and autonomous navigation, requires adherence to stringent performance standards for continuity and reliability of services. Faults or anomalies without timely detection can lead to hazardous consequences in such operations. Potential failure modes of GNSS for precise positioning are discussed in detail in [3].

While the performance of the global positioning system (GPS) has improved since 2012 with a lower number of fault counts [4], continuous monitoring is required by the constellation service provider as well as at the user level for high-reliability applications. Newer constellations (e.g., navigation with Indian constellation (NavIC) of India [5]) would possibly need several years to reach a state with few observed faults, thus requiring constant monitoring. Integrity is considered an important performance metric for operations that demand high reliability. It quantifies the level of trust placed in GNSS positioning. Receiver autonomous integrity monitoring (RAIM) is a GNSS receiver-based algorithm implemented in the user equipment to monitor integrity. It carries out consistency checks to detect faults and outputs position error bounds or protection levels (PL) when no anomaly is found. If an anomalous behavior is detected, the system needs to issue an alert to the user within a pre-determined time called the time to alert to avoid misleading information.

RAIM can be broadly classified into two categories depending on the standard implementations of the navigation processor of the receiver—weighted least squares (WLS) RAIM and Kalman filter (KF) RAIM. There exists a substantial amount of research on WLS RAIM [6,7,8]. In contrast, KF RAIM is much less explored [9]. KF is considered an integral component of advanced navigation algorithms. Examples are precise point positioning with GNSS [10], vector receiver architecture [11], and multi-sensor integration [12]. KF is a preferred choice in these systems for its ability to combine dynamics or a priori information with measurements and to estimate residual sensor errors [13]. The aforementioned advanced methods are shown to enhance performance by either improving accuracy, tracking high dynamics, or bridging GNSS outages in challenging environments. In order to ensure reliability along with these desirable features, integrity monitoring is an essential step. As a result, developing KF RAIM algorithms has received significant attention in recent years. A literature survey of contemporary work on KF-based integrity monitoring is provided next.

Reference [9] reviews different sensor integration architectures for intelligent transportation systems and state-of-the-art methods for fault detection and exclusion (FDE) and PL calculations. In this context, all KF-based methods up to the time of writing are also discussed. In [14], an FDE method is developed for the tight coupling of GNSS and inertial navigation system (INS) to detect multiple GNSS as well as GNSS and INS faults. It uses two detectors to separate GNSS and INS faults and a two-step fault exclusion process. A no subset advanced RAIM (ARAIM) algorithm is derived in [15] for the tight integration of GNSS and INS under some assumptions. It is shown to reduce the computational load to a great extent and provide global coverage for CAT I precision approaches. Reference [16] designs a computationally efficient KF RAIM algorithm for GNSS receivers. It checks measurement inconsistencies in the range domain, models time correlated errors, and requires Schmidt KF (SKF) as the navigation processor. A KF integrity monitoring algorithm is proposed in [17] for robot localization. It monitors faults within a time window with innovation-based test statistics and is also robust against previous faults. In [10] a navigation system of GNSS, INS, and an odometer is developed, and tight PLs are produced with ARAIM for precise point positioning. A parallel processing strategy coupled with the efficient software design of [18] shows that real-time operations of a substantially large number of KF banks are feasible on a GPU for integrity monitoring.

In [19], an FDE scheme is presented with the colored noise of KF represented as a first order autoregressive model. Fault detection tests for a real-time kinematic application show improved FDE performance and lower false alarm rates. An adaptive interactive multiple model framework is described in [12] to switch between Huber KF and RAIM for robustness and reliability. The proposed model offers better positioning performance than individual sub-models for a relative navigation system with GNSS and ultra wideband observations. Reference [20] discusses a Kalman integrated PL approach with empirically tuned parameters. It demonstrates better performance of the new approach over two other methods for an autonomous vehicle state estimation with GNSS, INS, and an odometer. A sliding window-based FDE method is proposed in [21] for integrated GNSS–INS. It is shown to improve positioning accuracy over that without fault exclusion. RAIM algorithms for advanced nonlinear filters [22] and vector receivers [23] are also described in the literature.

Five of the aforementioned references discuss KF-based PL calculations. Among them, [10,15] perform fault detection in the position domain, [16,17] follow range domain methods, and [20] assumes that fault detection and mitigation have been executed. This paper attempts to explore range domain methods further for their relatively low architectural complexity. To this end, it extends the work in [16], which uses an SKF [24] navigation processor and carrier-smoothed pseudoranges for absolute positioning. SKF accounts for measurement errors as ”consider parameters” without estimating them. As explained in the prior work, it allows the formation of two simultaneous fault detection tests, each with terms of a finite number of epochs. This ensures that the KF PLs do not grow with time. However, the methodology assumes a 100-second time constant for all errors, which is unrealistic, and hence, an important limitation of the work. In addition, the algorithm is validated with simulated range and range rate measurements, not generated by processing GNSS signals with a receiver. Therefore, the measurements do not properly represent the carrier smoothing of code. Fault detection performance is also not analyzed.

In order to address the above shortcomings, first, realistic error models of the carrier-smoothed pseudorange measurements are integrated into KF RAIM. More precisely, the error covariance matrix for fault detection is modified to capture the temporal variations of individual errors with different time constants. Uncertainties of the model parameters are also taken into account. Performance of the modified KF RAIM is then evaluated with the simulated signals of GPS and NavIC for aircraft navigation during different phases of flight. Intermediate frequency (IF) data are post-processed by a developed GPS L1-NavIC L5 software receiver to generate carrier-smoothed code and precise range rate measurements for RAIM. Some salient features of KF RAIM are also studied to gain insights into its working principles.

Performance analyses show that WLS RAIM [25] used for comparison purposes has lower PLs than the KF RAIM. However, WLS is not suitable for advanced navigation methods, such as multi-sensor integration and vector receivers, where KF is an integral part. Therefore, KF-based integrity monitors are essential. An important finding of the analyses is as follows. Innovation-based tests with an extended KF (EKF) navigation processor mistake slow ramp for the residual measurement error that the filter estimates, leading to missed detection. RAIM with SKF, on the other hand, can successfully detect such faults. It also completes processing in time on a low-end computer. Thus, it offers a promising solution to developing KF RAIM in the range domain. The current work does not consider constellation-wide faults. Measurements from different satellites at a particular time epoch are assumed to be uncorrelated with one another.

The remaining paper is organized as follows. First, the KF RAIM extended in this paper is briefly outlined for ease of understanding of the subsequent modifications. Next, extensions of the algorithm are presented in detail. Following this, simulation studies are performed. Finally, the paper concludes with a summary of the key findings and future work.

2. Overview of Existing KF RAIM

Both fault detection test statistics and PLs of KF RAIM presented in [16] are briefly discussed in this section for ease of understanding of the subsequent developments.

2.1. Fault Detection Test Statistics

Assuming only white noise in the GNSS measurements, three fault detection test statistics (${\alpha }_{j,k}$; j = $1,\dots ,3$) at $t_k$ are formulated in [26] as

$${\alpha }_{j,k}^2=\sum _{\ell =p_{j,k}}^{m_{j,k}}\left(\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }^T{\mathbf{W}}_{\ell }^{-1}(\mathbf{I}-{\mathbf{G}}_{\ell })\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }+\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }^T{{\mathbf{W}}_{rr,\ell }}^{-1}(\mathbf{I}-{\mathbf{G}}_{\ell }^{rr})\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }\right)$$

(1)

where $\mathbf{I}$ represents the identity matrix of dimension $n\times n$. n is the number of visible satellites at a given epoch and can change with time. $\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }$ and $\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }$ are pseudorange and pseudorange rate innovation vectors of the EKF, respectively, at time epoch ℓ. The limits of ℓ on the right hand side, $m_{j,k}$ and $p_{j,k}$, change with j and k. ${\mathbf{G}}_{\ell }$ and ${\mathbf{G}}_{\ell }^{rr}$ are given by WLS estimation as

$${\mathbf{G}}_{\ell }={\mathbf{H}}_{\ell }{\left({\mathbf{H}}_{\ell }^T{\mathbf{W}}_{\ell }^{-1}{\mathbf{H}}_{\ell }\right)}^{-1}{\mathbf{H}}_{\ell }^T{\mathbf{W}}_{\ell }^{-1}$$

(2)

$${\mathbf{G}}_{\ell }^{rr}={\mathbf{H}}_{\ell }{\left({\mathbf{H}}_{\ell }^T{{\mathbf{W}}_{rr,\ell }}^{-1}{\mathbf{H}}_{\ell }\right)}^{-1}{\mathbf{H}}_{\ell }^T{{\mathbf{W}}_{rr,\ell }}^{-1}$$

(3)

where ${\mathbf{H}}_{\ell }$ is the linearized measurement model matrix for pseudoranges (or pseudorange rates) and of dimension $n\times 4$ for a single constellation and $n\times 5$ for dual constellations. ${\mathbf{W}}_{\ell }$ is the error covariance matrix of the pseudorange measurements at $t_{\ell }$. ${\mathbf{W}}_{rr,\ell }$ is its counterpart for pseudorange rates.

WLS estimation-based fault detection in KF is used because it helps in the following ways. It results in the noise statistics of $\alpha$s being easily computed. Furthermore, each term under the summation of Equation (1) contributes the fault of its own epoch only, not of any other epochs. Thus, the computation of the failure mode slope (FMS) for mean position error bound is simplified. In addition, each $\alpha$ can have a separate FMS. Consequently, three tests can be formulated. With three suitably designed test statistics, the PL is shown not to increase with time in [26], as opposed to the PL with a single statistic.

The limitation of the above formulation is that measurement errors across time epochs are assumed to be uncorrelated, which is a gross simplification. To overcome this limitation, [16] assumes the following error models. Broadcast satellite orbit and clock error, residual tropospheric error, receiver noise, and multipath are together represented as a first order Gauss Markov (FOGM) random process with a 100 s time constant (carrier smoothing interval of code). From [27], individual error standard deviations are obtained for GPS carrier smoothed code for airborne users. The second constellation considered is NavIC. As for NavIC, the same equations are used with frequency bands properly changed. This is because equivalent error models are yet to be developed for this constellation. No ionospheric error is modeled, which can be eliminated with dual frequency measurements. Only white noise pertaining to 45 dB-Hz C/N${}_0$ is assumed to be present in the pseudorange rate measurements.

Considering time-correlated errors in the carrier-smoothed pseudorange measurements, and dropping time instant k for simplicity of notation, the first part of Equation (1) for an epoch changes to

$$\begin{array}{c}\hfill {\left({\alpha }_j^{\prime }\right)}^2={\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}}^T{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}\mathsf{\Theta }{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}\end{array}$$

(4)

where $\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}$ stacks all $\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }^{*}$ vectors given by $\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }^{*}={\mathbf{W}}_{\ell }^{-1/2}\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }$, ℓ varying from $p_j$ to $m_j$; $\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }$ denotes the carrier-smoothed pseudorange innovation vector at epoch ℓ; and ${\mathbf{W}}_{\ell }$ is its error covariance. ${\tilde{\mathbf{S}}}^{\prime }$ is a full rank block matrix, with its ℓth block along the diagonal given by ${\mathbf{S}}_{\ell }^{\prime }$, where $(\mathbf{I}-{\mathbf{G}}_{\ell }^{*})={\mathbf{S}}_{\ell }^{\prime }{{\mathbf{S}}_{\ell }^{\prime }}^T$ and ${\mathbf{G}}_{\ell }^{*}={\mathbf{W}}_{\ell }^{-1/2}{\mathbf{H}}_{\ell }{\left({\mathbf{H}}_{\ell }^T{\mathbf{W}}_{\ell }^{-1}{\mathbf{H}}_{\ell }\right)}^{-1}\times {\mathbf{H}}_{\ell }^T{\mathbf{W}}_{\ell }^{-1/2}$, and other elements are zero. $\mathsf{\Theta }$ is a square matrix with dimension $\left({\sum }_{\ell =p_j}^{m_j}n\right)\times \left({\sum }_{\ell =p_j}^{m_j}n\right)$. ’×’ represents the multiplication operation. ${\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}$ on either side of $\mathsf{\Theta }$ in Equation (4) ensures that only the ℓth epoch’s fault is included from ${\mathsf{\Delta }{\underline{\mathit{\rho }}}^{*}}_{\ell }$.

In order to perform fault detection tests with ${\left({\alpha }_j^{\prime }\right)}^2$, SKF implementation [28] is used in the navigation processor. SKF offers reasonably good performance, while being computationally efficient. It models time-correlated errors without estimating additional states. As a result, these errors are not changed in the pseudorange innovations. Thus, SKF along with the help of ${\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}$ allows a separate FMS for each of the test statistics. Hence, multiple test statistics can be formed. Each is with a given number of epochs.

Equation (4) is finally reduced to the form given by

$${\left({\alpha }_j^{\prime }\right)}^2={{\underline{\mathit{\eta }}}_r}^T{\mathsf{\Theta }}_2{\underline{\mathit{\eta }}}_r$$

(5)

where ${\underline{\mathit{\eta }}}_r$ consists of unit-variance independent elements. Matrix ${\mathsf{\Theta }}_2$ is determined by ${\tilde{\mathbf{S}}}^{\prime }$, $\mathsf{\Psi }$ (measurement error covariance matrix of $\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}$), and $\mathsf{\Theta }$. Since the maximum eigenvalue of ${\mathsf{\Theta }}_2$, ${\lambda }_{\mathsf{\Theta },2}^{\max }$ is designed to be less than or equal to one, fault detection thresholds of ${\left({\alpha }_j^{\prime }\right)}^2$ can be derived from the chi-squared distribution [29]. $\mathsf{\Theta }$ is formulated such that ${\lambda }_{\mathsf{\Theta },2}^{\max }$ is less than or equal to one. It should be noted that ${\lambda }_{\mathsf{\Theta },2}^{\max }$ is also equal to the maximum eigenvalue of ${\tilde{\mathbf{S}}}^{\prime T}\mathsf{\Psi }{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime \mathbf{T}}\mathsf{\Theta }{\tilde{\mathbf{S}}}^{\prime }$ (=$\mathsf{\Omega })$, ${\lambda }_{\mathsf{\Omega }}^{\max }$.

Thus, the test statistics shown in Equation (1) are changed with time-correlated errors in the carrier-smoothed pseudorange measurements as

$${\alpha }_j^2={\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}}^T{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}\mathsf{\Theta }{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime T}\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}+\sum _{\ell =p_j}^{m_j}\left(\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }^T{{\mathbf{W}}_{rr,\ell }}^{-1}(\mathbf{I}-{\mathbf{G}}_{\ell }^{rr})\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }\right)$$

(6)

where j = 1, 2. At an epoch, if ${\alpha }_1$ crosses its threshold $T_{th,1}$, or ${\alpha }_2$ is greater than its threshold $T_{th,2}$, or both occur, a fault is detected. The third test is not formulated. Necessary modifications for not performing the third test are made to PLs. ${\alpha }_1$ is formed with terms of M recent epochs including the current one, whereas ${\alpha }_2$ has terms of $(N-M)$ previous epochs. Both M and N are calculated adaptively at each time epoch. The matrices $\mathsf{\Psi }$ and $\mathsf{\Theta }$ are briefly described next.

$\mathsf{\Psi }$ can be re-arranged into a block diagonal matrix, with ith block for satellite i given as ${\mathsf{\Psi }}_{\mathrm{sub}}^i=\left[\begin{array}{ccccc}1& a_{11}& a_{11}{a}_{21}& a_{11}{a}_{21}{a}_{31}& \cdots \\ a_{11}& 1& a_{21}& a_{21}{a}_{31}& \cdots \\ a_{11}{a}_{21}& a_{21}& 1& a_{31}& \cdots \\ & & ⋮& \end{array}\right]$, where $a_{\ell 1}$ at an epoch $t_{\ell }$ takes a value between $a_{\min }$ and $a_{\max }$; $a_{\min }$ = 0.9890; and $a_{\max }$ = 0.9911 with a 5° elevation mask angle. In this context, another matrix ${\mathsf{\Psi }}_{\max }$ is also defined. It has the same structure as that of $\mathsf{\Psi }$, but its sub-matrix ${\mathsf{\Psi }}_{\mathrm{sub},\max }^i$ is formed with $a_{\ell 1}$ always set as $a_{\max }$. $\mathsf{\Theta }$ is given by ${(c\mathbf{I}+\mathsf{\Psi })}^{-1}$ for ${\alpha }_1^{\prime }$, with c being a positive multiplication factor calculated adaptively to ensure that ${\lambda }_{\mathsf{\Theta },2}^{\max }$ is less than or equal to one. For ${\alpha }_2^{\prime }$, $\mathsf{\Theta }$ is ${(c\mathbf{I}+{\mathsf{\Psi }}_{\max })}^{-1}$, with c as $({\lambda }_{\mathsf{\Psi },\max }^{\max }-{\lambda }_{\mathsf{\Psi },\max }^{\min })$. ${\lambda }_{\mathsf{\Psi },\max }^{\max }$ and ${\lambda }_{\mathsf{\Psi },\max }^{\min }$ are the maximum and minimum eigenvalues of ${\mathsf{\Psi }}_{\max }$.

2.2. Protection Levels (PLs)

The SKF estimates error in the a priori user position and velocity in the Earth Centered Earth Fixed (ECEF) frame, and the receiver clock offsets and drifts for GPS and NavIC from respective truths. Furthermore, the time correlated error of each carrier-smoothed pseudorange measurement is modeled as an FOGM process but not estimated by the SKF. First, mean position error bounds are discussed, and then PL equations are provided.

2.2.1. Mean Position Error Bounds

The mean position error, ${\underline{\mathit{\mu }}}_k^{\prime }$, of the navigation filter at epoch $t_k$ in the local navigation frame can be written as

$${\underline{\mathit{\mu }}}_k^{\prime }=-\underset{\mathrm{Ignored}\mathrm{terms}}{\underbrace{\sum _{j=N-1}^{k-(m+1)}{\mathbf{B}}_{k-j}^{\prime }{\underline{\mathit{f}}}_{k-m-j}}}-\underset{N\mathrm{terms}\mathrm{corresponding}\mathrm{to}{\alpha }_1\&{\alpha }_2}{\underbrace{\left(\sum _{j=0}^{N-2}{\mathbf{B}}_{k-j}^{\prime }{\underline{\mathit{f}}}_{k-m-j}+{\mathbf{K}}_{s,k}^{\prime }{\underline{\mathit{f}}}_{k-m+1}\right)}}$$

(7)

where ${\mathbf{K}}_{s,k}^{\prime }$ = ${\mathbf{C}}_e^n{\mathbf{K}}_{s,k}(1:2:5,:)$ and ${\mathbf{B}}_{k-j}^{\prime }$ = ${\mathbf{C}}_e^n{\mathbf{B}}_{k-j}(1:2:5,:)$. ${\mathbf{B}}_{k-j}$ = $\left[{\prod }_{\ell =k}^{k-j}{\mathbf{A}}_{\ell }\right]{\mathbf{K}}_{s,k-j-1}$, and ${\mathbf{A}}_{\ell }=(\mathbf{I}-{\mathbf{K}}_{s,\ell }{\mathbf{C}}_{s,\ell }){\mathbf{F}}_{s,\ell -1}$. ${\mathbf{F}}_s$, ${\mathbf{C}}_s$, and ${\mathbf{K}}_s$ are the state transition, measurement model, and Kalman gain matrices of SKF for the estimated states. For a matrix $\mathbf{J}$, $\mathbf{J}(1:2:5,:)$ comprises the first, third, and fifth rows of $\mathbf{J}$. In ${\mathbf{B}}_{k-j}$ and ${\mathbf{K}}_{s,k}$, these rows belong to the position states of the filter. ${\mathbf{C}}_e^n$ is the coordinate transformation matrix from ECEF to the north-east-down (NED) frame. The two test statistics, ${\alpha }_1$ and ${\alpha }_2$, together have terms of N recent epochs. Terms from other previous epochs are ignored. An upper bound of the contributions from these terms to the mean position error is derived mathematically. Fault is represented by a vector ${\underline{\mathit{f}}}_{k-m+1}$ at $t_k$; it started at epoch m. Three or more simultaneous faults are not required to be monitored for the integrity risk considered. If $s_1$ and $s_2$ represent any two faulty satellite indices, then ${\underline{\mathit{f}}}_{k-m+1}\left(s_1\right)=b_{1,k-m+1}$, ${\underline{\mathit{f}}}_{k-m+1}(s_1+n)={\dot{b}}_{1,k-m+1}$, ${\underline{\mathit{f}}}_{k-m+1}\left(s_2\right)=b_{2,k-m+1}$, and ${\underline{\mathit{f}}}_{k-m+1}(s_2+n)={\dot{b}}_{2,k-m+1}$, with all other elements being zero. For single fault (SF) cases, only the $s_{1}th$ satellite has a fault. The mean error bound of the vertical position for ${\alpha }_1$ under a fault mode i (which can be SF or dual faults (DF), with i = $1,\dots ,N_{fault}$), $VP{E}_{1,i}^U$ is obtained as follows.

For a fault mode, the square of the maximum FMS [25] for ${\alpha }_1$ (having M recent epoch terms from $t_k$ to $t_{k-M+1}$) is

$$(\mathrm{maximum}{\mathrm{FMS})}_{\alpha 1,i}^2=\underset{{\underline{\mathit{f}}}_w^T\mathsf{\Lambda }{\underline{\mathit{f}}}_w=\mathsf{\Gamma }}{max}{\displaystyle \frac{{\underline{\mathit{f}}}_w^T\mathbf{L}{\underline{\mathit{f}}}_w}{{\underline{\mathit{f}}}_w^T\mathsf{\Lambda }{\underline{\mathit{f}}}_w}}$$

(8)

The FMS is maximized with the worst-case fault vector, ${\underline{\mathit{f}}}_w$. $\mathbf{L}$ is formed with the appropriate elements of ${\mathbf{B}}^{\prime }$ matrices over M recent epochs for a fault mode. ${\underline{\mathit{f}}}_w$ constitutes the corresponding b and $\dot{b}$ terms. $\mathsf{\Lambda }$ is

$$\begin{array}{cc}\hfill \mathsf{\Lambda }=& \left[\begin{array}{cc}{\mathsf{\Theta }}_3& \left[\mathbf{0}\right]\\ \left[\mathbf{0}\right]& {\mathsf{\Theta }}_4\end{array}\right]\hfill \end{array}$$

(9)

where $\left[\mathbf{0}\right]$ denotes a matrix of all zero entries. ${\mathsf{\Theta }}_3$ is determined with the terms of Equation (4) corresponding to a fault mode. ${\mathsf{\Theta }}_4$ is obtained from the second term of Equation (6) on the right. It is a diagonal matrix for SF modes and a block diagonal matrix with each block of size 2 × 2 for DF modes. A computationally efficient method is used to calculate the inverse of each block of ${\mathsf{\Theta }}_4$ without inverting it. If ${\lambda }_{\mathbf{L}{\mathsf{\Lambda }}^{-1}}^{\max }$ is the maximum eigenvalue of $\mathbf{L}{\mathsf{\Lambda }}^{-1}$,

$${\mathrm{VPE}}_{1,i}^U=\sqrt{{\lambda }_{\mathbf{L}{\mathsf{\Lambda }}^{-1}}^{\max }\mathsf{\Gamma }}$$

(10)

where $\mathsf{\Gamma }$ is calculated as follows:

$$\begin{array}{cc}\hfill {\left(T_{th,1}^{*}\right)}^2& ={\displaystyle \frac{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}{{\left({\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }\right)}^2}}{T}_{th,1}^2;\sqrt{{\mathsf{\Gamma }}^{*}}={\displaystyle \frac{1}{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}}(\sqrt{\mathsf{\Gamma }}-T_{th,1})+T_{th,1}^{*}\hfill \end{array}$$

(11)

where ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ and ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ are the minimum and maximum eigenvalues of $\tilde{\mathsf{\Theta }}$. The block diagonal matrix, $\tilde{\mathsf{\Theta }}$, has two blocks. The first one is ${\mathsf{\Theta }}_2$ (see Equation (5)), and the second one is the identity matrix. Hence, ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ = 1, and ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ = ${\lambda }_{\mathsf{\Theta },2}^{\min }$ (minimum eigenvalue of ${\mathsf{\Theta }}_2$), which is shown to be equal to the minimum eigenvalue of $\mathsf{\Omega }$, ${\lambda }_{\mathsf{\Omega }}^{\min }$. It is also proved that all eigenvalues of $\mathsf{\Omega }$ are real and positive. First, ${\mathsf{\Gamma }}^{*}$ is determined such that the probability of the non-central chi-squared distribution $P_{ncx}({\left(T_{th,1}^{*}\right)}^2,\mathrm{DOF},{\mathsf{\Gamma }}^{*})$ = $P_{MD,i}^{alloc}$; DOF stands for the required degrees of freedom. $P_{MD,i}^{alloc}$ is the allocated probability of missed detection for the ith fault mode. Then, $\mathsf{\Gamma }$ is calculated from Equation (11) with ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ = 1 and saved in a lookup table for different values of DOF, $P_{MD,i}^{alloc}$, and $1/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$. Following [30], it is also justified that $\mathsf{\Gamma }$ increases as ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ decreases. $\mathsf{\Gamma }$ for ${\alpha }_1$ is obtained from the lookup table for a fault mode, the next entry of $1/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ higher than the calculated value, and the required DOF.

$VP{E}_{2,i}^U$ for ${\alpha }_2$ is calculated similarly to that of $VP{E}_{1,i}^U$ but using a more computationally efficient approach. This is because it involves terms of $(N-M)$ epochs, often found to be greater than M. ${VPE}_3^U$ is the mean vertical position error (VPE) upper bound for the ignored terms and the same for all fault modes. It is

$${VPE}_3^U=U_{\max }{N}_{\max }\times {\displaystyle \frac{1\times {10}^{-6}}{1-0.1}}({\widehat{b}}_{\max }+{\dot{\widehat{b}}}_{\max })$$

(12)

In order to have the north and east components, the factor $1\times {10}^{-6}$ is changed to $\sqrt{2}\times {10}^{-6}$ for the horizontal position. The maximum of all Ns up to the current epoch $t_k$ is denoted as $N_{\max }$. $U_{\max }$ is an upper bound obtained in a computationally efficient way at each epoch [26]. The maximum bias ${\widehat{b}}_{\max }$ at epoch k is given as

$$\begin{array}{cc}\hfill \left(\underset{\forall \ell =1,\dots ,k}{min}\underset{{\alpha }_1,{\alpha }_2}{min}\underset{i\in \mathrm{fault}\mathrm{modes}}{min}{\lambda }_{\mathsf{\Theta },3}^{\min }\right){\widehat{b}}_{\max }^2& ={\mathsf{\Gamma }}_{\max }\hfill \end{array}$$

(13)

where ${\lambda }_{\mathsf{\Theta },3}^{\min }$ is the minimum eigenvalue of ${\mathsf{\Theta }}_3$. Operation “min” represents the minimum of all terms underneath it. ${\mathsf{\Gamma }}_{\max }$ is the maximum of all calculated $\mathsf{\Gamma }$s. The maximum bias rate, ${\dot{\widehat{b}}}_{\max }$, can be determined in a similar way with the minimum eigenvalue of ${\mathsf{\Theta }}_4$. ${\mathsf{\Theta }}_4$ is a block-diagonal matrix, with a block having a dimension of one or two. It can be shown that for each two-dimensional block, the minimum eigenvalue is more than the ratio of its determinant to the sum of the diagonal elements. This lower bound of the minimum eigenvalue is used to find ${\dot{\widehat{b}}}_{\max }$. Next, the PL equations are briefly mentioned.

2.2.2. PLs

Equations for the vertical PL (VPL) are given as follows. They can be easily modified for the horizontal PL (HPL). Under the SF modes, VPL is

$${\mathrm{VPL}}_{SF}=\sum _{j=1}^2\left(\underset{i\in \mathrm{SF}\mathrm{modes}}{max}VP{E}_{j,i}^U\right)+VP{E}_3^U+{\gamma }_s{\sigma }_v$$

(14)

The operation “max” represents the maximum of all terms underneath it. The probability $P(\mathrm{absolute}\mathrm{VPE}>{\gamma }_s{\sigma }_v)$ = $P_{MD,sf}^{alloc}$. ${\sigma }_v$ is the standard deviation of the error in the estimated vertical position. $P_{MD,sf}^{alloc}$ is the allocated probability of missed detection for SF modes. Similarly, the VPL under DF modes is calculated. The final VPL is

$$\mathrm{VPL}=max({\mathrm{VPL}}_{NF},{\mathrm{VPL}}_{SF},{\mathrm{VPL}}_{DF})$$

(15)

In the preceding equation, VPL${}_{NF}$ stands for the VPL under no fault conditions.

2.3. Limitations of Existing Approach and Modifications

The limitations of the described KF RAIM algorithm are as follows:

(1)

SKF does not estimate user acceleration needed for moderate to high dynamics.

(2)

Time constants of all time-correlated errors are assumed to be 100 s. More appropriate error models having different time constants for different error sources must be adopted.

(3)

$a_{\ell 1}$ entries of ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ are calculated with modeled values of measurement error standard deviations. By definition, they should be obtained with true error standard deviations, which are not known. A more accurate model of ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ would therefore account for uncertainties in the knowledge of actual $a_{\ell 1}$.

(4)

In the beginning of the receiver operation, relatively large transient position errors of the navigation filter can be handled by including artificial white noise in the carrier-smoothed pseudorange measurement error model of the SKF. While this can be easily incorporated in the SKF, the $\mathsf{\Psi }$ matrix of $\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}$ does not model white noise.

Thus, $\mathsf{\Psi }$ needs to be appropriately modeled, considering more realistic error models, initial artificial white noise, and uncertainties in the knowledge of $a_{\ell 1}$. The next two sections discuss the adopted error models in the current work and appropriate bounds of the $\mathsf{\Psi }$ matrix. SKF formulation also estimates error in the a priori user acceleration from truth. N is adaptively calculated the same way as before with $\mathbf{B}$ matrices. M is also calculated as before, but with ${\mathbf{B}}^{″}$ matrices, where ${\mathbf{B}}_k^{″}$(1:2:5,:) = ${\mathbf{B}}_k^{\prime }$ = ${\mathbf{C}}_e^n{\mathbf{B}}_k$(1:3:7,:), and ${\mathbf{B}}_k^{″}$(2:2:6,:) = ${\mathbf{C}}_e^n{\mathbf{B}}_k$(2:3:8,:). ${\mathbf{B}}_k$(1:3:7,:) has rows of ${\mathbf{B}}_k$ corresponding to the position states of the filter, and ${\mathbf{B}}_k$(2:3:8,:) has rows for velocity states. This is found to reduce M and hence PLs. The effect of overbounding the pseudorange rate (or precise range rate) error covariance is also included, and a tighter value of $\mathsf{\Gamma }$ in Equation (11) is determined.

3. Modified Pseudorange Error Models and $\mathsf{\Psi }$

In this section, first, the adopted error models of the carrier-smoothed pseudorange measurements with different time constants are explained. Following this, it is briefly discussed how they are incorporated in the SKF navigation filter. Finally, necessary modifications of ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ for satellite i are presented. It is assumed that the Ionospheric error is eliminated with dual frequency measurements. The following models are adopted for other error sources with dual frequency correction terms, when needed.

Residual zenith tropospheric delay after applying corrections of the MOPS model [31] has a latitude- and longitude-dependent bias, as noted in [32]. It can be corrected in the pseudorange measurements using a lookup table as a function of latitude and longitude. The remaining zenith error, ${\underline{\mathit{\eta }}}_{\mathrm{tr}}\left(i\right)$, is a zero-mean FOGM process with a standard deviation (${\sigma }_{\mathrm{tr}}$) of 0.2 m and time constant (${\tau }_{\mathrm{tr}}$) of 25 h for any visible satellite i. The broadcast ephemeris and satellite clock error, which is overbounded by the user range accuracy (URA), is denoted by ${\underline{\mathit{\eta }}}_{\mathrm{ec}}\left(i\right)$. It is modeled as an FOGM process, having a zero mean, unit standard deviation (${\sigma }_{\mathrm{ura}}$), and 20 h time constant (${\tau }_{\mathrm{ura}}$). The time between two effectively independent range error samples for satellite clock and orbit error is considered 50 h [33], which is taken as approximately 2.5 × the time constant. It is justified in [34] that common errors of code and carrier measurements remain unchanged in the carrier-smoothed code, thus allowing the same models to be used. As multipath and receiver noise are application specific, they are modeled for safety critical aviation applications in this paper. Their standard deviations, denoted as ${\sigma }_{\mathrm{ml},i}$ and ${\sigma }_{\mathrm{n},i}$, respectively, are obtained from [27] for airborne users. In addition, the constant term of the receiver noise standard deviation is inflated further to bound the autocorrelation function of carrier-smoothed pseudoranges due to multipath and noise derived in [35] for airborne users. The inflated value is 0.23 m in place of 0.15 m. Finally, multipath and receiver noise together (${\underline{\mathit{\eta }}}_{\mathrm{ml}-\mathrm{n}}\left(i\right)$) are represented as a zero-mean FOGM process with a time constant (${\tau }_{\mathrm{ml}-\mathrm{n}}$) of 100 s and standard deviation $\sqrt{{\sigma }_{\mathrm{ml},i}^2+{\sigma }_{\mathrm{n},i}^2}$. In the current work, NavIC is considered the second constellation. The same models as those of GPS are adopted for it, but with frequency bands appropriately changed. The approach described in the paper also holds for other constellations with model parameters suitably changed.

The SKF navigation filter estimates error in the a priori position, velocity and acceleration of the user in ECEF coordinates, and receiver clock offsets and drifts for GPS and NavIC from respective truths. The states of the aforementioned FOGM processes for each visible satellite are the “consider parameters” of the SKF. The zenith tropospheric error is projected onto the line of sight (LOS) using the mapping function of the MOPS model in the measurement model matrix. In the simulation studies, the position error of the SKF is compared with that of an EKF that estimates all FOGM states to depict comparable performance. Since the SKF takes into account the effects of the “consider parameters” in its estimates, no additional terms other than those in Equation (14) are needed for PL calculation.

With modified error models, the $a_{\min }$ and $a_{\max }$ mentioned after Equation (6) for test statistics are re-calculated as follows. If $\underline{\mathit{\eta }}\left(i\right)$ is the carrier-smoothed pseudorange measurement error of satellite i, then it constitutes the following at an epoch ℓ:

$${\underline{\mathit{\eta }}}^{(\ell )}\left(i\right)={\mathrm{mp}}^{(\ell )}\left({\theta }^{(\ell )}\left(i\right)\right){\underline{\mathit{\eta }}}_{\mathrm{tr}}\left(i\right)+{\underline{\mathit{\eta }}}_{\mathrm{ec}}\left(i\right)+{\underline{\mathit{\eta }}}_{\mathrm{ml}-\mathrm{n}}^{(\ell )}\left(i\right)$$

(16)

where the superscript with parentheses denotes time epochs, and mp($\theta$) is the mapping function of the MOPS model for elevation angle $\theta$. The autocorrelation between ${\underline{\mathit{\eta }}}^{(\ell )}\left(i\right)$ and ${\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)$ is

$$\begin{array}{cc}E\left({\underline{\mathit{\eta }}}^{(\ell )}\left(i\right){\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)\right)\hfill & ={\left({\sigma }^{(\ell )}\left(i\right)\right)}^2{exp}^{(-\beta q\mathsf{\Delta }t)}={\mathrm{mp}}^{(\ell )}\left({\theta }^{(\ell )}\left(i\right)\right){\mathrm{mp}}^{(\ell +q)}\left({\theta }^{(\ell +q)}\left(i\right)\right){\sigma }_{\mathrm{tr}}^2{exp}^{(-q\mathsf{\Delta }t/{\tau }_{\mathrm{tr}})}+\hfill \\ \hfill & {\sigma }_{\mathrm{ura}}^2{exp}^{(-q\mathsf{\Delta }t/{\tau }_{\mathrm{ura}})}+({\left({\sigma }_{\mathrm{ml},i}^{(\ell )}\right)}^2+{\left({\sigma }_{\mathrm{n},i}^{\ell )}\right)}^2){exp}^{(-q\mathsf{\Delta }t/{\tau }_{\mathrm{ml}-\mathrm{n}})}\hfill \end{array}$$

(17)

where $\mathsf{\Delta }t$ = 1 s. E represents the expectation operation. ${\sigma }^{(\ell )}\left(i\right)$ is the standard deviation of ${\underline{\mathit{\eta }}}^{(\ell )}\left(i\right)$. $\beta$ denotes the reciprocal of the time constant of an equivalent FOGM process for $\underline{\mathit{\eta }}\left(i\right)$. ${\left({\sigma }^{(\ell )}\left(i\right)\right)}^2$ is

$${\left({\sigma }^{(\ell )}\left(i\right)\right)}^2={\left({\mathrm{mp}}^{(\ell )}\left({\theta }^{(\ell )}\left(i\right)\right)\right)}^2{\sigma }_{\mathrm{tr}}^2+{\sigma }_{\mathrm{ura}}^2+({\left({\sigma }_{\mathrm{ml},i}^{(\ell )}\right)}^2+\left({\sigma }_{\mathrm{n},i}^{(\ell )}{)}^2\right)$$

(18)

Thus, based on the correlation between epochs ℓ and $(\ell +q)$, $\beta$ (denoted as ${\beta }_{\ell q}$) is calculated as

$$\begin{array}{c}\hfill {\beta }_{\ell q}=-log\left({\displaystyle \frac{E\left({\underline{\mathit{\eta }}}^{(\ell )}\left(i\right){\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)\right)}{{\left({\sigma }^{(\ell )}\left(i\right)\right)}^2}}\right)/\left(q\mathsf{\Delta }t\right)\end{array}$$

(19)

It is evident from Equations (17) and (19) that ${\beta }_{\ell q}$ is a function of q and $\theta$. It should be noted that ${\left({\sigma }_{\mathrm{ml},i}^{(\ell )}{)}^2+({\sigma }_{\mathrm{n},i}^{(\ell )}\right)}^2$ is also a function of $\theta$.

Following the definition of $\mathsf{\Psi }$ for normalized innovations $\mathsf{\Delta }{\underline{\tilde{\mathit{\rho }}}}^{*}$, ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ for satellite i mentioned after Equation (6) can be re-written as $\left[\begin{array}{ccccc}1& a_{11}& a_{12}& a_{13}& \cdots \\ a_{11}& 1& a_{21}& a_{22}& \cdots \\ a_{12}& a_{21}& 1& a_{31}& \cdots \\ & & ⋮& \end{array}\right]$, where $a_{\ell q}$ for an epoch ℓ with $q=1,2,\dots$ is defined as

$$a_{\ell q}={\displaystyle \frac{E\left({\underline{\mathit{\eta }}}^{(\ell )}\left(i\right){\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)\right)}{{\sigma }^{(\ell )}\left(i\right){\sigma }^{(\ell +q)}\left(i\right)}}={\left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}\right)}^q{\displaystyle \frac{{\left({\sigma }^{(\ell )}\left(i\right)\right)}^2}{{\sigma }^{(\ell )}\left(i\right){\sigma }^{(\ell +q)}\left(i\right)}}$$

(20)

$a_{\ell q}$ can be rearranged as

$$a_{\ell q}=\left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }^{(\ell +1)}\left(i\right)}}\right)\left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +1)}\left(i\right)}{{\sigma }^{(\ell +2)}\left(i\right)}}\right)\cdots \left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +q-1)}\left(i\right)}{{\sigma }^{(\ell +q)}\left(i\right)}}\right)$$

(21)

Representing the terms within the parentheses as $a_{\iota 1,\ell q}$, $a_{\ell q}$ = ${\prod }_{\iota =\ell }^{\iota =\ell +q-1}{a}_{\iota 1,\ell q}$. $\ell q$ in the subscript of $a_{\iota 1,\ell q}$ denotes that $\beta$ is obtained by considering the correlation between two epochs ℓ and $\ell +q$. Maximum $a_{\iota 1,\ell q}$ over $\iota =\ell ,\dots ,\ell +q-1$ is denoted as $a_{\iota 1,\ell q}^{\max }$, and the minimum $a_{\iota 1,\ell q}$ is $a_{\iota 1,\ell q}^{\min }$.

$a_{\iota 1,\ell q}$s are calculated with ${\beta }_{\ell q}$ of Equation (19). Change in $\theta$ between two consecutive epochs, $({\theta }^{(\iota +1)}\left(i\right)-{\theta }^{\left(\iota \right)}\left(i\right))$, lies between $\pm {\displaystyle \frac{v_{\mathrm{m}}^{sat}+v_{\mathrm{m}}^{user}}{r_{su}^{\left(\iota \right)}}}\times {\displaystyle \frac{180}{\pi }}$, where $v_{\mathrm{m}}^{sat}$ = maximum satellite speed, $v_{\mathrm{m}}^{user}$ = maximum user speed = Mach 1, $r_{su}^{\left(\iota \right)}$ = range between the user and satellite at ${\theta }^{\left(\iota \right)}$ calculated using the approach of [36], and maximum user altitude = 15 km. With M and $(N-M)$ less than or equal to 20, three different values of q are chosen for the calculation of ${\beta }_{\ell q}$. With maximum satellite and user speeds and user altitude, equivalent time constants, $1/{\beta }_{\ell q}$, and $a_{\iota 1,\ell q}^{\max }$ are plotted as a function of ${\theta }^{(\ell )}$ with increasing and decreasing $\theta$ over $t_{\ell },\dots ,t_{\ell +q}$ and presented in Figure 1 and Figure 2, respectively, for both GPS and NavIC. $a_{\iota 1,\ell q}^{\min }$ is not shown to avoid clutter. It can be obtained along with the time constants and $a_{\iota 1,\ell q}^{\max }$ using the Matlab program matlab_program_figures_1_2.m provided in the Supplementary Materials. It is apparent that at a given ${\theta }^{(\ell )}$, the variation of $a_{\iota 1,\ell q}^{\max }$ with q is very small. The maximum difference is 0.0004 between q = 1 and q = 19. For $a_{\iota 1,\ell q}^{\min }$, it is 0.0002. The selected $a_{\max }(={max}_{q,\theta }{a}_{\iota 1,\ell q}^{\max })$ and $a_{\min }(={min}_{q,\theta }{a}_{\iota 1,\ell q}^{\min })$ are shown with thick black lines in Figure 2.

Exhaustive simulations are also performed to validate $a_{\min }$ and $a_{\max }$ of Figure 2. This is described next. Uniformly random user speeds within [$-v_{\mathrm{m}}^{user}$$v_{\mathrm{m}}^{user}$] are selected over q epochs. An initial satellite speed is chosen by a uniformly random sample within [$-v_{\mathrm{m}}^{sat}$$v_{\mathrm{m}}^{sat}$]. Then, the satellite speed is propagated over q epochs with an acceleration randomly picked between its minimum and maximum limits. This way, four hundred thousand sets of user and satellite speeds are obtained. Half of these are for GPS, and the other half are for NavIC. For GPS, with each ${\theta }^{(\ell )}$ between 5° and 90° and each of three qs, each of the two hundred thousand sets of speeds is used to calculate ${\theta }^{(\ell +1)},\dots ,{\theta }^{(\ell +q)}$ for a given user altitude, and the corresponding $a_{\iota 1,\ell q}^{\min }$ and $a_{\iota 1,\ell q}^{\max }$ are noted. With half of the speeds, ${\theta }^{(\iota +1)}={\theta }^{\left(\iota \right)}+{\displaystyle \frac{v_{sat}^{\left(\iota \right)}+v_{user}^{{}^{(}\iota )}}{r_{su}^{\left(\iota \right)}}}\times {\displaystyle \frac{180}{\pi }},\iota =\ell ,\dots \ell +q-1$, where $v_{sat}^{\left(\iota \right)}$ and $v_{user}^{\left(\iota \right)}$ are satellite and user speeds at epoch $\iota$, and ${\theta }^{(\iota +1)}$ = ${\theta }^{\left(\iota \right)}-{\displaystyle \frac{v_{sat}^{\left(\iota \right)}+v_{user}^{\left(\iota \right)}}{r_{su}^{\left(\iota \right)}}}\times {\displaystyle \frac{180}{\pi }}$ with the other one hundred thousand sets. Considering all sets of data, $\left({max}_{q,\theta }{a}_{\iota 1,\ell q}^{\max }\right)$ is found to be less than the selected $a_{\max }$ and $\left({min}_{q,\theta }{a}_{\iota 1,\ell q}^{\min }\right)>a_{\min }$. The same is repeated for NavIC with the other two hundred thousand sets. This is performed numerous times with the Matlab program matlab_program_figures_1_2.m of the Supplementary Materials. Next, it is repeated for each user altitude from 0 to 15 km in steps of 5 km. Finally, in each simulation, the user altitude is varied randomly with uniform distribution between 0 and 15 km over q epochs. In all simulations, $a_{\max }$ and $a_{\min }$ of Figure 2 are found to be the correct upper and lower bounds.

The structure of the corresponding ${\mathsf{\Psi }}_{\mathrm{sub},\max }^i$ remains the same. That is, ${\mathsf{\Psi }}_{\mathrm{sub},\max =}^i\left[\begin{array}{ccccc}1& a_{\max }& a_{\max }^2& a_{\max }^3& \cdots \\ a_{\max }& 1& a_{\max }& a_{\max }^2& \cdots \\ a_{\max }^2& a_{\max }& 1& a_{\max }& \cdots \\ & & ⋮& \end{array}\right]$.

Next, modifications of $\mathsf{\Psi }$ due to uncertainties in the knowledge of actual $a_{\ell q}$ and initial artificial white noise are provided.

4. $\mathsf{\Psi }$ with Uncertain Parameters and White Noise

In this section, $\mathsf{\Psi }$ (or ${\mathsf{\Psi }}_{\mathrm{sub}}^i$) is first appropriately modeled with underlying parameters within a range to include the effect of uncertainties in their knowledge, and to account for initial artificial white noise. Subsequently, with the help of extensive simulations, it is shown that ${\mathsf{\Psi }}_{\max }$ can be used in place of $\mathsf{\Psi }$ to form test statistics and PLs.

4.1. Uncertain $a_{\ell q}$

Uncertainties in the knowledge of $a_{\ell q}$ arise because actual values of the parameters needed to calculate it are not known. Suppose the actual value of ${\tau }_{\mathrm{tr}}$ for a satellite is ${\varrho }_{\mathrm{tr}}$ times 25 h (modeled value), where ${\varrho }_{\mathrm{tr},\min }\le {\varrho }_{\mathrm{tr}}\le 1$. Similarly, the scale factor for ${\tau }_{\mathrm{ura}}$ is ${\varrho }_{\mathrm{ura}}$, which lies between ${\varrho }_{\mathrm{ura},\min }$ and one. Actual ${\sigma }_{\mathrm{tr}}^2$, ${\sigma }_{\mathrm{ura}}^2$ and $({\sigma }_{\mathrm{ml}}^2+{\sigma }_{\mathrm{n}}^2)$ also take on values below their respective modeled values. Their actual values can be obtained by multiplying scale factors ${\mathcal{K}}_{\mathrm{tr}}$, ${\mathcal{K}}_{\mathrm{ura}}$, and ${\mathcal{K}}_{\mathrm{ml}-\mathrm{n}}$, respectively, to the respective modeled values. The lower limits of these scale factors are ${\mathcal{K}}_{\mathrm{tr},\min }$, ${\mathcal{K}}_{\mathrm{ura},\min }$, and ${\mathcal{K}}_{\mathrm{ml}-\mathrm{n},\min }$, respectively, and the upper limits are one. All lower limits should be determined by prior analysis. A possible approach would be to analyze several months of data, following [32,33] for the tropospheric error and URA, respectively. The autocorrelation functions, power spectral densities, and error distributions of the data can be formed to determine both upper and lower bounds of the time constants and variances. The ratios of the lower to upper bounds will give the lower limits of the scale factors, and the upper bounds will be the modeled values discussed earlier. Although not considered in this paper, modeled values may change across satellites. Similar analyses may be done for multipath and receiver noise, but this will be application specific. Equation (17) is now re-written as

$$\begin{array}{cc}E\left({\underline{\mathit{\eta }}}^{(\ell )}\left(i\right){\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)\right)\hfill & ={\left({\sigma }^{(\ell )}\left(i\right)\right)}^2{exp}^{(-\beta q\mathsf{\Delta }t)}={\mathrm{mp}}^{(\ell )}\left({\theta }^{(\ell )}\left(i\right)\right){\mathrm{mp}}^{(\ell +q)}\left({\theta }^{(\ell +q)}\left(i\right)\right){\mathcal{K}}_{\mathrm{tr}}{\sigma }_{\mathrm{tr}}^2{exp}^{(-q\mathsf{\Delta }t/\left({\varrho }_{\mathrm{tr}}{\tau }_{\mathrm{tr}}\right))}+\hfill \\ \hfill & {\mathcal{K}}_{\mathrm{ura}}{\sigma }_{\mathrm{ura}}^2{exp}^{(-q\mathsf{\Delta }t/\left({\varrho }_{\mathrm{ura}}{\tau }_{\mathrm{ura}}\right))}+{\mathcal{K}}_{\mathrm{ml}-\mathrm{n}}({\left({\sigma }_{\mathrm{ml},i}^{(\ell )}\right)}^2+{\left({\sigma }_{\mathrm{n},i}^{\ell )}\right)}^2){exp}^{(-q\mathsf{\Delta }t/{\tau }_{\mathrm{ml}-\mathrm{n}})}\hfill \end{array}$$

(22)

In this work, no uncertainty in the time constant of multipath and noise is assumed. However, this can be relaxed if needed. Using Equation (22), modified ${\beta }_{\ell q}$, ${\beta }_{\ell q}^{\mathrm{mod}}$ is

$${\beta }_{\ell q}^{\mathrm{mod}}=-log\left({\displaystyle \frac{E\left({\underline{\mathit{\eta }}}^{(\ell )}\left(i\right){\underline{\mathit{\eta }}}^{(\ell +q)}\left(i\right)\right)}{{\left({\sigma }^{(\ell )}\left(i\right)\right)}^2}}\right)/\left(q\mathsf{\Delta }t\right)$$

(23)

Comparing Equations (17) and (22), it is evident that ${exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}q\mathsf{\Delta }t)}$≤${exp}^{(-{\beta }_{\ell q}q\mathsf{\Delta }t)}$. The actual ${\sigma }^{(\ell )}{\left(i\right))}^2$, ${\sigma }_{\mathrm{act}}^{(\ell )}{\left(i\right))}^2$ is given by

$${\sigma }_{\mathrm{act}}^{(\ell )}{\left(i\right))}^2={\left({\mathrm{mp}}^{(\ell )}\left({\theta }^{(\ell )}\left(i\right)\right)\right)}^2{\mathcal{K}}_{\mathrm{tr}}{\sigma }_{\mathrm{tr}}^2+{\mathcal{K}}_{\mathrm{ura}}{\sigma }_{\mathrm{ura}}^2+{\mathcal{K}}_{\mathrm{ml}-\mathrm{n}}\left({\left({\sigma }_{\mathrm{ml},i}^{(\ell )}\right)}^2+{\left({\sigma }_{\mathrm{n},i}^{\ell )}\right)}^2\right)$$

(24)

As Equation (20) has the modeled values of $\sigma$ in the denominator, it reduces to

$$a_{\ell q}={\left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}\right)}^q{\displaystyle \frac{{\left({\sigma }^{(\ell )}\left(i\right)\right)}^2}{{\sigma }^{(\ell )}\left(i\right){\sigma }^{(\ell +q)}\left(i\right)}}$$

(25)

With the rearrangement of terms as before,

$$a_{\ell q}=\left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }^{(\ell +1)}\left(i\right)}}\right)\cdots \left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +q-1)}\left(i\right)}{{\sigma }^{(\ell +q)}\left(i\right)}}\right)$$

(26)

The ℓth diagonal term of ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ is

$$a_{\ell 0}={\left({\displaystyle \frac{{\sigma }_{\mathrm{act}}^{(\ell )}\left(i\right)}{{\sigma }^{(\ell )}\left(i\right)}}\right)}^2$$

(27)

It should be noted that the term $({\sigma }_{\mathrm{act}}^{(\ell )}\left(i\right)/{\sigma }^{(\ell )}\left(i\right))$ is less than or equal to one as ${\sigma }^{(\ell )}\left(i\right)$≥${\sigma }_{\mathrm{act}}^{(\ell )}\left(i\right)$. If ${\mathcal{K}}_{\min }^{\prime }$ is the minimum of ${\mathcal{K}}_{\mathrm{tr},\min }$, ${\mathcal{K}}_{\mathrm{ura},\min }$, and ${\mathcal{K}}_{\mathrm{ml}-\mathrm{n},\min }$, then using Equation (24), ${\mathcal{K}}_{\min }^{\prime }$ ≤ $a_{\ell 0}$ ≤ 1.

Multiplying and dividing ${\left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}\right)}^q$ in Equation (26) and rearranging terms,

$$a_{\ell q}=\left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }^{(\ell +1)}\left(i\right)}}\right)\cdots \left({exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +q-1)}\left(i\right)}{{\sigma }^{(\ell +q)}\left(i\right)}}\right){\left({\displaystyle \frac{{exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}}{{exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}}}\right)}^q$$

(28)

Using $a_{\min }$ and $a_{\max }$, and ${exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}q\mathsf{\Delta }t)}$≤${exp}^{(-{\beta }_{\ell q}q\mathsf{\Delta }t)}$ noted after Equation (23),

$$a_{\min }^q\left(\underset{q,\theta \mathcal{K},\varrho }{min}{\left({\displaystyle \frac{{exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}}{{exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}}}\right)}^q\right)\le a_{\ell q}\le a_{\max }^q$$

(29)

where ${min}_{q,\theta \mathcal{K},\varrho }{\left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}/{exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)}\right)}^q$ is the minimum value obtained over all values of $\mathcal{K}$, $\varrho$, $\theta$, and q (maximum q is $(M-1)$ for ${\alpha }_1^{\prime }$ and $(N-M-1)$ for ${\alpha }_2^{\prime }$). For a given q and $\theta$, ${({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}/{exp}^{(-{\beta }_{\ell q}\mathsf{\Delta }t)})}^q$ is minimum when all $\mathcal{K}$s and $\varrho$s are set to their respective minimum values (see Equations (17), (19), (22), and (23)). This is greater than or equal to the value ${\mathcal{K}}_{\min }^{″}$ obtained by setting all $\mathcal{K}$s to ${\mathcal{K}}_{\min }^{\prime }$. ${\mathcal{K}}_{\min }^{″}$ is close but less than ${\mathcal{K}}_{\min }^{\prime }$. The minimum of all ${\mathcal{K}}_{\min }^{″}$s over all values of $\theta$ and q is set as the lower bound, and denoted as ${\mathcal{K}}_{\min }$ (0 < ${\mathcal{K}}_{\min }$ < 1). It can be determined by following the instructions in matlab_program_figures_1_2.m and then running it. Next, a provision is made to include artificial white noise in $a_{\ell q}$.

4.2. White Noise in $a_{\ell q}$

White noise is injected initially in the carrier-smoothed pseudorange measurement error model of the SKF to deal with transient errors in position estimates. During this initial period, the modeled ${\sigma }^{(\ell )}\left(i\right)$ is given by

$${\left({\sigma }_{\mathrm{init}}^{(\ell )}\left(i\right)\right)}^2={\left({\sigma }^{(\ell )}\left(i\right)\right)}^2+{\left({\sigma }_{\mathrm{wh}}^{(\ell )}\left(i\right)\right)}^2$$

(30)

If the injected white noise variance, ${\left({\sigma }_{\mathrm{wh}}^{(\ell )}\left(i\right)\right)}^2$, is ${\mathcal{K}}_{\mathrm{wh}}{\sigma }^{(\ell )}{\left(i\right))}^2$, ${\left({\sigma }_{\mathrm{init}}^{(\ell )}\left(i\right)\right)}^2$ becomes equal to $({\mathcal{K}}_{\mathrm{wh}}+1){\left({\sigma }^{(\ell )}\left(i\right)\right)}^2$. ${\mathcal{K}}_{\mathrm{wh}}$ is a known scale factor as white noise is artificially included. With this, Equation (26) is changed to

$$a_{\ell q}=\left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }^{(\ell +1)}\left(i\right)}}\right)\cdots \left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +q-1)}\left(i\right)}{{\sigma }^{(\ell +q)}\left(i\right)}}\right)\left({\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }_{\mathrm{init}}^{(\ell )}\left(i\right)}}\right)\left({\displaystyle \frac{{\sigma }^{(\ell +q)}\left(i\right)}{{\sigma }_{\mathrm{init}}^{(\ell +q)}\left(i\right)}}\right)$$

(31)

Replacing the last two terms on the right,

$$a_{\ell q}=\left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell )}\left(i\right)}{{\sigma }^{(\ell +1)}\left(i\right)}}\right)\cdots \left({exp}^{(-{\beta }_{\ell q}^{\mathrm{mod}}\mathsf{\Delta }t)}{\displaystyle \frac{{\sigma }^{(\ell +q-1)}\left(i\right)}{{\sigma }^{(\ell +q)}\left(i\right)}}\right){\displaystyle \frac{1}{({\mathcal{K}}_{\mathrm{wh}}+1)}}$$

(32)

Similarly, the ℓth diagonal term of ${\mathsf{\Psi }}_{\mathrm{sub}}^i$ is given by

$$a_{\ell 0}={\left({\displaystyle \frac{{\sigma }_{\mathrm{act}}^{(\ell )}\left(i\right)}{{\sigma }^{(\ell )}\left(i\right)}}\right)}^2{\displaystyle \frac{1}{({\mathcal{K}}_{\mathrm{wh}}+1)}}\le 1$$

(33)

Its lower bound is ${\mathcal{K}}_{\min }^{\prime }/({\mathcal{K}}_{\mathrm{wh}}+1)$. Using Equation (29), $a_{\ell q}$ is bounded by

$$a_{\min }^q{\displaystyle \frac{{\mathcal{K}}_{\min }}{({\mathcal{K}}_{\mathrm{wh}}+1)}}\le a_{\ell q}\le a_{\max }^q$$

(34)

4.3. Simulation Results

Having modeled $\mathsf{\Psi }$ with underlying parameters residing within a range, results of extensive simulations are shown to justify that ${\mathsf{\Psi }}_{\max }$ can be used in place of $\mathsf{\Psi }$ to determine test statistics and PLs. In other words, it is illustrated that ${\lambda }_{\mathsf{\Theta },2}^{\max }$ or ${\lambda }_{\mathsf{\Omega }}^{\max }$ with $\mathsf{\Theta }$ formed as ${(c\mathbf{I}+{\mathsf{\Psi }}_{\max })}^{-1}$ is less than or equal to one (see after Equation (5)). In addition, ${\lambda }_{\mathsf{\Theta },2}^{\min }$ or ${\lambda }_{\mathsf{\Omega }}^{\min }$ is shown to be more than the minimum eigenvalue of ${\tilde{\mathbf{S}}}^{\prime T}{\mathsf{\Psi }}_{\max }{\tilde{\mathbf{S}}}^{\prime }{\tilde{\mathbf{S}}}^{\prime \mathbf{T}}\mathsf{\Theta }{\tilde{\mathbf{S}}}^{\prime }$ (=${\mathsf{\Omega }}_1$) times a known scale factor. This is needed to ensure that the calculated $\mathsf{\Gamma }$ of Equation (11) is larger than its actual value (see Equation (10) and the paragraph after Equation (11) for justification of this). The actual value of $\mathsf{\Gamma }$ would need the actual $\mathsf{\Psi }$, which is not known.

For the simulation results of this section, lower limits of the scale factors used to account for uncertain parameters are: ${\varrho }_{\mathrm{tr},\min }$ = 0.5; ${\varrho }_{\mathrm{ura},\min }$ = 0.5; ${\mathcal{K}}_{\mathrm{tr},\min }$ = 0.25; ${\mathcal{K}}_{\mathrm{ura},\min }$ = 0.25; and ${\mathcal{K}}_{\mathrm{ml}-\mathrm{n},\min }$ = 0.25. For the selected lower limits, ${\mathcal{K}}_{\min }$ can be obtained a priori by running matlab_program_figures_1_2.m. ${\mathcal{K}}_{\mathrm{wh}}$ is 5 for white noise overbound. Scale factors are uniformly distributed between their lower and upper limits. While they are changed between simulations, for a given satellite they are assumed to be constant over M epochs (for ${\alpha }_1^{\prime }$) or $(N-M)$ epochs (for ${\alpha }_2^{\prime }$). M = 10, N = 20, and the time interval between two epochs $\mathsf{\Delta }t$ = 1 s. Fourteen GPS and seven NavIC satellites are considered. A brief description of the algorithm needed for simulations is provided in Appendix A. Matlab implementation of the algorithm can be found in the Supplementary Materials (matlab_program_algorithm$\_1\_\mathrm{ts}1\_$figure_3.m and matlab_program_algorithm$\_1\_\mathrm{ts}2\_$figure_3.m), where all inputs can be changed as needed, and results with the changed values can be generated.

The algorithm also has a provision to remove a satellite at the Mth or $(N-M)$th epoch to check the change in visibility and remove white noise injection from a random intermediate time epoch. Satellites are excluded when the elevation angle is below the 5° mask angle. A satellite can be chosen to come into visibility after a random number of epochs.

LOS vectors of $\mathbf{H}$ are randomly selected without considering actual GPS and NavIC satellite orbits. Hence, results of these simulations are likely to hold for other constellations if error model parameters, frequency bands, and $a_{\min }$ and $a_{\max }$ are suitably changed. Although there is an inherent relationship between elevation angles and $\mathbf{H}$, it is not important here as long as the structure of the simulated $\mathbf{H}$ resembles that of the actual measurement model matrix. That is, the LOS vector constraint of each row is satisfied, the clock bias columns are properly formed, and the rank of $\mathbf{H}$ is 5. With this, valid measurement model matrices would just be a subset of the wide range of $\mathbf{H}$ matrices considered in the simulations.

Figure 3 depicts the difference between the maximum eigenvalues of ${\mathsf{\Omega }}_1$ and $\mathsf{\Omega }$ for both ${\alpha }_1^{\prime }$ and ${\alpha }_2^{\prime }$ when $\theta$ decreases with time (see subplots 1 and 3—one below the other on the left). ${\lambda }_{\mathsf{\Omega },1}^{\min }$ and ${\lambda }_{\mathsf{\Omega },1}^{\max }$ are the minimum and maximum eigenvalues of ${\mathsf{\Omega }}_1$, respectively. In each simulation, the initial $\theta$ of a satellite is chosen with a uniform random variable between 5.1° and 90°. Satellite speed is constant at its maximum value. Five hundred thousand simulations are divided into ten batches for ease of plotting in Matlab, with each batch having 50,000 runs. The minimum and maximum differences over 50,000 simulations in a batch are shown as the end points of a vertical line. Within a batch, 10,000 runs have white noise injected up to some random intermediate epochs. Since the minimum difference between ${\lambda }_{\mathsf{\Omega },1}^{\max }$ and ${\lambda }_{\mathsf{\Omega }}^{\max }$ is greater than zero in all batches, ${\lambda }_{\mathsf{\Omega }}^{\max }$ is less than ${\lambda }_{\mathsf{\Omega },1}^{\max }$, which in turn is less than or equal to one by proper calculation of c. The differences between $1/({\lambda }_{\mathsf{\Omega },1}^{\min }{\mathcal{K}}_{\min })$ and $1/{\lambda }_{\mathsf{\Omega }}^{\min }$ for ${\alpha }_1^{\prime }$ (and between $1/({\lambda }_{\mathsf{\Omega },1}^{\min }{\displaystyle \frac{{\mathcal{K}}_{\min }}{({\mathcal{K}}_{\mathrm{wh}}+1)}})$ and $1/{\lambda }_{\mathsf{\Omega }}^{\min }$ for ${\alpha }_2^{\prime }$) are illustrated in the second (and fourth) subplot of the same figure. In this case, the longer vertical line of each batch belongs to the 10,000 simulations with artificial white noise. It is evident that $1/{\lambda }_{\mathsf{\Omega }}^{\min }$ < $1/({\lambda }_{\mathsf{\Omega },1}^{\min }{\mathcal{K}}_{\min })$ for ${\alpha }_1^{\prime }$, and $1/{\lambda }_{\mathsf{\Omega }}^{\min }$ < $1/({\lambda }_{\mathsf{\Omega },1}^{\min }{\displaystyle \frac{{\mathcal{K}}_{\min }}{({\mathcal{K}}_{\mathrm{wh}}+1)}})$ for ${\alpha }_2^{\prime }$. ${\lambda }_{\mathsf{\Omega }}^{\min }$ is verified to be always positive. It can also be concluded that $\sqrt{{\lambda }_{\mathsf{\Omega }}^{\max }}/{\lambda }_{\mathsf{\Omega }}^{\min }$ < $\sqrt{{\lambda }_{\mathsf{\Omega },1}^{\max }}/({\lambda }_{\mathsf{\Omega },1}^{\min }{\mathcal{K}}_{\min })$ for ${\alpha }_1^{\prime }$. Similarly, for ${\alpha }_2^{\prime }$, $\sqrt{{\lambda }_{\mathsf{\Omega }}^{\max }}/{\lambda }_{\mathsf{\Omega }}^{\min }$ < $\sqrt{{\lambda }_{\mathsf{\Omega },1}^{\max }}/({\lambda }_{\mathsf{\Omega },1}^{\min }{\displaystyle \frac{{\mathcal{K}}_{\min }}{({\mathcal{K}}_{\mathrm{wh}}+1)}})$. Hence, using ${\mathsf{\Psi }}_{\max }$ in place of $\mathsf{\Psi }$, ${\alpha }_1^{\prime }$, ${\alpha }_2^{\prime }$ and PLs can be calculated. Although not shown, the same conclusion can be derived with increasing elevation angle and changing satellite speeds, and verified with the provided Matlab programs.

It should be noted that when c is adaptively calculated for ${\alpha }_1^{\prime }$, ${\lambda }_{\mathsf{\Omega }}^{\max }$ is less than ${\lambda }_{\mathsf{\Omega },1}^{\max }$ and one only if c is more than or equal to a minimum value, $c_{\min }$. It is found that $c_{\min }$ is a function of M and elevation angle. Appendix B explains how $c_{\min }$ is obtained. During adaptive calculation, the initial guess of c is taken as $c_{\min }$ or higher. Next, a tighter value of $\mathsf{\Gamma }$ is determined. The effect of ${\mathbf{W}}_{rr,\ell }$ overbounding the actual pseudorange rate (or precise range rate) error covariance is also discussed.

5. Determination of $\mathsf{\Gamma }$

Substituting Equation (5) in the first term of Equation (6) on the right and having a similar representation for the second term, one can write

$${\alpha }_j^2={{\underline{\mathit{\eta }}}_r}^T{\mathsf{\Theta }}_2{\underline{\mathit{\eta }}}_r+{{\underline{\mathit{\eta }}}_{rr}}^T{\mathsf{\Theta }}_5{\underline{\mathit{\eta }}}_{rr}$$

(35)

where ${\underline{\mathit{\eta }}}_{rr}$ has unit-variance independent elements. ${\mathsf{\Theta }}_5$ would be the identity matrix only if ${\mathbf{W}}_{rr,\ell }$ of Equation (6) is equal to the actual pseudorange rate error covariance matrix. However, ${\mathbf{W}}_{rr,\ell }$ generally overbounds the actual error covariance. As a result, the maximum eigenvalue of ${\mathsf{\Theta }}_5$ is less than or equal to one, and the minimum eigenvalue is at some value greater than zero [29]. They are denoted as ${\lambda }_{\mathsf{\Theta },5}^{\max }$ and ${\lambda }_{\mathsf{\Theta },5}^{\min }$, respectively, and determined using prior knowledge of overbounding different satellite range rate errors. $\tilde{\mathsf{\Theta }}$ of Equation (11) thus changes to

$$\tilde{\mathsf{\Theta }}=\left[\begin{array}{cc}{\mathsf{\Theta }}_2& \left[\mathbf{0}\right]\\ \left[\mathbf{0}\right]& {\mathsf{\Theta }}_5\end{array}\right]$$

(36)

It is apparent that ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ = $max\{{\lambda }_{\mathsf{\Theta },2}^{\max },{\lambda }_{\mathsf{\Theta },5}^{\max }\}$ and ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ = $min\{{\lambda }_{\mathsf{\Theta },2}^{\min },{\lambda }_{\mathsf{\Theta },5}^{\min }\}$. The upper and lower bounds of ${\lambda }_{\mathsf{\Theta },2}^{\max }$ and ${\lambda }_{\mathsf{\Theta },2}^{\min }$, respectively, as determined in the previous section, are used in place of their actual, unknown values.

With the formulation of $\tilde{\mathsf{\Theta }}$ after Equation (11), where ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ = 1, $\mathsf{\Gamma }$ can be directly obtained from the lookup table. As $\mathsf{\Gamma }$ was shown to increase with decreasing ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$, the next entry of the lookup table immediately higher than the calculated $1/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ is considered. With $\tilde{\mathsf{\Theta }}$ of Equation (36), and the minimum and maximum eigenvalues defined after it, $\mathsf{\Gamma }$ is now calculated in the following way. In this case, ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$ is not necessarily equal to one. First, $\mathsf{\Gamma }$ is obtained from the lookup table the same way as before, using $P_{MD,i}^{alloc}$, DOF, and the next ratio entry higher than the calculated $\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$. This $\mathsf{\Gamma }$ is not the right one when ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$≠ 1, and therefore is denoted by ${\mathsf{\Gamma }}^{\prime }$. In addition, for the lookup table, the detection threshold, $T_{th,1}$, is now determined from the non-central chi-squared distribution with the required DOF and some predetermined natural bias [37]. Using ${\mathsf{\Gamma }}^{\prime }$, the ${\mathsf{\Gamma }}^{*}$ of Equation (11) is calculated as

$$\sqrt{{\mathsf{\Gamma }}^{*}}=\sqrt{{\mathsf{\Gamma }}^{\prime }}-T_{th,1}+T_{th,1}{\left({\displaystyle \frac{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }}}\right)}_{\mathrm{lookup}}$$

(37)

where ${\left({\displaystyle \frac{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }}}\right)}_{\mathrm{lookup}}$ is the lookup table entry. Using Equation (11), the correct $\mathsf{\Gamma }$ is given as

$$\sqrt{\mathsf{\Gamma }}=\left(\sqrt{{\mathsf{\Gamma }}^{*}}-T_{th,1}{\left({\displaystyle \frac{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }}}\right)}_{\mathrm{lookup}}\right)\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}+T_{th,1}$$

(38)

Substituting Equation (37) in the preceding equation

$$\sqrt{\mathsf{\Gamma }}=(\sqrt{{\mathsf{\Gamma }}^{\prime }}-T_{th,1})\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}+T_{th,1}$$

(39)

If $\sqrt{{\mathsf{\Gamma }}^{\prime }}$<$T_{th,1}$, $\sqrt{{\mathsf{\Gamma }}^{\prime }}$ is taken as $T_{th,1}$, which may occur at a very large DOF. As ${\left({\displaystyle \frac{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }}}\right)}_{\mathrm{lookup}}$ > calculated $\left({\displaystyle \frac{\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}}{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }}}\right)$, ${\mathsf{\Gamma }}^{\prime }$ is greater than the actual ${\mathsf{\Gamma }}^{\prime }$. Hence, the $\mathsf{\Gamma }$ from Equation (39) is more than the actual $\mathsf{\Gamma }$. The $\mathsf{\Gamma }$ corresponding to ${\alpha }_2$ is also obtained in the same way.

Performance and salient features of the modified KF RAIM are studied next for simulated scenarios of aircraft flight.

6. Simulation Studies

In this section, various aspects of KF RAIM are analyzed using simulated GPS L1 and NavIC L5 signals for aircraft flight. For this purpose, a GNSS simulator from Accord Software and Systems is used [38]. It has a provision to load GPS and NavIC Yuma files, and configure user motion profiles and required error models (e.g., ionospheric and tropospheric errors). In the simulation, the ionospheric error is modeled by the Klobuchar model and the tropospheric error by RTCA98 [31]. Simulated signals are received by an IFEN receiver front end and digitized at 20 MHz. Stored 20 MHz IF data are processed by a developed GPS L1-NavIC L5 software receiver. Figure 4 shows the signal simulation and data collection setup. Performance is first analyzed in detail for a segment of an aircraft flight path during descent and approach, which is illustrated in Figure 5 along with GPS-NavIC satellite visibility.

The IF data processing block diagram is depicted in Figure 6. A GPS L1-NavIC L5 software radio is developed in Matlab to perform all receiver operations. Following acquisition, it implements tracking of the visible satellites in parallel processing mode to speed up computation time. Tracking outputs are then sent to the navigation filter to estimate user position, velocity, and time. Relevant RAIM inputs are passed from the filter to the multi-constellation KF or WLS RAIM algorithm to carry out fault detection and PL computation.

The receiver removes the ionospheric error using the Klobuchar model but uses a different model (modified Hopfield [39]) for the tropospheric error. The residual zenith tropospheric error is represented by an FOGM process in the navigation filter with a zero mean, standard deviation of 0.2 m, and time constant of 25 h. It bounds the actual error in the simulated measurements. There is no residual ionospheric error. Satellite orbit and clock error and multipath were not included in the simulation. They are, however introduced in the carrier-smoothed pseudorange measurements following the respective FOGM models. As a result, with each run, different positioning error results can be obtained. Range rate measurements, on the other hand, are generated from the Doppler frequencies of the carrier phase tracking loops. Only white noise with standard deviations obtained from the tracking loops of the receiver is modeled for range rate measurements. ${\mathcal{K}}_{\min }$ for KF RAIM is taken as 0.75${}^2$ to be over conservative (see after Equation (34)).

There is a provision for three different navigation filters in the receiver: WLS, SKF, and full-state EKF. WLS estimates error in the a priori ECEF position and velocity of the aircraft, and the receiver clock offsets and drifts for GPS and NavIC from respective truths. It calculates the measurement error standard deviation considering the residual tropospheric error, URA, multipath, and noise. Full-state EKF and SKF estimate error in the a priori acceleration states of the aircraft in ECEF apart from those mentioned for WLS. In addition, the full-state EKF estimates three states for the measurement error of each visible satellite (residual zenith tropospheric error, satellite clock and orbit error, and multipath and noise) using the corresponding FOGM models. The measurement update interval is 1 s. EKF and SKF position estimation errors in the NED frame and carrier-smoothed pseudorange measurement innovations are shown in Figure 7 for a simulation run. As the SKF treats residual errors, multipath, and noise as “consider parameters”, its innovations remain biased, as opposed to those of full-state EKF.

6.1. Analyses of KF RAIM

In this paper, the following parameters are considered for analyses of KF RAIM. The continuity budget for false alert, $P_{FA}$ is $3.33\times {10}^{-7}$/sample [31]. The allocated probability of hazardously misleading information (HMI), $P_{HMI}^{alloc}$ is ${10}^{-7}$/h during aircraft descent and initial and intermediate approach, and ${10}^{-7}$/approach during the final approach [40]. The prior probability of narrow fault of GPS and NavIC satellites is ${10}^{-4}$/satellite/h. Assuming a one-hour mean time to notify by the constellation service provider, the prior probability of fault is ${10}^{-4}$/satellite as the simulation duration is less than an hour [41]. Constellation-wide faults are not taken into account. $P_{HMI}^{alloc}$ for the horizontal position, $P_{HMI,H}^{alloc}$, is 10% of $P_{HMI}^{alloc}$. $P_{HMI}^{alloc}$ for the vertical position, $P_{HMI,V}^{alloc}$, is the remaining 90%. For the prior fault probability and allocated integrity budget, at most two simultaneous independent satellite failures need to be monitored, as justified in [27]. The probability of three or more non-monitored faults $P_{nm}$, being less than the $P_{HMI}^{alloc}$, is subtracted from the latter. During the final approach phase, vertical and horizontal alert limits are considered to be 50 m and 40 m.

Next, the test statistics and PLs of KF RAIM obtained using the SKF are shown alongside those of WLS RAIM in Figure 8. WLS RAIM follows the method of [25] but accounts for the parameter uncertainty of the measurement error model, as discussed in [29]. It is used for comparison purposes for the following reasons: First, it performs fault detection in the range domain, as is with the developed KF RAIM. Second, it is computationally simple. The presented KF RAIM takes 5.5 min in Matlab on an old, low-end IBM Thinkpad laptop for the 19.48 min simulation. The laptop runs on Linux and has 1 GB RAM, 2 GB swap memory, and a Core 2 Duo processor (CPU @ 2 GHz). WLS RAIM, on the other hand, requires 18 s.

In order to gain a good understanding of the properties of KF RAIM, a few salient features are also analyzed. In this regard, Figure 9 depicts the following: M (number of terms of the first test statistic ${\alpha }_1$), N of Equation (7), maximum eigenvalues, ${\lambda }_{\mathsf{\Omega },1}^{\max }$ and DOF of the test statistics, and the eigenvalue ratios, $\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$. It is evident that ${\lambda }_{\mathsf{\Omega },1}^{\max }$ of the two statistics is less than one, as desired. Initially, eigenvalue ratios are large because of the addition of white noise during the transient period. The initial transient period is determined by checking the filter innovation sequence. As the eigenvalue ratios have different orders of magnitude, they are subsequently plotted on different scales, each given on either side. The figure also shows adaptive calculations of M and N and the corresponding changes in the DOF. It is apparent that M and N must be calculated adaptively as they may change over time.

Mean VPE bounds under faults are presented in Figure 10 for ${\alpha }_1$ and ${\alpha }_2$. They are plotted on different y-axes, each on either side. This shows separately the contribution of mean VPE bound corresponding to each test statistic. While N is generally much larger than M, the second bound with $(N-M)$ terms is much less than the first with M terms. Although not shown in the figure, the third bound from ignored terms before N recent epochs has an even smaller contribution. It is on the order of 0.02 m. That said, the contributions of all epoch measurements, including those before N, to the mean position error bound must be taken into account for PL calculation.

Next, $\sqrt{{\mathsf{\Gamma }}^{\prime }}$ defined before Equation (37) is plotted in Figure 11 from the saved lookup table against $\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ with DOF as a parameter. $\sqrt{{\mathsf{\Gamma }}^{\prime }}$ for $P_{MD}^{alloc}$ of both vertical and horizontal positions and dual fault modes is shown. For a given DOF, it levels off at a moderate eigenvalue ratio. Although the figure is plotted up to a ratio of 5000, the trend continues even afterwards. Similar patterns are also seen for the single fault mode. Thus, for a given DOF (or threshold) and ${\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }$, $\mathsf{\Gamma }$ calculated from ${\mathsf{\Gamma }}^{\prime }$ using Equation (39) would not change significantly with large ratios. Therefore, a large $\sqrt{{\lambda }_{\tilde{\mathsf{\Theta }}}^{\max }}/{\lambda }_{\tilde{\mathsf{\Theta }}}^{\min }$ does not imply that the PL would increase to a great extent.

In what follows, fault detection performance with three different navigation filters is discussed.

6.2. Fault Detection Performance

Three canonical faults—step, ramp, and sinusoidal—are chosen to analyze fault detection performance. It is explained in [3] that potential failure modes for GNSS precise positioning as complied from the existing literature can be largely modeled as one of the three fault types.

With full-state EKF, the test statistic ${\alpha }_{EKF,k}$ at time step k is formed as

$$\begin{array}{cc}\hfill {\alpha }_{EKF,k}& =\sum _{\ell =k}^{k-N+1}\left[\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }^T\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }^T\right]{\mathbf{R}}_{\ell }^{-1}\left[\begin{array}{c}\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }\\ \mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }\end{array}\right]\hfill \end{array}$$

(40)

where ${\mathbf{R}}_{\ell }$ is the innovation covariance matrix at time step ℓ. $\mathsf{\Delta }{\underline{\mathit{\rho }}}_{\ell }$ and $\mathsf{\Delta }{\underline{\dot{\mathit{\rho }}}}_{\ell }$ are carrier-smoothed pseudorange and range rate innovation vectors of the full-state EKF. N is determined the same way as that in Equation (7), but for the EKF. Under no fault conditions, ${\alpha }_{EKF,k}$ is chi-squared distributed. Its threshold is determined using $P_{FA}$ and DOF. As for WLS and SKF, thresholds are calculated from the non-central chi-squared distribution with the non-centrality parameter 25 for natural biases, following [37]. Test statistics ${\alpha }_1$ and ${\alpha }_2$ each have the probability of false alarm $P_{FA}/2$.

From the PL calculation of WLS RAIM and KF RAIM with SKF, it is noted that simultaneous faults in PRN 28 of GPS (G28 in the skyplot) and PRN 2 of NavIC (N2 in the skyplot) result in maximum mean VPE (or FMS). Hence, faults are injected into these satellites (either N2 or (G28 and N2)) at 100 s into the simulation. In the case of dual fault modes, both satellites have the same fault.

Table 1. Fault detection performance of SKF and full-state EKF. Fault detection time is calculated from the beginning of simulation. PE stands for position error. It can be horizontal or vertical PE. The corresponding PL is HPL or VPL. Maximum of horizontal PE and VPE is noted under PE. “sat” stands for satellite.

Fault Type	SKF			Full-State EKF
	Fault Detection Time		If PE > PL	Fault Detection Time		Change in PE
	1 sat (N2)	2 sat (G28 and N2)		1 sat (N2)	2 sat (G28 and N2)
−5 m step	No	No	No; constant offset	100 s	100 s	Constant offset
	detection	detection	1 sat: PE < 5 m			1 sat: PE < 6.5 m
			2 sat: PE < 10 m			2 sat: PE < 8.55 m
Ramp			Yes	No	No	PE
with slope	563 s	614 s	1 sat: from 1100 s	detection	detection	grows with
−0.05 m/s			2 sat: from 661 s			time
Sinusoidal	No	No	No; sinusoidal PE	No	No	Sinusoidal PE
5$sin(2\pi t/3600)$ m	detection	detection	1 sat: PE < 5 m	detection	detection	1 sat: PE < 6.5 m
			2 sat: PE < 10 m			2 sat: PE < 10 m
Sinusoidal	No	No	No; sinusoidal PE	100 s	100 s	Sinusoidal PE
5$cos(2\pi t/3600)$ m	detection	detection	1 sat: PE < 5 m			1 sat: PE < 6.5 m
			2 sat: PE < 10 m			2 sat: PE < 10 m

and

Table 2. Fault detection performance of WLS. Fault detection time is calculated from the beginning of simulation.

Fault Type	WLS
	Fault Detection Time		If PE > PL
	1 sat (N2)	2 sat (G28 and N2)
−5 m step	No	No	No; constant offset;
	detection	detection	1 sat: PE < 5 m
			2 sat: PE < 10 m
Ramp			Yes
with slope	495 s	534 s	1 sat: from 943 s
−0.05 m/s			2 sat: from 584 s
Sinusoidal	No	No	No; sinusoidal PE
5$sin(2\pi t/3600)$ m	detection	detection	1 sat: PE < 5 m
			2 sat: PE < 10 m
Sinusoidal	No	No	No; sinusoidal PE
5$cos(2\pi t/3600)$ m	detection	detection	1 sat: PE < 5 m
			2 sat: PE < 8 m

summarize the fault detection performance for a particular simulation run. Since PLs are not computed for full-state EKF, its change in position error under faulty conditions is only noted. For the –5 m step fault and cosine fault of 5 m amplitude, both ${\alpha }_1$ with SKF and ${\alpha }_{EKF}$ exhibit an initial peak, but the test statistic of WLS does not. However, the peak of the EKF test statistic crosses its threshold, whereas that with the SKF does not. As a result, no detection of faults occurs with the test statistics of SKF. A possible reason for non-detection is that the magnitude of the faults is close to that of the innovation sequences for SKF, as is evident from Figure 7. It is also important to note that the ramp fault is detected later by both WLS and SKF test statistics when two satellites are faulty as opposed to a single satellite. This is probably because faults in G28 and N2 maximize the FMS for mean VPE, and therefore they take longer times to become detected. The WLS test statistic detects the ramp fault before that of SKF. However, the position error of the latter crosses its PL later.

Next, it is checked if the SKF test statistics can detect step faults in G28 and N2, for which the VPE just crosses the VPL in different simulation runs and at different time instants. It is observed that its first statistic ${\alpha }_1$ exhibits a sharp peak at the onset time of faults. The magnitude of the peak does cross the threshold and is near 45–50, resulting in instant detection. Similar analysis is also performed with WLS. It detects the corresponding step faults within 10 s of PE crossing the PL during initial aircraft descent and before the crossing afterwards. Thus, detection occurs within the required time to alert [41]. The same observation is made with cosine faults.

Table 1 indicates that the full-state EKF does not detect the ramp faults, but its position error grows with time. With dual faults, the position error increases faster. Although faults were introduced into two satellites of GPS and NavIC—one each—as they maximize the FMS for mean VPE, it is unlikely that satellites from different constellations would fail simultaneously. Therefore, it is investigated what would occur when two satellites of NavIC fail. With NavIC being a new constellation, such a possibility cannot be ruled out. It is noted that if N2 and N6 fail at 10 s into simulation, both with a ramp fault of −0.05 m/s, the EKF VPE crosses the vertical alert limit during the final approach phase without any detection. This is illustrated in Figure 12. The performance of KF RAIM with SKF is also shown for comparison purposes. The EKF test statistic does not display any abnormal behavior because the filter confuses such a slowly growing fault with a residual measurement error that it is supposed to estimate. SKF can detect such faults because it does not estimate measurement errors, but treats them as “consider parameters” instead. Missed detection with full-state EKF is also observed in all of 100 different simulation runs. This occurs with positive slopes as well as different slopes in two satellites.

In order to check if the full-state EKF can detect any other ramp faults, Figure 13 plots ${\alpha }_{EKF}$ as the magnitude of the fault slope in N2 and N6 is increased. Detection occurs when the slope is −0.5 m/s. Next, to see if detection occurs before the VPE crosses its alert limit or not, faults are injected into N2 and N6 from 895 s in a simulation run. It should be noted that the final approach of the aircraft starts at 911 s and ends around 975 s when the minimum decision altitude of 350 ft is reached. During this period, VPE can be compared against the alert limit of 50 m for EKF.

Table 3. Ramp fault detection performance of full-state EKF and SKF with different slopes. Time is noted from the beginning of simulation.

Fault Slope (m/s)	Full-State EKF		SKF
	Fault Detection Time (s)	Time When VPE > Alert Limit (s)	Fault Detection Time (s)	Time When VPE > VPL (s)
−0.5	1055	≥975	925	≥999
−0.75	937	≥948	914	≥963
−1	904	≥934.8	909	≥946

shows the performance of EKF and SKF. With −0.5 m/s, ${\alpha }_{EKF}$ detects faults after the VPE crosses the alert limit, but with the other two slopes, detection occurs before exceeding the alert limit, as desired. Thus, full-state EKF can detect ramp faults successfully when the slope magnitude is greater than a minimum value, which in this case is more than 0.5 m/s. It is also interesting to note that the VPE of EKF grows faster than that of SKF in the presence of ramp faults. As a result, it crosses its alert limit earlier than its counterpart in SKF.

Test results of KF RAIM for a 54.48 min aircraft trajectory are also presented in the Supplementary Materials. The results demonstrate the following: (1) adaptive calculations of all KF RAIM parameters of Figure 9 for long duration data, and (2) the ability of the algorithm to respond to an anomalous behavior and recover from it when the source of the anomaly is gone. The low-end laptop needs 16.3 min to run the KF RAIM algorithm in Matlab. WLS and full-state EKF results are also included. Details of all figures are provided in the description of the Supplementary Materials.

7. Conclusions

This paper integrates realistic error models of the carrier-smoothed pseudorange measurements into the RAIM algorithm with an SKF navigation filter. Error models include appropriate time constants of individual error sources, thus overcoming an important limitation of the prior work on KF RAIM. Uncertainties of the model parameters are accounted for. A provision is also made for adding artificial white noise to deal with initial transients in the position error. The suitability of the adopted models for RAIM with SKF is justified through extensive simulations. Performance of the modified KF RAIM is then analyzed with simulated GPS L1 and NavIC L5 signals for aircraft navigation during different phases of flight. The following important conclusions are drawn from the analyses of the results. WLS RAIM has lower PLs than the developed KF RAIM, but it is not appropriate for advanced navigation systems such as vector tracking and sensor integration, where KF is an integral part. A full-state EKF mistakes slow ramp faults for the residual measurement errors it needs to estimate, which results in missed detection. An SKF implementation, on the other hand, can successfully detect such faults. Hence, it offers a promising solution to developing KF RAIM algorithms in the range domain. With realistic error models, the modified KF RAIM is shown to complete processing well within time. Some salient features are also studied to gain insights into its working principles. It is noted that the mean position error bound is mainly determined by a few recent epoch measurements. However, contributions of all epochs must be taken into account for rigorous calculation of the PL for high-reliability applications.

Future work aims at determining error model parameters specific to NavIC measurements, and subsequently validating the KF RAIM with real GPS and NavIC signals and different geometries. Constellation-wide faults need to be handled in the range-based formulation. Furthermore, a fault isolation methodology is required to successfully exclude faulty satellites upon detection. The performance of KF RAIM in weak signal environments will be studied. Its suitability for sensor integration applications will be explored. As ionospheric scintillation prevalent in low latitude regions such as the Indian subcontinent can degrade GNSS measurements and introduce non-Gaussian errors [42,43], it is also important to analyze KF RAIM performance in such scenarios.