RAIM and Failure Mode Slope: Effects of Increased Number of Measurements and Number of Faults

Abstract

This article provides a comprehensive analysis of the impact of the increasing number of measurements and the possible increase in the number of faults in multi-constellation Global Navigation Satellite System (GNSS) Receiver Autonomous Integrity Monitoring (RAIM). Residual-based fault detection and integrity monitoring techniques are ubiquitous in linear over-determined sensing systems. An important application is RAIM, as used in multi-constellation GNSS-based positioning. This is a field in which the number of measurements, m, available per epoch is rapidly increasing due to new satellite systems and modernization. Spoofing, multipath, and non-line of sight signals could potentially affect a large number of these signals. This article fully characterizes the impact of measurement faults on the estimation (i.e., position) error, the residual, and their ratio (i.e., the failure mode slope) by analyzing the range space of the measurement matrix and its orthogonal complement. For any fault scenario affecting h measurements, the eigenvalue problem that defines the worst-case fault is expressed and analyzed in terms of these orthogonal subspaces, which enables further analysis. For $h>(m-n)$, where n is the number of estimated variables, it is known that there always exist faults that are undetectable from the residual vector, yielding an infinite value for the failure mode slope. This article uses the range space and its complement to explain: (1) why, for fixed h and n, the failure mode slope decreases with m; (2) why, for a fixed n and m, the failure mode slope increases toward infinity as h increases; (3) why a failure mode slope can become infinite for $h\le (m-n)$. A set of examples demonstrate the results of the paper.

1. Introduction

Multi-constellation Global Navigation Satellite Systems (GNSSs) have reached a maturity level where they are routinely relied upon as critical positioning systems that require a high degree of integrity [1,2]. The integrity of a positioning system is a measure of how much trust can be allocated to the correctness of the position solution [3]. This measure includes the ability of that system to be monitored (by itself or otherwise) and to provide timely warnings to the user when the system should not be used for positioning. Integrity monitoring, thus, requires the system to detect faulty measurements before they cause out-of-specification performance [4].

Receiver Autonomous Integrity Monitoring (RAIM) is one well-established method of ensuring the consistency of the positioning solution and assessing the integrity risk posed by available measurements [4,5,6]. Integrity risk is defined as the probability of undetected faults causing unacceptably large positioning errors [7]. Current implementations of residual-based RAIM (RB RAIM) are designed to detect faulty measurements and to evaluate the safety risk posed by possibly undetected faults (see Section 23.7 in [8]).

Integrity monitoring is applicable to a wide array of sensors that provide redundant measurements. Early approaches to RAIM emphasized its use to supplement various modes of aerospace navigation [9,10,11] or the ability of RAIM to enhance satellite-based augmentation services [12], such as the Wide Area Augmentation System (WAAS) [13]. The research in [13,14] proposed an implementation of weighted RAIM to compute Vertical Protection Levels (VPLs) for precision approach and landing supported by information from WAAS. Notable is that all the above methods assumed only single-measurement faults. This was justified by the assumption that measurement faults were caused by partial or complete satellite failure and that the likelihood of more than one satellite failing simultaneously is small.

Multi-measurement fault scenarios can be computationally challenging. In a scenario where the total number of measurements is denoted by m and the number of faulty measurements is denoted by h, the number of combinations of sensor faults that must be considered is $\left(\genfrac{}{}{0pt}{}{m}{h}\right)$, which grows rapidly with both m and h.

The growing availability of multiple GNSS constellations with multi-frequency measurements available per system means that the number of measurements available per epoch is rapidly increasing; therefore, m entering the 50–100 range is not unreasonable even in low-cost receivers [15,16,17,18].

In emerging applications, such as urban automotive navigation, healthy satellites could have faulty measurements (i.e., out-of-specification range errors) due to non-line-of-sight signals [19,20], strong multipath [21,22], ionospheric effects [23,24,25], and other environmental factors [26]. In low-latitude regions, ionospheric scintillation can be more pronounced [27], and persistent satellite oscillator anomalies resembling ionospheric scintillation have been recently reported to affect multiple GPS satellites [28,29]. Spoofing is an important threat to GNSS, which is achieved by simultaneously altering all satellite measurements in a targeted region [30,31,32,33]. These examples demonstrate that multi-measurement faults (i.e., $h\gg 1$) are a reasonable possibility that cannot be ignored.

In such multi-fault hypotheses, the fault direction is unknown; therefore, an upper bound on the integrity risk is derived for the worst-case fault direction in each h-fault scenario, using the worst-case failure mode slope [7,34]. This slope is defined as the largest ratio of the fault-induced estimation error to the fault-induced residual [34]. The failure mode slope is infinite when a fault is undetectable from the residual vector. For example, when the fault is in the range space of the measurement matrix, it is undetectable and the failure mode slope is infinite (see p. 110 in [10]). Such a fault can always be found when $h>(m-n)$, where n is the dimension of the vector that is being estimated. Single (i.e., $h=1$) and double measurement (i.e., $h=2$) fault scenarios are considered in, e.g., refs. [13,34], respectively. Multi-measurement fault scenarios ($1\le h$) are considered in [7,35], where the worst-case fault vector was determined by solving an eigenvalue problem. The worst-case slope is the maximum eigenvalue and the worst-case fault direction corresponds to its eigenvector. The analysis of [7] considers only one component of the state vector, which would be important in applications such as the vertical error in precision aircraft landing.

However, in other applications, such as land and sea navigation, there is interest in how multiple faults affect vector quantities, i.e., the horizontal position [36,37].

This paper studies how the RB RAIM failure mode slope changes as a function of the number of measurements, m, and the number of faulty measurements, h. The analysis uses the SVD as a tool: (1) to explicitly characterize the fault-induced estimation error and the fault-induced residual; (2) to find and characterize the linear space of faults that affect the residual, but not the estimate; (3) to find and characterize the linear space of (undetectable) faults that affect the estimate, but not the residual; and (4) to provide an analytic expression for the worst-case fault direction for all $1\le h\le m$. In addition, the expression that results from the SVD approach provides additional physical insight into the solution. In the case where $h>(m-n)$, this article provides an orthogonal basis for the $(m-n)$ dimensional linear space of faults that have infinite failure mode slope and are, hence, undetectable by residual-based detectors. Finally, this paper shows that, even in the case that $h\le (m-n)$, it is possible to have an infinite failure mode slope and provides analytic expressions that determine when this is the case. The examples in Section 10.4 demonstrate the main points of the article, including an example of undetectable faults when $h=m-n$.

This paper is organized as follows: Section 2 introduces the measurement model, including assumptions and models of the noise and fault. Section 3 defines the notation and symbols adopted. Section 4 uses the SVD decomposition of the measurement matrix to analyze the effects of noise and faults on the estimation error and residual. Section 5 reviews concepts related to integrity risk, while Section 6 reviews the concept of the failure mode slope. Section 7 presents new results on the best and worst-case fault direction vectors for general fault scenarios, showing them to be intrinsically related to certain singular vectors of the measurement matrix. An analytical expression for the upper bound on the estimation error is also provided. Section 8 analyzes worst-case fault modes for different fault scenarios, beginning with single-measurement faults, then double-measurement faults, and, ultimately, multi-measurement faults. Section 9 discusses the generalization of the approach when only subspaces of the estimated state are of interest. Section 10 presents and discusses GNSS examples demonstrating the results herein. Section 11 provides concluding remarks.

2. Problem Formulation

Consider the following over-determined linear measurement model, which applies to multi-constellation GNSS:

$$\mathbf{y}=\mathbf{H}\mathit{x}+\mathsf{\eta }+\mathbf{f},$$

(1)

where $\mathbf{y}\in {\mathbb{R}}^{m\times 1}$ is the measurement vector; $\mathbf{H}\in {\mathbb{R}}^{m\times n}$ is the measurement matrix, with $m>n$ and $\mathrm{rank}\left(\mathbf{H}\right)=n$; $\mathit{x}\in {\mathbb{R}}^n$ is a vector to be estimated; $\mathsf{\eta }\in {\mathbb{R}}^m$ is measurement noise that is assumed to be Gaussian with $\mathsf{\eta }\sim \mathcal{N}(\mathbf{0},{\mathbf{C}}_{\eta })$, where ${\mathbf{C}}_{\eta }={\sigma }_{\eta }^2\mathbf{I}$; and $\mathbf{f}\in {\mathbb{R}}^m$ is the fault (sometimes also referred to as the bias or failure). The fault $\mathbf{f}$ is defined as:

$$\mathbf{f}\triangleq \mu \overrightarrow{\mathbf{f}},$$

(2)

where $\mu \in \mathbb{R}$ is the magnitude of the fault, and the unit vector $\overrightarrow{\mathbf{f}}$ is the fault direction. In the absence of a fault, $\mu =0$.

This papers investigates the impact of $\mathbf{f}$ on the estimation error and the fault detector, within the context of evaluating Hazardously Misleading Information (HMI). This paper uses the term fault scenario to mean a combination of faults affecting h measurements (i.e., the vector $\overrightarrow{\mathbf{f}}$ has h non-zero components). The failure mode is defined to be a specific combination of measurements within a fault scenario. Note that there are $M_h=\left(\genfrac{}{}{0pt}{}{m}{h}\right)$ different failure modes within each fault scenario. For each failure mode, one of the goals of this paper is to explicitly define the worst-case fault direction and its failure-mode slope.

In many applications, the analyst is only interested in a rotated subvector of $\mathit{x}$. Such situations are easily addressed within the context of this article by defining $\mathbf{z}=\mathbf{M}\mathit{x}$ and applying the methods herein to the vector $\mathbf{z}$ instead of $\mathit{x}$. In this context, $\mathbf{M}$ combines a rotation and projection matrix. For example, in an automotive application, the horizontal position error is usually of primary interest. Assume that GNSS provides an estimate of $\mathit{x}={[x,y,z,b]}^{\top }\in {\mathbb{R}}^4$, where $x,y,z$ are the position coordinates in an Earth-Centered Earth-Fixed (ECEF) reference frame, and b is the receiver clock bias. Then, $\mathbf{M}\in {\mathbb{R}}^{2\times 4}$ can be defined as:

$$\begin{array}{c}\hfill \mathbf{M}\triangleq \left[\begin{array}{c}{\mathbf{I}}_{2\times 2}|{\mathbf{0}}_{2\times 2}\end{array}\right]{}_E^T\mathbf{R},\end{array}$$

(3)

where ${}_E^T\mathbf{R}\in {\mathbb{R}}^{4\times 4}$ is a rotation matrix from ECEF to tangent plane [38]. This definition of $\mathbf{M}$ transforms $\mathit{x}$ from ECEF to tangent plane and only retains the horizontal position components, removing the vertical position and the clock bias components.

In the analysis that follows, the derivations will be performed on the general case, estimating the full estimation vector $\mathit{x}$. Later in the paper, the horizontal position special case will be considered.

3. Notation

Different equality symbols will be used to distinguish between definitions, computations, and theoretical models used for analysis. The symbol ‘≜’ indicates a definition of a specific model or quantity. The symbol ‘≐’ indicates that the quantity on the left-hand side can computed from the quantities on the right-hand side. The symbol ‘=’ indicates that an equation is a theoretical model, not a definition or a computation used in implementations. Such models are used for analysis and physical interpretation of quantities.

The analysis of this article will use the singular value decomposition (SVD). The SVD of a matrix $\mathbf{H}\in {\mathbb{R}}^{m\times n}$ is

$$\mathbf{H}\triangleq \left[{\mathbf{U}}_1,{\mathbf{U}}_2\right]\left[\begin{array}{c}\Sigma \\ \mathbf{0}\end{array}\right]{\mathbf{V}}^{\top },$$

(4)

where $\mathbf{U}=\left[{\mathbf{U}}_1,{\mathbf{U}}_2\right]\in {\mathbb{R}}^{m\times m}$ is a unitary matrix with two orthogonal submatrices ${\mathbf{U}}_1=[{\mathbf{u}}_1^1,\cdots ,{\mathbf{u}}_1^n]\in {\mathbb{R}}^{m\times n}$ and ${\mathbf{U}}_2=[{\mathbf{u}}_2^1,\cdots ,{\mathbf{u}}_2^{m-n}]\in {\mathbb{R}}^{m\times (m-n)}$; $\Sigma \in {\mathbb{R}}^{n\times n}$ is positive definite and diagonal, where the diagonal elements ${\alpha }_1\ge {\alpha }_2\ge \dots \ge {\alpha }_n>0$ are the singular values of $\mathbf{H}$; and $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ is unitary. These matrices have interpretations that are useful for the analysis: $\mathbf{V}$ is an orthonormal basis for the domain of $\mathbf{H}$; ${\mathbf{U}}_1$ is an orthonormal basis for the n dimensional linear space that is the range of $\mathbf{H}$; and ${\mathbf{U}}_2$ is an orthonormal basis for the $(m-n)$ dimensional linear space that is orthogonal to the range of $\mathbf{H}$.

The term eigenvalue problem will refer to the constrained quadratic optimization problem of the form:

$$\begin{array}{cc}\hfill {\mathit{x}}^{*}& =arg\underset{{\parallel \mathit{x}\parallel }^2=1}{max}{\mathit{x}}^{\top }\mathbf{A}\mathit{x}=arg\underset{\mathit{x}\ne \mathbf{0}}{max}\frac{{\mathit{x}}^{\top }\mathbf{A}\mathit{x}}{{\mathit{x}}^{\top }\mathit{x}},\hfill \end{array}$$

(5)

where $\mathbf{A}\in {\mathbb{R}}^{n\times n}$ is a real symmetric matrix. The solution ${\mathit{x}}^{*}\in {\mathbb{R}}^n$ is the eigenvector corresponding to the largest eigenvalue of $\mathbf{A}$ [39]. That eigenvalue is ${max}_{{\parallel \mathit{x}\parallel }^2=1}{\mathit{x}}^{\top }\mathbf{A}\mathit{x}$.

When necessary, the actual and estimated values of a quantity are distinguished by x and $\widehat{x}$, respectively. Vectors and matrices are written in boldface, with vectors in lowercase and matrices in uppercase. For example, $\mathbf{v}\in {\mathbb{R}}^n$ is an n-element vector and $\mathbf{V}\in {\mathbb{R}}^{n\times n}$ is an $n\times n$ matrix. Positive definiteness (i.e., ${\mathbf{v}}^{\top }\mathbf{V}\mathbf{v}>0,\forall \mathbf{v}\ne \mathbf{0}$) and semi-definiteness (i.e., ${\mathbf{v}}^{\top }\mathbf{V}\mathbf{v}\ge 0,\forall \mathbf{v}\ne \mathbf{0}$) of a matrix is indicated by $\mathbf{V}\succ 0$ and $\mathbf{V}⪰0$, respectively. The matrix ${\mathbf{I}}_{n\times n}={\mathbf{I}}_n\in {\mathbb{R}}^{n\times n}$ is the $n\times n$ identity matrix. The symbol $\overrightarrow{\mathbf{v}}$ indicates the unit vector in the direction of vector $\mathbf{v}$. The vector ${\mathbf{e}}_i$ is the i-th column of $\mathbf{I}$, which is the i-th standard basis vector of ${\mathbb{R}}^n$, where $i\le n$. Throughout the article, $\mathrm{tr}\left(·\right)$ is the trace operator and $\parallel ·\parallel$ is the $L^2$ norm. For random vector $\zeta$, its mean squared value will be denoted as ${\parallel \zeta \parallel }_M^2=E\langle {\zeta }^{\top }\zeta \rangle$. Other notations will be defined when used.

4. Analysis of the Estimation Error and Residual

For measurements described by Equation (1), the minimum-variance unbiased (MVU) estimate of $\mathit{x}$ is computed as [40]:

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\widehat{\mathit{x}}& \doteq {\left({\mathbf{H}}^{\top }{\mathbf{C}}_{\eta }^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }{\mathbf{C}}_{\eta }^{-1}\mathbf{y}\hfill \end{array}$$

(6)

$$\begin{array}{cc}& \doteq {\mathbf{H}}^{*}\mathbf{y},\hfill \end{array}$$

(7)

where ${\mathbf{H}}^{*}\triangleq {\left({\mathbf{H}}^{\top }{\mathbf{C}}_{\eta }^{-1}\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }{\mathbf{C}}_{\eta }^{-1}$. In the special case where ${\mathbf{C}}_{\eta }={\sigma }_{\eta }^2\mathbf{I}$, this simplifies to ${\mathbf{H}}^{*}\triangleq {\left({\mathbf{H}}^{\top }\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }$. This special case is the main focus of this article.

The analyses that follow use the SVD of $\mathbf{H}$ to analyze the impact of faults on the estimation of $\widehat{\mathit{x}}$. Substituting Equation (4) into the definition of ${\mathbf{H}}^{*}$ in Equation (7) yields ${\displaystyle {\mathbf{H}}^{*}=\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }}$, so that Equation (7) becomes:

$$\begin{array}{cc}\hfill \widehat{\mathit{x}}& \doteq \mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{y}.\hfill \end{array}$$

(8)

Equation (8) is derived in Equation (S5) of Section S1 in the Supplementary Materials [41].

4.1. Effects of Noise and Faults on the Estimate

To analyze the effect of the measurement noise $\eta$ and fault $\mathbf{f}$, substitute (1) into (7):

$$\begin{array}{cc}\hfill \widehat{\mathit{x}}& =\left({\left({\mathbf{H}}^{\top }\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }\right)(\mathbf{H}\mathit{x}+\eta +\mathbf{f})\hfill \\ \hfill & ={\left({\mathbf{H}}^{\top }\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }\mathbf{H}\mathit{x}+{\left({\mathbf{H}}^{\top }\mathbf{H}\right)}^{-1}{\mathbf{H}}^{\top }(\eta +\mathbf{f})\hfill \\ \hfill & =\mathit{x}+{\mathbf{H}}^{*}(\eta +\mathbf{f}).\hfill \end{array}$$

(9)

The estimate $\widehat{\mathit{x}}$ is a Gaussian random variable with

$$\mathrm{E}〈\widehat{\mathit{x}}〉=\mathit{x}+{\mathbf{H}}^{*}\mathbf{f} \mathrm{and} {\mathbf{C}}_{\mathit{x}}=\mathrm{cov}\left(\widehat{\mathit{x}}\right)={\sigma }_{\eta }^2{\left({\mathbf{H}}^{\top }\mathbf{H}\right)}^{-1}={\sigma }_{\eta }^2\mathbf{V}{\Sigma }^{-2}{\mathbf{V}}^{\top }.$$

The estimation error defined as $\delta \mathit{x}=\widehat{\mathit{x}}-\mathit{x}$ relates to the noise and fault vector as:

$$\begin{array}{cc}\hfill \delta \mathit{x}& ={\mathbf{H}}^{*}(\eta +\mathbf{f})=\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\eta +\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{f}\hfill \\ \hfill & =\delta {\mathit{x}}_{\eta }+\delta {\mathit{x}}_{\mathbf{f}}.\hfill \end{array}$$

(10)

with noise-induced error $\delta {\mathit{x}}_{\eta }\triangleq \mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\eta$ and fault-induced bias $\delta {\mathit{x}}_{\mathbf{f}}\triangleq \mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{f}$. The expected value and covariance of the estimation error are:

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} \mathrm{E}〈\delta \mathit{x}〉& ={\mathbf{H}}^{*}\mathrm{E}〈\eta +\mathbf{f}〉=\delta {\mathit{x}}_{\mathbf{f}} \mathrm{and}\hfill \end{array}$$

(11)

$$\begin{array}{cc}\hfill \mathrm{cov}\left(\delta \mathit{x}\right)& ={\sigma }_{\eta }^2\mathbf{V}{\Sigma }^{-2}{\mathbf{V}}^{\top }.\hfill \end{array}$$

(12)

The fault induces an unknown bias $\delta {\mathit{x}}_{\mathbf{f}}$ on the estimate but does not affect its covariance, which remains ${\mathbf{C}}_{\mathit{x}}$. The estimation error distribution is Gaussian:

$$\delta \mathit{x}\sim \mathcal{N}(\delta {\mathit{x}}_{\mathbf{f}},{\mathbf{C}}_{\mathit{x}}).$$

(13)

Section S2 of the Supplementary Materials [41] shows that the mean-squared-error (MSE) of the estimate is:

$$\begin{array}{cc}\hfill {\parallel \delta \mathit{x}\parallel }_M^2& \triangleq \mathrm{E}〈{\left(\delta \mathit{x}\right)}^{\top }\left(\delta \mathit{x}\right)〉={\parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel }^2+{\kappa }_1,\hfill \end{array}$$

(14)

where

$$\begin{array}{cc}\hfill \parallel \delta {\mathit{x}}_{\mathbf{f}}{\parallel }^2={\parallel {\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{f}\parallel }^2 \mathrm{and} {\kappa }_1& \triangleq \mathrm{tr}\left({\mathbf{C}}_{\mathit{x}}\right).\hfill \end{array}$$

(15)

Equation (14) shows that the fault and noise have independent effects on the norm of the estimation error.

The first term allows the evaluation of the effect of specific faults on the estimate, given any specific instances of the measurement matrix $\mathbf{H}$ and fault vector $\mathbf{f}$. This will allow the analyst to find bounds on the estimation error under different fault scenarios.

It is convenient to rewrite the fault $\mathbf{f}$ using the orthogonal basis of ${\mathbb{R}}^m$ defined by the columns of $\mathbf{U}$:

$$\begin{array}{c}\hfill \mathbf{f}=\sum _{k=1}^n{a}_k{\mathbf{u}}_1^k+\sum _{\ell =1}^{m-n}{b}_{\ell }{\mathbf{u}}_2^{\ell }={\mathbf{U}}_1\mathbf{a}+{\mathbf{U}}_2\mathbf{b}=\left[\begin{array}{cc}{\mathbf{U}}_1& {\mathbf{U}}_2\end{array}\right]\left[\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right],\end{array}$$

(16)

where $\mathbf{a}\in {\mathbb{R}}^n$ and $\mathbf{b}\in {\mathbb{R}}^{(m-n)}$ are coefficient vectors. Specifically,

$$\begin{array}{cc}\hfill {\left[\begin{array}{cc}\mathbf{a}& \mathbf{b}\end{array}\right]}^{\top }& ={\left[\begin{array}{cc}{\mathbf{U}}_1& {\mathbf{U}}_2\end{array}\right]}^{\top }\mathbf{f}.\hfill \end{array}$$

(17)

Substituting the orthogonal decomposition of (16) into (10) yields:

$$\begin{array}{cc}\hfill \delta \mathit{x}& =\left(\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\right)\left(\left[\begin{array}{cc}{\mathbf{U}}_1& {\mathbf{U}}_2\end{array}\right]\left[\begin{array}{c}\mathbf{a}\\ \mathbf{b}\end{array}\right]+\eta \right)=\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }{\mathbf{U}}_1\mathbf{a}+\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }{\mathbf{U}}_2\mathbf{b}+\delta {\mathit{x}}_{\eta }\hfill \\ \hfill & =\mathbf{V}{\Sigma }^{-1}\mathbf{a}+\delta {\mathit{x}}_{\eta }.\hfill \end{array}$$

(18)

Equation (18) shows that $\delta {\mathit{x}}_{\mathbf{f}}\triangleq \mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{f}=\mathbf{V}{\Sigma }^{-1}\mathbf{a}$, which has the following interpretation.

Fact 1. 

Any portion of a fault $\mathbf{f}$ that is within the span of the $(m-n)$ dimensional linear space defined by the columns of ${\mathbf{U}}_2$ has no effect on the estimate. △

4.2. Effects of Noise and Faults on the Measurement Residual

The measurement residual is

$$\mathbf{r}\doteq \mathbf{y}-\widehat{\mathbf{y}}=(\mathbf{I}-\mathbf{P})\mathbf{y},$$

(19)

where $\widehat{\mathbf{y}}\doteq \mathbf{H}\widehat{\mathit{x}}\doteq \mathbf{P}\mathbf{y}$ with $\mathbf{P}\triangleq \mathbf{H}{\mathbf{H}}^{*}={\mathbf{U}}_1{\mathbf{U}}_1^{\top }\in {\mathbb{R}}^{m\times m}$. The matrix $\mathbf{P}$ is the projection matrix onto the range space of $\mathbf{H}$. It is symmetric, positive semi-definite, and idempotent with $\mathrm{rank}\left(\mathbf{P}\right)=n$. See Lemma S1 of Section S1 in the Supplementary Materials [41].

The effect of noise $\eta$ and fault $\mathbf{f}$ on the measurement residual is analyzed by substituting (1) into (19):

$$\begin{array}{cc}\hfill \mathbf{r}& =(\mathbf{I}-\mathbf{P})(\mathbf{H}\mathit{x}+\eta +\mathbf{f})=(\mathbf{I}-\mathbf{P})\left(\eta +\mathbf{f}\right)\hfill \\ \hfill & =\mathbf{Q}\eta +\mathbf{Q}\mathbf{f}.\hfill \end{array}$$

(20)

Equation (20) shows that $\mathbf{r}={\mathbf{r}}_{\mathbf{f}}+{\mathbf{r}}_{\eta }$ where the fault-induced residual is ${\mathbf{r}}_{\mathbf{f}}=\mathbf{Q}\mathbf{f}$ and the noise-induced residual is ${\mathbf{r}}_{\eta }=\mathbf{Q}\eta .$ The matrix $\mathbf{Q}\triangleq (\mathbf{I}-\mathbf{P})={\mathbf{U}}_2{\mathbf{U}}_2^{\top }\in {\mathbb{R}}^{m\times m}$ is the projection matrix onto the complement of the range space of $\mathbf{H}$. Therefore, $\mathbf{Q}\mathbf{P}=\mathbf{P}\mathbf{Q}=\mathbf{0}$ and $\mathbf{Q}\mathbf{H}=\mathbf{0}$. The properties of $\mathbf{Q}$ are discussed in Lemma S2 of Section S1 of the Supplementary Materials [41]. It is a symmetric, positive semi-definite, and idempotent matrix with $\mathrm{rank}\left(\mathbf{Q}\right)=(m-n)$. Substituting Equation (16) into (20) yields

$$\begin{array}{cc}\hfill \mathbf{r}& =\mathbf{Q}\eta +\mathbf{Q}{\mathbf{U}}_2\mathbf{b}\hfill \\ \hfill & =\mathbf{Q}\eta +{\mathbf{U}}_2\mathbf{b}.\hfill \end{array}$$

(21)

Equations (20) and (21) have the following interpretation.

Fact 2. 

Any portion of a fault $\mathbf{f}$ that is within the span of the n dimensional linear space defined by the columns of ${\mathbf{U}}_1$ has no effect on the residual. △

The expected value and covariance of $\mathbf{r}$ are:

$$\begin{array}{cc}\hfill \mathrm{E}\langle \mathbf{r}\rangle & =\mathbf{Q}\mathbf{f}={\mathbf{U}}_2{\mathbf{U}}_2^{\top }\mathbf{f}\mathrm{and}{\mathbf{C}}_{\mathbf{r}}\doteq {\sigma }_{\eta }^2\mathbf{Q}={\sigma }_{\eta }^2{\mathbf{U}}_2{\mathbf{U}}_2^{\top }.\hfill \end{array}$$

(22)

The covariance matrix ${\mathbf{C}}_{\mathbf{r}}$ is only positive semi-definite, so not invertible [42]. The distribution of $\mathbf{r}$ is Gaussian:

$$\begin{array}{c}\hfill \mathbf{r}\sim \mathcal{N}(\mathbf{Q}\mathbf{f},{\mathbf{C}}_{\mathbf{r}}).\end{array}$$

(23)

Section S3 of the Supplementary Materials [41] shows that the MSE of the residual is:

$$\begin{array}{c}\hfill {\parallel \mathbf{r}\parallel }_M^2={\parallel {\mathbf{r}}_{\mathbf{f}}\parallel }^2+{\kappa }_2,\end{array}$$

(24)

where

$$\begin{array}{c}\hfill \parallel {\mathbf{r}}_{\mathbf{f}}{\parallel }^2=\parallel {\mathbf{U}}_2^{\top }{\mathbf{f}\parallel }^2\mathrm{and}{\kappa }_2\triangleq {\parallel {\mathbf{r}}_{\eta }\parallel }_M^2=(m-n){\sigma }_{\eta }^2=\mathrm{tr}\left({\mathbf{C}}_r\right).\end{array}$$

(25)

The residual MSE defined in (24) is the sum of a fault-dependent component, $\parallel {\mathbf{r}}_{\mathbf{f}}{\parallel }^2$, and a noise-dependent component, $\parallel {\mathbf{r}}_{\eta }{\parallel }_M^2$.

4.3. Fault Decisions

The residual-based test statistic is

$$T_r\doteq {\mathbf{r}}^{\top }\mathbf{r}={\parallel {\mathbf{r}}^{\top }\mathbf{r}\parallel }^2,$$

(26)

with faults declared when $T_r\ge T^{*}$. (Alternatively, RAIM methods can be defined based on a parity vector [7,10]. By the definition of a parity vector, it is the case that any parity vector must satisfy $\mathbf{p}\triangleq {\mathbf{R}}^{\top }{\mathbf{U}}_2^{\top }\mathbf{y}$, where $\mathbf{R}$ is a unitary matrix. Therefore, $\mathbf{r}={\mathbf{U}}_2\mathbf{R}\mathbf{p}$ and $\parallel \mathbf{r}\parallel =\parallel \mathbf{p}\parallel$ because ${\mathbf{U}}_2^{\top }{\mathbf{U}}_2=\mathbf{I}$. Because $\parallel \mathbf{r}\parallel =\parallel \mathbf{p}\parallel$, the two test statistics are completely equivalent. This article chooses to focus on the residual vector. All results herein extend directly to the parity vector.)

Based on Equation (19), the test statistic can also be computed as $T_r\doteq {\mathbf{y}}^{\top }\mathbf{Q}\mathbf{y}$. The procedures to select $T^{*}$ based on the probability of false alarm and continuity risk specifications are described in [13,43].

5. Integrity Risk Evaluation

Integrity risk evaluates the probability of having Hazardously Misleading Information (HMI) [8]. HMI refers to the case when the estimation error exceeds a predefined threshold referred to as alarm limit L, while the fault detector $T_r$ does not detect a fault. It is written as:

$$\begin{array}{c}\hfill HMI\triangleq \left(\parallel \delta \mathit{x}\parallel >L\right)\cap \left(T_r<T^{🟉}\right).\end{array}$$

(27)

The probability of HMI, $\mathrm{P}\left(HMI\right)$, is evaluated under different fault hypotheses $H_{h,i}$. The number of fault hypotheses for m measurements, of which h are fault-affected, is $M_h=\left(\genfrac{}{}{0pt}{}{m}{h}\right)$ (i.e., the number of permutations of m items taken h at a time). Therefore, the integrity risk can then be expressed as:

$$\begin{array}{cc}\hfill \mathrm{P}\left(HMI\right)& =\sum _{h=1}^m{\mathrm{P}}_h\left(HMI\right),\mathrm{where}\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill \hspace{1em}\hspace{1em} {\mathrm{P}}_h\left(HMI\right)& =\sum _{i=0}^{M_h}\mathrm{P}\left(HM{I}_{h,i}\right|H_{h,i})\mathrm{P}\left(H_{h,i}\right),\hfill \end{array}$$

(29)

where the possible $M_h$ hypotheses for each $h=1,\dots ,m$ are assumed mutually exclusive and jointly exhaustive.

Full evaluation of integrity risk requires evaluating each term in the summation of (28). However, this incurs a heavy computational cost when $M_h$ becomes large; therefore, some references employ hypothesis reduction approaches that only evaluate a subset of hypotheses, while accounting for the remaining hypotheses by assigning them a probability of occurrence [44]. For example, if a bound $P_R$ can be defined such that

$$P_R\ge \sum _{h=h^{🟉}+1}^m{\mathrm{P}}_h\left(HMI\right),$$

then,

$$\begin{array}{cc}\hfill \mathrm{P}\left(HMI\right)& \le \sum _{h=1}^{h^{🟉}}{\mathrm{P}}_h\left(HMI\right)+P_R.\hfill \end{array}$$

(30)

In this case, the analyst only needs to evaluate ${\mathrm{P}}_h\left(HMI\right)$ for $h=1,\dots ,h^{🟉}$. The procedures to determine $h^{🟉}$ and the corresponding $P_R$ are described in [45,46]. These papers presents reasonable values for $h^{🟉}$ in the context of the historical data of single satellite failure. They do not consider situations with large numbers of faulty measurements due to spoofing, non-line-of-sight signals, etc.

5.1. Hypothesis Probabilities

Let $\mathrm{P}\left(H_h\right)$ denote the probability of one of the h-fault scenarios occurring. Usually the probabilities of occurrence of all fault hypotheses $H_{h,i}$ are assumed to be equal: $\mathrm{P}\left(H_{h,i}\right)=\mathrm{P}\left(H_{h,j}\right)for1\le i,j\le M_h.$ Then, the probability of occurrence of the h-fault scenario is

$$\mathrm{P}\left(H_h\right)=\sum _{i=1}^{M_h}\mathrm{P}\left(H_{h,i}\right)=M_h\mathrm{P}\left(H_{h,1}\right)\mathrm{for} h=1,\cdots ,m.$$

The allocation of the nominal probabilities for each fault scenario is such that

$$\begin{array}{c}\hfill \sum _{h=1}^m{M}_h\mathrm{P}\left(H_{h,1}\right)=1-P_0,\end{array}$$

(31)

where $P_0$ is the probability of the fault-free scenario. The value of each $\mathrm{P}\left(H_{h,1}\right)$ is small and determined empirically. See, e.g., ref. [7].

5.2. Evaluating $\mathrm{P}\left(HM{I}_{h,i}\right|H_{h,i})$

The probability of HMI under the i-th of the h-fault hypotheses is

$$\begin{array}{cc}\hfill \mathrm{P}\left(HM{I}_{h,i}\right|H_{h,i})& =\mathrm{P}\left(\parallel \delta \widehat{\mathit{x}}{\parallel }^2>L^2,T_r<T^{🟉}|H_{h,i}\right)\hfill \\ \hfill & =\mathrm{P}\left((\parallel \delta \widehat{\mathit{x}}{\parallel }^2>L^2)|(T_r<T^{🟉},H_{h,i})\right)P\left(T_r<T^{🟉}|H_{h,i}\right).\hfill \end{array}$$

The failure mode slope, which is reviewed in Section 6, provides a means of predicting the probability that $\mathrm{P}\left((\parallel \delta \widehat{\mathit{x}}{\parallel }^2>L^2)|(T_r<T^{🟉},H_{h,i})\right)$ based on $T_r$, as defined in Equation (26).

6. Failure Mode Slope

The literature defines the failure mode slope $g_f$ as the ratio between the norm of the fault-induced bias of the estimate, $\parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel$ and the norm of fault-induced residual $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$:

$$\begin{array}{c}\hfill g_f\triangleq \frac{\parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel }{\parallel {\mathbf{r}}_{\mathbf{f}}\parallel },\end{array}$$

(32)

as in refs. [4,7,35,43,47,48,49,50]. Equation (32) can be reorganized as:

$$\begin{array}{c}\hfill \parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel =g_f\parallel {\mathbf{r}}_{\mathbf{f}}\parallel ,\end{array}$$

(33)

which shows that, for a given $g_f$, the norm of the fault-induced estimation bias grows proportionately with norm of the fault-induced residual. The largest estimation error $\parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel$ induced by the fault $\mathbf{f}$ that is not detectable for the residual test $T_r\le T^{🟉}$ is:

$$\begin{array}{c}\hfill \parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel =g_f\sqrt{T^{🟉}}.\end{array}$$

(34)

When the maximum failure mode slope is bounded, $g_f\le {\overline{g}}_f$, over a specified set of possible fault scenarios (e.g., $1\le h\le h^{🟉}$), then selecting $T^{🟉}$ as a function of $\frac{L^2}{{\overline{g}}_f^2}$ allows the RAIM designer to control the probability of HMI. If $T_r>T^{🟉}$ or $g_f>{\overline{g}}_f$ for the current set of available satellites, RAIM can declare the navigation solution to be unavailable [11,48].

Based on Equations (15) and (25), for any given fault direction $\overrightarrow{\mathbf{f}}$, the failure mode slope can be computed as:

$$\begin{array}{c}\hfill g_f=\frac{\parallel {\Sigma }^{-1}{\mathbf{U}}_1^{\top }\overrightarrow{\mathbf{f}}{\parallel }^2}{\parallel {\mathbf{U}}_2^{\top }\overrightarrow{\mathbf{f}}{\parallel }^2}=\frac{\parallel {\Sigma }^{-1}{\mathbf{a}\parallel }^2}{{\parallel \mathbf{b}\parallel }^2},\end{array}$$

(35)

where the fault direction $\overrightarrow{\mathbf{f}}$ defines the coordinate vectors $\mathbf{a}$ and $\mathbf{b}$, as in Equation (17). The failure mode slope depends on the satellite constellation geometry (i.e., $\mathbf{H}$ through $\Sigma$, ${\mathbf{U}}_1$ and ${\mathbf{U}}_2$), and on fault direction through $\overrightarrow{\mathbf{f}}$, but not on the fault magnitude.

7. Best and Worst-Case Faults

Equation (35) shows that it is the direction of the fault, not the magnitude, that is of interest. Because $\Sigma$ is diagonal, Equation (35) is equivalent to

$$\begin{array}{cc}\hfill g_f& \doteq \frac{{\sum }_{i=1}^n{\left(\frac{a_i}{{\alpha }_i}\right)}^2}{{\sum }_{j=1}^{m-n}{\left(b_j\right)}^2}.\hfill \end{array}$$

(36)

where ${\alpha }_i$ are the singular values of $\mathbf{H}$. Based on Equation (36), with $\mathbf{a}$ and $\mathbf{b}$ defined as in Equation (17), the following conclusions are straightforward:

Best Case: 

For any fault that has $\parallel \mathbf{a}\parallel =0$ and $\parallel \mathbf{b}\parallel =1$, the fault direction $\overrightarrow{\mathbf{f}}\in span\left({\mathbf{U}}_2\right)$. In this case, the numerator is zero and the failure mode slope $g_f=0.$ Physically, this means that the fault has absolutely no impact on the state estimate.

Worst Case: 

For any fault that has $\parallel \mathbf{a}\parallel =1$ and $\parallel \mathbf{b}\parallel =0$, the fault direction $\overrightarrow{\mathbf{f}}\in span\left({\mathbf{U}}_1\right)$. In this case, the denominator is zero and the failure mode slope $g_f=\infty .$ Physically this means that the fault has no impact on the residual. Therefore, the residual test cannot detect it.

Note that, in both cases, the fault direction is not unique. Rather, each case allows an uncountably infinite number of fault directions within a finite dimensional subspace of ${\mathbb{R}}^m$.

Within the worst-case subspace, certain directions can still be considered worse than others. The problem of finding a unit vector $\mathbf{a}$ that maximizes the numerator (i.e., the estimation error for $\mu =1$) is formulated as:

$$\begin{array}{cc}\hfill \mathbf{a}& =arg\underset{a}{max}\sum _{i=1}^n{\left(\frac{a_i}{{\alpha }_i}\right)}^2,\mathrm{subject}\mathrm{to}\sum _{i=1}^n{a}_i^2=1\hfill \\ \hfill & =arg\underset{{\parallel \mathbf{a}\parallel }^2=1}{max}\left({\mathbf{a}}^{\top }{\Sigma }^{-2}\mathbf{a}\right).\hfill \end{array}$$

(37)

The problem in Equation (37) is a constrained quadratic maximization problem as in (5), of which the solution is the eigenvector corresponding to the largest eigenvalue of ${\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }$.

Therefore, the worst-case fault direction that has the greatest impact on the estimate per unit change in fault magnitude is the vector ${\overrightarrow{\mathbf{f}}}_w={\mathbf{u}}_1^n$, making the worst-case fault vector $\mathbf{f}=\mu {\mathbf{u}}_1^n$. Using Equation (17), the coefficient vector corresponding to this fault direction is:

$${\mathbf{a}}_w=\mu {\mathbf{U}}_1^{\top }{\mathbf{u}}_1^n=\mu {\mathbf{e}}_n,$$

(38)

where ${\mathbf{e}}_n\in {\mathbb{R}}^n$ is the n-th standard n dimensional unit direction vector. Using Equation (14), the corresponding mean-squared estimation error is:

$$\begin{array}{cc}\hfill {\parallel \delta \mathit{x}\parallel }_M^2& ={\mu }^2\left({\mathbf{e}}_n^{\top }{\Sigma }^{-2}{\mathbf{e}}_n\right)+{\kappa }_1={\left(\frac{\mu }{{\alpha }_n}\right)}^2+{\kappa }_1.\hfill \end{array}$$

(39)

This shows that the growth in the estimation error caused by the fault $\mathbf{f}=\mu {\mathbf{u}}_1^n$ is directly proportional to the inverse of the smallest singular value of $\mathbf{H}$ and the residual vector is not affected.

This section discussed the best and worst-case fault directions in the general case. Typically, these worst-case faults have non-zero contributions from all sensors; therefore, they do not correspond to any specific fault scenario with $h<m$. When a fault has a non-zero component in the linear space spanned by ${\mathbf{U}}_2$ (i.e., $\parallel \mathbf{b}\parallel \ne 0$), then $g_f$ is finite. The next section analyzes and determines the worst-case fault for general h-fault scenarios.

8. Single, Double, Multi-Measurement Faults

This section considers different h-fault scenarios. The failure mode slope of Equation (32) can be used as the metric for fault impact; therefore, the worst-case fault direction is defined as:

$$\begin{array}{cc}\hfill {\overrightarrow{\mathbf{f}}}_w& =arg\underset{{\parallel \mathbf{f}\parallel }^2=1}{max}{g}_f\left(\mathbf{f}\right)=arg\underset{{\parallel \mathbf{f}\parallel }^2=1}{max}\frac{\parallel {\Sigma }^{-1}{\mathbf{U}}_1^{\top }\overrightarrow{\mathbf{f}}{\parallel }^2}{\parallel {\mathbf{U}}_2^{\top }\overrightarrow{\mathbf{f}}{\parallel }^2}.\hfill \end{array}$$

(40)

This problem will be solved for different h-fault scenarios, beginning with faults affecting single measurements (i.e., $h=1$), then two measurements (i.e., $h=2$), and finally culminating with a general discussion of multi-measurements faults (i.e., $h\ge 1$). For each scenario, the analysis determines conditions when there are faults that cause the failure mode slope to be infinite. When the failure mode slope is finite the worst-case fault direction is determined.

8.1. Single-Measurement Faults

Single-measurement faults have the form $\mathbf{f}=\mu {\mathbf{e}}_j$ for $j=1,\cdots ,m$. For each such single-measurement fault, Equation (40) becomes

$$\begin{array}{cc}\hfill g_f\left({\mathbf{e}}_j\right)& =\frac{{\mathbf{e}}_j^{\top }{\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }{\mathbf{e}}_j}{{\mathbf{e}}_j^{\top }{\mathbf{U}}_2{\mathbf{U}}_2^{\top }{\mathbf{e}}_j}=\frac{{\sum }_{i=1}^n{\left(\frac{{\mathbf{u}}_1^{i\left(j\right)}}{{\alpha }_i}\right)}^2}{{\sum }_{i=1}^{m-n}{\left({\mathbf{u}}_2^{i\left(j\right)}\right)}^2}.\hfill \end{array}$$

(41)

where ${\mathbf{u}}_1^{i\left(j\right)}$ is the j-th element of vector ${\mathbf{u}}_1^i$ for $i=1,\dots ,n$ and ${\mathbf{u}}_2^{i\left(j\right)}$ is the j-th element of vector ${\mathbf{u}}_2^i$ for $i=1,\dots ,(m-n)$. For each j, the value of $g_f\left({\mathbf{e}}_j\right)$ can be determined. The one with the largest value determines the worst-case single sensor failure (i.e., fault direction).

In the single fault scenario, $g_f$ can be infinite if an entire row of ${\mathbf{U}}_2$ is zero. For example, the rank 4 matrix $\mathbf{H}$ defined as:

$$\mathbf{H}=\frac{1}{2}\left[\begin{array}{cccc}1& -1& \sqrt{2}& 2\\ \sqrt{2}& \sqrt{2}& 0& 2\\ -1& 1& \sqrt{2}& 2\\ \sqrt{2}-1& \sqrt{2}+1& \sqrt{2}& 2\\ 1& -1& \sqrt{2}& 2\end{array}\right]\mathrm{has} {\mathbf{U}}_2=\left[\begin{array}{c}-\sqrt{2}\\ 0\\ 0\\ 0\\ \sqrt{2}\end{array}\right],$$

which means that a single-sensor fault on any one of measurements 2, 3, or 4 cannot be detected by the residual (or parity) test statistic.

8.2. Double-Measurement Faults

For a double-measurement fault, two of the measurements are affected, not necessarily to the same extent. The fault is a linear combination of the individual measurement directions:

$${\mathbf{f}}_2=\mu \left(s_j{\mathbf{e}}_j+s_k{\mathbf{e}}_k\right)=\mu {\mathbf{e}}_{jk},$$

(42)

where, for $j\ne k$, ${\mathbf{e}}_{jk}\triangleq \left(s_j{\mathbf{e}}_j+s_k{\mathbf{e}}_k\right)$ is the overall fault direction, and $s_j,s_k\in [-1,1]$ are scaling coefficients such that $\parallel {\mathbf{e}}_{jk}{\parallel }^2=s_j^2+s_k^2=1$.

The fault direction vector ${\mathbf{e}}_{jk}$ can be written as:

$$\begin{array}{c}\hfill {\mathbf{e}}_{jk}={\mathbf{D}}_{jk}\mathbf{s},\end{array}$$

(43)

where $\mathbf{s}\in {\mathbb{R}}^2$, and ${\mathbf{D}}_{jk}\in {\mathbb{R}}^{m\times 2}$ are a binary matrix of which the only non-zero elements are in the j-th and k-th positions of the two columns, respectively, i.e.:

$$\begin{array}{c}\hfill {\mathbf{D}}_{jk}=\left[\begin{array}{cc}|& |\\ {\mathbf{e}}_j& {\mathbf{e}}_k\\ |& |\end{array}\right].\end{array}$$

(44)

For each such two-fault mode, Equation (40) becomes

$$\begin{array}{cc}\hfill g_f\left({\mathbf{e}}_{jk}\right)& =\frac{{\mathbf{e}}_{jk}^{\top }{\mathbf{U}}_1{\Sigma }^{-1}{\mathbf{U}}_1^{\top }{\mathbf{e}}_{jk}}{{\mathbf{e}}_{jk}^{\top }{\mathbf{U}}_2{\mathbf{U}}_2^{\top }{\mathbf{e}}_{jk}}\mathrm{for}\parallel {\mathbf{e}}_{jk}\parallel \ne 0\hfill \\ \hfill & =\frac{{\mathbf{s}}^{\top }{\Gamma }_{jk}\mathbf{s}}{{\mathbf{s}}^{\top }{\Delta }_{jk}\mathbf{s}}\mathrm{for}\parallel \mathbf{s}\parallel \ne 0,\hfill \end{array}$$

(45)

with ${\Delta }_{jk}={\mathbf{D}}_{jk}^{\top }\mathbf{Q}{\mathbf{D}}_{jk}\in {\mathbb{R}}^{2\times 2}$, $\mathbf{Q}={\mathbf{U}}_2{\mathbf{U}}_2^{\top }$, and ${\Gamma }_{jk}={\mathbf{D}}_{jk}^{\top }{\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }{\mathbf{D}}_{jk}\in {\mathbb{R}}^{2\times 2}.$ Given these definitions,

$$\begin{array}{ccc}\hfill {\Gamma }_{jk}& =\left[\begin{array}{cc}\sum _{i=1}^n\frac{{\left({\mathbf{u}}_1^{i\left(j\right)}\right)}^2}{{\alpha }_i^2}& \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i\left(j\right)}{\mathbf{u}}_1^{i\left(k\right)}}{{\alpha }_i^2}\\ \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i\left(j\right)}{\mathbf{u}}_1^{i\left(k\right)}}{{\alpha }_i^2}& \sum _{i=1}^n\frac{{\left({\mathbf{u}}_1^{i\left(k\right)}\right)}^2}{{\alpha }_i^2}\end{array}\right]\mathrm{and}{\Delta }_{jk}\hfill & \hfill =\left[\begin{array}{cc}\sum _{i=1}^{m-n}{\left({\mathbf{u}}_2^{i\left(j\right)}\right)}^2& \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i\left(j\right)}{\mathbf{u}}_2^{i\left(k\right)}\\ \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i\left(j\right)}{\mathbf{u}}_2^{i\left(k\right)}& \sum _{i=1}^{m-n}{\left({\mathbf{u}}_2^{i\left(k\right)}\right)}^2\end{array}\right].\end{array}$$

Both ${\Gamma }_{jk}$ and ${\Delta }_{jk}$ are symmetric and, at least, positive semi-definite.

The failure mode slope for each pair $(j,k)$ is determined by the value of $\mathbf{s}$ that maximizes $g_f\left({\mathbf{e}}_{jk}\right)$, as defined in Equation (45). This optimization problem is:

$$\begin{array}{cc}\hfill {\mathbf{s}}_{jk}^{*}& =arg\underset{\mathbf{s}\ne \mathbf{0}}{max}\frac{{\mathbf{s}}^{\top }{\Gamma }_{jk}\mathbf{s}}{{\mathbf{s}}^{\top }{\Delta }_{jk}\mathbf{s}}.\hfill \end{array}$$

(46)

When ${\Delta }_{jk}$ is positive definite (i.e., nonsingular and invertible), then it has a square root matrix, denoted as ${\Delta }_{jk}^{\frac{1}{2}}$, which is, itself, positive definite and symmetric [42], and ${\Delta }_{jk}$ can be expressed as ${\Delta }_{jk}={\Delta }_{jk}^{\frac{1}{2}}{\Delta }_{jk}^{\frac{1}{2}}$. As discussed in [35], the problem in Equation (46) is equivalent to

$$\begin{array}{cc}\hfill {\mathbf{s}}_{jk}^{*}& =arg\underset{\mathbf{s}\ne \mathbf{0}}{max}\frac{{\mathbf{s}}^{\top }{\Delta }_{jk}^{\frac{1}{2}}{\Delta }_{jk}^{-\frac{1}{2}}{\Gamma }_{jk}{\Delta }_{jk}^{-\frac{1}{2}}{\Delta }_{jk}^{\frac{1}{2}}\mathbf{s}}{{\mathbf{s}}^{\top }{\Delta }_{jk}^{\frac{1}{2}}{\Delta }_{jk}^{\frac{1}{2}}\mathbf{s}}.\hfill \end{array}$$

(47)

Letting $\mathbf{q}={\Delta }_{jk}^{\frac{1}{2}}\mathbf{s}$, the problem becomes:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathbf{q}}_{jk}^{*}=arg\underset{\mathbf{q}\ne \mathbf{0}}{max}\frac{{\mathbf{q}}^{\top }{\Psi }_{jk}\mathbf{q}}{{\mathbf{q}}^{\top }\mathbf{q}},}\hfill \end{array}\end{array}$$

(48)

where ${\Psi }_{jk}={\Delta }_{jk}^{-\frac{1}{2}}{\Gamma }_{jk}{\Delta }_{jk}^{-\frac{1}{2}}\in {\mathbb{R}}^{2\times 2}$. The problem in (48) is an eigenvalue problem of which the solution ${\mathbf{q}}_{jk}^{*}$ is the eigenvector corresponding to the largest eigenvalue, denoted as ${\lambda }_{jk}^{*}$, of ${\Psi }_{jk}$. Therefore, ${\mathbf{s}}^{*}={\Delta }^{-\frac{1}{2}}{\mathbf{q}}^{*}$ and $g_f^{*}={\lambda }_{jk}^{*}$.

Each pair of satellites yields a value for ${\lambda }_{jk}^{*}$. There are $\left(\genfrac{}{}{0pt}{}{m}{2}\right)$ possible combinations that must be evaluated to quantify the worst-case impact of double-measurement faults on the estimation error and solution integrity.

8.3. Multi-Measurement Faults

This section generalizes the approach of Section 8.2 to the general case of h sensor failures.

A fault that affects h measurements, for $1\le h\le m$, can be written as:

$$\begin{array}{c}\hfill {\mathbf{f}}_h=\mu {\mathbf{e}}_{\pi },\end{array}$$

(49)

where $\pi$ denotes a permutation of the set of integers between 1 and m that has cardinality $\left|\pi \right|=h$. The fault direction vector ${\mathbf{e}}_{\pi }$ can be written as:

$$\begin{array}{c}\hfill {\mathbf{e}}_{\pi }={\mathbf{D}}_{\pi }\mathbf{s},\end{array}$$

(50)

where $\parallel {\mathbf{e}}_{\pi }{\parallel }^2=1$ and $\mathbf{s}\in {\mathbb{R}}^h$ with ${\parallel \mathbf{s}\parallel }^2=1$. The matrix ${\mathbf{D}}_{\pi }\in {\mathbb{R}}^{m\times h}$ is as defined in (44), but with h unique columns.

For the h-fault scenario, the $g_f\left(\mathbf{f}\right)$ in Equation (40) is equivalent to

$$\begin{array}{cc}\hfill g_f\left({\mathbf{e}}_{\pi }\right)& =\frac{{\mathbf{e}}_{\pi }^{\top }{\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }{\mathbf{e}}_{\pi }}{{\mathbf{e}}_{\pi }^{\top }{\mathbf{U}}_2{\mathbf{U}}_2^{\top }{\mathbf{e}}_{\pi }}for\parallel {\mathbf{e}}_{\pi }\parallel \ne 0\hfill \end{array}$$

(51)

$$\begin{array}{cc}\hfill & =\frac{{\mathbf{s}}^{\top }{\Gamma }_{\pi }\mathbf{s}}{{\mathbf{s}}^{\top }{\Delta }_{\pi }\mathbf{s}}for\parallel \mathbf{s}\parallel \ne 0.\hfill \end{array}$$

(52)

Letting ${\pi }_j$ denote the j-th element of $\pi$ for $j\in \{1,\dots ,h\}$ and defining $\Phi ={\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }$, the matrix ${\Gamma }_{\pi }={\mathbf{D}}_{\pi }^{\top }\Phi {\mathbf{D}}_{\pi }\in {\mathbb{R}}^{h\times h}$ has the form:

$$\begin{array}{c}\hfill {\Gamma }_{\pi }=\left[\begin{array}{cccc}\sum _{i=1}^n{\left(\frac{{\mathbf{u}}_1^{i,{\pi }_1}}{{\alpha }_i}\right)}^2& \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_1}{\mathbf{u}}_1^{i,{\pi }_2}}{{\alpha }_i^2}& \cdots & \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_1}{\mathbf{u}}_1^{i,{\pi }_h}}{{\alpha }_i^2}\\ \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_2}{\mathbf{u}}_1^{i,{\pi }_1}}{{\alpha }_i^2}& \sum _{i=1}^n{\left(\frac{{\mathbf{u}}_1^{i,{\pi }_2}}{{\alpha }_i}\right)}^2& \cdots & \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_2}{\mathbf{u}}_1^{i,{\pi }_h}}{{\alpha }_i^2}\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_h}{\mathbf{u}}_1^{i,{\pi }_1}}{{\alpha }_i^2}& \sum _{i=1}^n\frac{{\mathbf{u}}_1^{i,{\pi }_h}{\mathbf{u}}_1^{i,{\pi }_2}}{{\alpha }_i^2}& \cdots & \sum _{i=1}^n{\left(\frac{{\mathbf{u}}_1^{i,{\pi }_h}}{{\alpha }_i}\right)}^2\end{array}\right].\end{array}$$

where ${\mathbf{u}}_1^{i,{\pi }_j}$ denotes the element in row i and column ${\pi }_j$ of ${\mathbf{U}}_1$ (similarly for ${\mathbf{u}}_2^{i,{\pi }_j}$). The matrix ${\Delta }_{\pi }={\mathbf{D}}_{\pi }^{\top }\mathbf{Q}{\mathbf{D}}_{\pi }\in {\mathbb{R}}^{h\times h}$ has the form:

$$\begin{array}{c}\hfill {\Delta }_{\pi }=\left[\begin{array}{cccc}\sum _{i=1}^{m-n}{\left({\mathbf{u}}_2^{i,{\pi }_1}\right)}^2& \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_1}{\mathbf{u}}_2^{i,{\pi }_2}& \cdots & \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_1}{\mathbf{u}}_2^{i,{\pi }_h}\\ \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_2}{\mathbf{u}}_2^{i,{\pi }_1}& \sum _{i=1}^{m-n}{\left({\mathbf{u}}_2^{i,{\pi }_2}\right)}^2& \cdots & \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_2}{\mathbf{u}}_2^{i,{\pi }_h}\\ ⋮& ⋮& \ddots & ⋮\\ \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_h}{\mathbf{u}}_2^{i,{\pi }_1}& \sum _{i=1}^{m-n}{\mathbf{u}}_2^{i,{\pi }_h}{\mathbf{u}}_2^{i,{\pi }_2}& \cdots & \sum _{i=1}^{m-n}{\left({\mathbf{u}}_2^{i,{\pi }_h}\right)}^2\end{array}\right].\end{array}$$

The matrix ${\Delta }_{\pi }$ is singular if any vector in the span of ${\mathbf{D}}_{\pi }$ lies entirely in the span of ${\mathbf{U}}_1$ or if at least one of the rows of ${\mathbf{U}}_2$ is zero. Additionally, Section S4 of the Supplementary Materials [41] shows that ${\Delta }_{\pi }$ is singular whenever $h>m-n$ but does not say anything about the case $h\le (m-n)$.

As was conducted to convert from Equation (46) to (48), the square root of ${\Delta }_{\pi }$ defines ${\Psi }_{\pi }={\Delta }_{\pi }^{-\frac{1}{2}}{\Gamma }_{\pi }{\Delta }_{\pi }^{-\frac{1}{2}}\in {\mathbb{R}}^{h\times h}$, where $\mathbf{q}={\Delta }_{\pi }^{\frac{1}{2}}\mathbf{s}$. The worst-case fault direction problem from Equation (40) for the h-fault scenario is equivalent to:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle {\mathbf{q}}_{\pi }^{*}=arg\underset{\mathbf{q}\ne \mathbf{0}}{max}\frac{{\mathbf{q}}^{\top }{\Psi }_{\pi }\mathbf{q}}{{\mathbf{q}}^{\top }\mathbf{q}}.}\hfill \end{array}\end{array}$$

(53)

For some values of h, ${\Delta }_{\pi }$ can be singular for some $\pi$. Therefore, there are two cases to be considered: ${\Delta }_{\pi }$ is singular and ${\Delta }_{\pi }$ is nonsingular.

When ${\Delta }_{\pi }$ is nonsingular, the solution to Equation (53) is the eigenvector corresponding to the largest eigenvalue ${\lambda }_{\pi }^{*}$ of ${\Psi }_{\pi }$. Then, ${\mathbf{s}}^{*}={\Delta }_{\pi }^{-\frac{1}{2}}{\mathbf{q}}_{\pi }^{*}$ and the worst-case fault direction is:

$$\begin{array}{c}\hfill {\overrightarrow{\mathbf{f}}}_w={\mathbf{D}}_{\pi }{\Delta }_{\pi }^{-\frac{1}{2}}{\mathbf{q}}_{\pi }^{*}.\end{array}$$

(54)

The failure mode slope $g_f^{*}={\lambda }_{\pi }^{*}$ is finite.

When ${\Delta }_{\pi }$ is singular, the matrix ${\Delta }_{\pi }^{-\frac{1}{2}}$ does not exist, so the problem in Equation (53) cannot be formulated and solved. This case is discussed in Section 8.4.

8.4. Undetectable Faults

This section further discusses the case that ${\Delta }_{\pi }$ is singular. Let q denote the number of zero eigenvalues of ${\Delta }_{\pi }$, with $0<q<h$. Consider an eigendecomposition ${\Delta }_{\pi }=\mathbf{W}\mathbf{A}{\mathbf{W}}^{\top }$, where $\mathbf{W}=[{\mathbf{w}}_1,\cdots ,{\mathbf{w}}_q,\cdots ,{\mathbf{w}}_h]\in {\mathbb{R}}^{h\times h}$ and $\mathbf{A}=\mathrm{diag}({\lambda }_1,\cdots ,{\lambda }_h)\in {\mathbb{R}}^{h\times h}$, with ${\left\{{\lambda }_i\right\}}_{i=1,\cdots ,q}=0$. The eigenvectors corresponding to the zero eigenvalues span a linear subspace of which the basis is ${\left\{{\mathbf{w}}_i\right\}}_{i=1,\cdots ,q}$. For $\mathbf{s}$ in this linear subspace:

$$\begin{array}{c}\hfill \mathbf{s}=\sum _{i=1}^q{c}_i{\mathbf{w}}_i={\mathbf{W}}_q\mathbf{c},\end{array}$$

where ${\mathbf{W}}_q=[{\mathbf{w}}_1,\cdots ,{\mathbf{w}}_q]\in {\mathbb{R}}^{h\times q}$ and $\mathbf{c}\in {\mathbb{R}}^q$, with ${\parallel \mathbf{c}\parallel }^2=1$. Any $\mathbf{s}$ in this subspace will be in the null space of ${\Delta }_{\pi }$, meaning that ${\mathbf{r}}_{\mathbf{f}}=\mathbf{0}$. Therefore, $g_f=\infty$. The fault-induced estimation error is $\delta {\mathit{x}}_{\mathbf{f}}={\Sigma }^{-1}{\mathbf{U}}_1^{\top }{\mathbf{D}}_{\pi }{\mathbf{W}}_q\mathbf{c}$, with norm

$$\begin{array}{cc}\hfill \parallel \delta {\mathit{x}}_{\mathbf{f}}{\parallel }^2& ={\mathbf{c}}^{\top }{\rm Y}\mathbf{c},\hfill \end{array}$$

(55)

where ${\rm Y}={\mathbf{W}}_q^{\top }{\mathbf{D}}_{\pi }^{\top }{\mathbf{U}}_1{\Sigma }^{-2}{\mathbf{U}}_1^{\top }{\mathbf{D}}_{\pi }{\mathbf{W}}_q$ is symmetric and positive semi-definite. The vector that maximizes the quadratic form in Equation (55) is:

$$\begin{array}{c}\hfill {\mathbf{c}}^{*}=arg\underset{{\parallel \mathbf{c}\parallel }^2=1}{max}{\mathbf{c}}^{\top }{\rm Y}\mathbf{c}.\end{array}$$

(56)

Equation (56) is an eigenvalue problem; therefore, ${\mathbf{c}}^{*}$ is the eigenvector associated with the largest eigenvalue ${\upsilon }^{*}$ of ${\rm Y}$. Of all the undetectable fault directions, the fault direction that causes the largest estimation error is

$$\begin{array}{c}\hfill {\overrightarrow{\mathbf{f}}}_w={\mathbf{D}}_{\pi }{\mathbf{W}}_q{\mathbf{c}}^{*};\end{array}$$

(57)

and the worst-case fault-induced estimation error (per unit fault) is: ${max}_{\mathbf{c}}(\parallel \delta {\mathit{x}}_{\mathbf{f}}{\parallel }^2)={\upsilon }^{*}.$

8.5. Effect of Number of Faults on Failure Mode Slope

Whether $g_f$ is finite (i.e., ${\Delta }_{\pi }$ is singular) depends on (1) how many measurements have faults, h, and (2) the structure of the measurement matrix $\mathbf{H}$. The discussion below summarizes how the failure mode slope changes with the number of faults h.

(a)

When $1\le h\le (m-n)$, ${\Delta }_{\pi }$ can be either singular or nonsingular. When nonsingular, $\mathrm{rank}\left({\Delta }_{\pi }\right)=h$, the worst-case fault direction is given by Equation (54), and $g_f<\infty$. When singular, $\mathrm{rank}\left({\Delta }_{\pi }\right)<h$, the worst-case fault direction is given by Equation (57) and the failure mode slope is unbounded, i.e., $g_f=\infty$.

(b)

As h increases toward $(m-n)$, the numerator in Equation (40) is bounded above by the squared reciprocal of the smallest singular value (i.e., $\frac{1}{{\alpha }_n^2}$), while the (worst-case) denominator can decrease toward zero, which causes the failure mode slope to increase toward infinity.

(c)

When $h>(m-n)$, ${\Delta }_{\pi }$ is singular. Therefore, there is at least one fault direction, as defined in Equation (57), that will make $\parallel {\mathbf{r}}_{\mathbf{f}}{\parallel }^2=0$. As a result, $g_f=\infty$.

(d)

For $h=m$, ${\mathbf{D}}_{\pi }={\mathbf{I}}_m$ and ${\Delta }_{\pi }=\mathbf{Q}$. This ${\Delta }_{\pi }$ has $(m-n)$ eigenvalues values that are one and n eigenvalues that are zero. The eigenvectors corresponding to the zero eigenvalues are in $\mathrm{span}\left({\mathbf{U}}_1\right)$, which is the null space of ${\mathbf{U}}_2$. As stated in Section 6, any fault in $\mathrm{span}\left({\mathbf{U}}_1\right)$ is not detectable from the residual and has $g_f=\infty$. In particular, the worst-case fault direction is ${\overrightarrow{\mathbf{f}}}_w={\mathbf{u}}_1^n$, which is undetectable and affects the state estimation error the most. This is the same solution as that in Equation (57).

9. Effect of Fault on Horizontal Position

The previous sections analyzed the general case for the full estimation vector $\mathit{x}$. In some applications, the analyst might only be interested in a portion of the estimation vector $\mathbf{z}$ that is linearly related to $\mathit{x}$. For example, in automotive transportation applications, there is often a focus on the horizontal position. In this case, the goal is to analyze how the worst-case fault affects the horizontal position, which satisfies $\mathbf{z}=\mathbf{M}{}_E^T\mathbf{R}\mathit{x}$, where $\mathbf{M}$ and ${}_E^T\mathbf{R}$ are defined in (3).

The failure mode slope associated with the horizontal position is defined as:

$$\begin{array}{c}\hfill g_f\triangleq \frac{\parallel \delta {\mathbf{z}}_{\mathbf{f}}{\parallel }^2}{\parallel {\mathbf{r}}_{\mathbf{f}}{\parallel }^2}.\end{array}$$

(58)

Because both $\mathbf{y}$ and $\widehat{\mathbf{y}}$ are unchanged despite the transformation on the estimation vector, the measurement residual $\mathbf{r}$, defined in (19), remains unchanged. The horizontal position error is:

$$\begin{array}{c}\hfill \delta \mathbf{z}=\mathbf{M}{}_E^T\mathbf{R}\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\left(\mathbf{f}+\eta \right),\end{array}$$

(59)

which makes the fault-induced and noise-dependent portions of the horizontal position error, respectively:

$$\begin{array}{c}\hfill \delta {\mathbf{z}}_{\mathbf{f}}=\mathbf{M}{}_E^T\mathbf{R}\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\mathbf{f}\mathrm{and}\delta {\mathbf{z}}_{\eta }=\mathbf{M}{}_E^T\mathbf{R}\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }\eta .\end{array}$$

The worst-case fault direction for the horizontal position error can be found by the same methods used in the previous subsections. The only difference is that the matrix $\Phi$ used to compute $\Gamma$ in Section 8.2 becomes:

$$\begin{array}{c}\hfill \Phi ={\mathbf{U}}_1{\Sigma }^{-1}{\mathbf{V}}^{\top }{}_E^T{\mathbf{R}}^{\top }{\mathbf{M}}^{\top }\mathbf{M}{}_E^T\mathbf{R}\mathbf{V}{\Sigma }^{-1}{\mathbf{U}}_1^{\top }.\end{array}$$

(60)

The norm of the fault-induced horizontal position error is

$$\begin{array}{c}\hfill \parallel \delta {\mathbf{z}}_{\mathbf{f}}{\parallel }^2={\mathbf{f}}^{\top }\Phi \mathbf{f},\end{array}$$

(61)

and the analysis in Section 8 for $g_f$ defined in (32) could be repeated with $g_f$ defined as in Equation (58).

10. Example Discussion

This section provides examples demonstrating selected topics from the paper in the context of GNSS applications.

10.1. Fixed Number of Measurements

Consider the measurement matrix $\mathbf{H}$ used in [43,49]:

$$\begin{array}{c}\hfill \mathbf{H}=\left[\begin{array}{cccc}0.7460266527& -0.4689257437& 0.4728137904& 1\\ -0.8607445743& -0.3446039300& 0.3746557209& 1\\ 0.2109370676& 0.3502943374& 0.9125784518& 1\\ -0.0619331319& -0.4967359072& 0.8656891623& 1\\ -0.7248969588& 0.4759681238& 0.4979746422& 1\\ -0.4009266538& 0.1274180997& 0.9072058455& 1\end{array}\right].\end{array}$$

(62)

The number of measurement is $m=6$. The number of variables to be estimated is $n=4$. The pseudorange noise vector $\eta$ is considered to be zero-mean Gaussian with standard deviation ${\sigma }_{\eta }=3.30$ m, i.e., $\eta \sim \mathcal{N}(\mathbf{0},{\sigma }_{\eta }^2)$. The position vector portion of $\mathit{x}$ is represented in a local tangent plane.

For single-measurement faults, $h=1$, the fault is written as ${\mathbf{f}}_i=\mu {\mathbf{e}}_i\in {\mathbb{R}}^6$, where ${\mathbf{e}}_i$ is the i-th standard basis vector in the ${\mathbb{R}}^6$ for $i=1,\cdots ,6$. The squared failure mode slope $g_f^i$ relating $\delta {\mathbf{z}}_{\mathbf{f}}$ to ${\mathbf{r}}_{\mathbf{f}}$ for each i is computed by Equation (41), as adapted for the horizontal case defined in Equation (58). The slope values are:

$$\sqrt{g_f^i}=\{2.144,1.099,0.917,1.228,1.199,0.275\},$$

for $i=1,\cdots ,6$. Figure 1 plots the horizontal position error norm versus the norm of the residual vector. This axes format matches that used in, e.g., refs. [10,11,43]. For each satellite, Figure 1 displays two items: (1) lines with slope $\sqrt{g_f^i}$, and (2) samples of $\parallel \delta \mathbf{z}\parallel$ computed for three values of $\mu$ and $N=1000$ instances of noise $\eta$. For each instance, the $\eta ,\mu$, and ${\mathbf{e}}_i$ are used to construct a measurement vector $\mathbf{y}$ using (1), which is solved for $\parallel \delta \mathbf{z}\parallel$ using (59) and the residual-based test statistic $\parallel \mathbf{r}\parallel$ is computed using (24). Because $\eta$ is a Gaussian random vector, simulating N measurements with the same ${\mathbf{f}}_i$ produces a point cloud of $(\parallel \mathbf{r}\parallel ,\parallel \delta \mathbf{z}\parallel )$ values, each depicted by a dot. Each point cloud is an ellipsoid of which the center is the mean $(\parallel \delta {\mathbf{r}}_{\mathbf{f}}\parallel ,\parallel \delta {\mathbf{z}}_{\mathbf{f}}\parallel )$ (for a given $\mu$). For the smallest value of $\mu$, the six clusters are very near the origin. As $\mu$ increases, the cluster corresponding to fault direction ${\mathbf{e}}_i$ moves away from the origin along the line with slope $\sqrt{g_f^i}$, as the theory predicts.

The largest failure mode slope occurs for $i=1$ (i.e., green). This means that a fault on the measurements from satellite 1 will cause the largest change in the estimated vector per unit increase in $\parallel \delta {\mathbf{r}}_{\mathbf{f}}\parallel$.

This does not, however, mean that this fault will produce the largest estimated position error for a given value of $\mu$, because its effect on both $\parallel \delta {\mathbf{z}}_{\mathbf{f}}\parallel$ and $\parallel \delta {\mathbf{r}}_{\mathbf{f}}\parallel$ can be small, while $g_f$ is their ratio.

Table 1. Values for means of $\parallel \delta {\mathbf{z}}_{\mathbf{f}}\parallel$ and $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$, and $g_f$ for single-measurement faults, $h=1$, and $\mu =1$ using $\mathbf{H}$ in (62).

i	$\parallel \mathit{\delta }{\mathit{z}}_{\mathit{f}}{\parallel }^2$	$\parallel {\mathit{r}}_{\mathit{f}}{\parallel }^2$	${\mathit{g}}_{\mathit{f}}$
1	$0.3496$	$0.0761$	$4.5955$
2	$0.3330$	$0.2755$	$1.2087$
3	$0.3479$	$0.4139$	$0.8405$
4	$0.5270$	$0.3496$	$1.5078$
5	$0.4367$	$0.3036$	$1.4382$
6	$0.0441$	$0.5813$	$0.0758$

shows the values of $\parallel \delta {\mathbf{z}}_{\mathbf{f}}\parallel$, $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$, and $g_f$ for $\mu =100$. A fault on a measurement from satellite 4 produces the largest position error, but it also produces a proportionally larger test statistic. The value of $\mu$ affects both the numerator and denominator of $g_f$ but not their ratio and, hence, not the value of $g_f$.

Table 2. Values of $\parallel \delta {\mathbf{z}}_{\mathbf{f}}\parallel$, $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$, and $g_f$ for worst-case fault modes when $h=1,\cdots ,6$ and $\mu =1$ for $\mathbf{H}$ defined in Equation (62). The matrix $\Delta$ is singular for fault scenarios marked with an asterisk (${}^{*}$).

h	$\parallel \mathit{\delta }{\mathit{z}}_{\mathit{f}}{\parallel }^2$	$\parallel {\mathit{r}}_{\mathit{f}}{\parallel }^2$	${\mathit{g}}_{\mathit{f}}$	Faulty Measurements
$0$	0	0	$\mathrm{Undefined}$	$\left[\right]$
$1$	$0.3496$	$0.0761$	$4.5955$	$\left[1\right]$
$2$	$0.3925$	$0.0085$	$46.2977$	$[1,6]$
3 *	$1.1456$	0	∞	$[3,4,5]$
4 *	$1.4856$	0	∞	$[2,3,4,5]$
5 *	$1.5028$	0	∞	$[1,2,3,4,5]$
6 *	$1.5254$	0	∞	$[1,2,3,4,5,6]$

shows the worst-case squared failure mode slopes for different h-fault scenarios. In addition, it shows the fault-induced horizontal position error and residual squared norms, computed for $\mu =1$. Due to the fact that the largest single measurement fault slopes in Table 1 correspond to measurements $i=1$ and $i=4$, one might assume that when $h=2$, the worst-case failure mode slope would occur when measurements from these two satellites are affected. This turns out not to be the case. The largest $g_f$ corresponds to faults on satellites 1 and 6 with coefficients $\mathbf{s}=\left[\begin{array}{cc}0.9352& -0.3541\end{array}\right]$ and the overall fault direction, as defined in (43).

10.2. Multi-Measurement Faults: Increasing m

Consider a matrix $\mathbf{H}$ consisting of $m=20$ GNSS measurements. This matrix is created using randomly selected elevation and azimuth angles, denoted by $E_i$ and $A_i$, respectively. To reflect practical GNSS situations, $E_i\in [{15}^{\circ },{90}^{\circ }]$ and $A_i\in [0^{\circ },{360}^{\circ }]$ for $i=1,\dots ,m$. Each row of $\mathbf{H}$ represents a new measurement with ${\mathbf{h}}_i=\left[\begin{array}{cc}{\mathbb{1}}_i& 1\end{array}\right]$, where the unit line-of-sight vector is

$${\mathbb{1}}_i=\left[\begin{array}{ccc}cos{E}_{i}sin{A}_i& cos{E}_{i}cos{A}_i& sin{E}_i\end{array}\right].$$

The matrix $\mathbf{H}$ is expressed as:

$$\begin{array}{c}\hfill \mathbf{H}={\left[{\mathbf{h}}_1^{\top },\cdots ,{\mathbf{h}}_m^{\top }\right]}^{\top }.\end{array}$$

(63)

As m is increased, the first $(m-1)$ rows remain the same, while one new row is added.

Using this $\mathbf{H}$, for the first six fault scenarios (i.e., $h=1,\dots ,6$), each column of

Table 3. Values of the worst-case $g_f$ for the first six fault scenarios. As m increases, the worst-case $g_f$ decreases uniformly for a fixed h. As h increases, the worst-case $g_f$ increases uniformly for a fixed m.

m	$\mathit{h}=1$	$\mathit{h}=2$	$\mathit{h}=3$	$\mathit{h}=4$	$\mathit{h}=5$	$\mathit{h}=6$
6	29.78	1813.29	∞	∞	∞	∞
7	1.85	54.32	22,545.14	∞	∞	∞
8	1.22	3.5	53.41	22,729.26	∞	∞
9	0.94	3.14	24.59	63.98	4.86 × 10${}^5$	∞
10	0.73	2.95	10.34	48.22	803.44	4.69 × 10${}^5$
11	0.64	2.83	7.59	21.47	51.57	9597.89
12	0.43	1.06	2.94	10.11	21.46	594.47
13	0.43	1	1.89	7.63	12.29	26.46
14	0.43	1	1.9	7.32	11.17	23.49
15	0.45	0.91	1.69	5.8	9.82	22.46
16	0.33	0.63	1.24	3.71	5.89	11.52
17	0.33	0.59	1.21	3.21	5.31	11.12
18	0.11	0.42	0.68	1.31	3.32	5.4
19	0.11	0.22	0.52	0.78	1.42	3.43
20	0.1	0.21	0.49	0.75	1.27	3.02

shows the worst-case $g_f$ for a fixed value of h as a function of m. Each column of Table 3 is graphed in Figure 2 as a dotted line of a unique color. The color code corresponding to different values of h is defined along the right side of the figure.

Note that, for each fixed h (i.e., any column), the value of $g_f$ decreases as m increases. Section 6 discussed that the detection threshold $T^{*}$ can be set as a function of the worst-case failure mode slope $g_f^{*}$. The fact that $g_f$ decreases for fixed h as m increases allows the designer to reduce integrity risk for the h-fault scenario by increasing the number of measurements used.

10.3. Multi-Measurement Faults: Increasing h

Each row of Table 3 shows the worst-case $g_f$ for a fixed value of m as a function of h. Note that for each fixed m (i.e., any row), the value of $g_f$ increases with h. The fact that $g_f$ increase for fixed m as h increases shows that integrity risk increases with higher h-fault scenarios.

For a given $\mathbf{H}$, m is constant. The reason why the worst-case $g_f$ increases as a function of h is explained as follows. As h increases to $h^{\prime }=h+1$, the matrix ${\mathbf{D}}_{\pi }$ gains an additional column and the vector $\mathbf{s}$ gains an additional row. This allows more flexibility in choosing the worst-case fault direction, which maximizes $g_f$. Whatever direction caused the worst-case fault for h is still a viable fault for $h^{\prime }$ (just set the coefficient for the new row of $\mathbf{s}$ to zero); however, the additional column of ${\mathbf{D}}_{\pi }$ may allow the optimization defined in Section 8.3 to find a new fault direction such that $g_f\left(h^{\prime }\right)\ge g_f\left(h\right)$. The fact that the worst-case $g_f$ increases with h is shown by Table 2 for the specific $\mathbf{H}$ in Equation (62); each row in Table 3 for the $\mathbf{H}$ defined in Equation (63); and the colored dots at any fixed value of m in Figure 2.

Figure 3 separately graphs the numerator $\parallel {\mathit{x}}_{{\overrightarrow{f}}_w}\parallel$ (blue) and denominator $\parallel {\mathbf{r}}_{{\overrightarrow{f}}_w}\parallel$ (red) of Equation (40) as a function of h for three values of m (dashed for $m=20$, dotted for $m=21$, and solid for $m=22$). For each m, as h increases, the value of $\parallel {\mathit{x}}_{{\overrightarrow{f}}_w}\parallel$ tends toward, but is bounded above by $\frac{1}{{\alpha }_n}$, where, for this example, ${\alpha }_n=1.938$ for $m=22$. The increase toward the upper bound is not monotonic. The worst-case fault ${\overrightarrow{f}}_w$ is defined based on the ratio $g_f$ defined in Equation (32). This ratio can increase even if the numerator decreases with h, as long as the denominator decreases at least a proportionate amount. For each m, as h increases, the value of $\parallel {\mathbf{r}}_{{\overrightarrow{f}}_w}\parallel$ tends (not monotonically) toward zero. As the $\parallel {\mathbf{r}}_{{\overrightarrow{f}}_w}\parallel$ approaches zero, $g_f$ approaches infinity. The fact that $g_f$ increases toward infinity is not due to $\parallel \delta {\mathit{x}}_{\mathbf{f}}\parallel$ becoming very large but due to $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$ approaching zero.

10.4. How Does $g_f$ Become Infinite for $h\le (m-n)$?

Define $h^{*}$ as the value of h at which ${\Delta }_{\pi }$ becomes nonsingular for a given $\mathbf{H}$. For $h\ge h^{*}$, there exist faults that do not affect the residual and, therefore, are undetectable from it for any threshold $T^{🟉}$.

Consider the $\mathbf{H}$ in Equation (63) that produced the results in Table 3. For $m=6$ and $(m-n)=2$, $g_f$ is finite for $h=1$ and 2 and infinite for $h>2$; therefore, $h^{*}=3$. For $m=7$ and $(m-n)=3$, $h^{*}=4$ and there are 3 finite values of $g_f$. For all values of m, there are, at most, $(m-n)$ finite values of $g_f$ since ${\Delta }_{\pi }$ is singular for all $h>(m-n)$. From Table 3, for m = 6 to 10, $h^{*}=(m-n)+1$, which is the known case [10].

As shown in Section 8.3, it is possible to have $h^{*}\le (m-n)$. The gray line in Figure 2a corresponds to $h=(m-n)$. When $m\ge 11$, there are only $(m-n)-1$ finite values of $g_f$ and, hence, $h^{*}=(m-n)$. The brown curve in Figure 2a corresponds to $g_f$ values for $h=(m-n)-1$. The general (though not absolute) trend of $g_f$ values on these curves is increasing with m, suggesting that as m increases the dimension of the subspace of undetectable faults may further increase, thus making $h^{*}=(m-n)-1$ as m increases.

This shows that increasing the number of measurements m results in a trade-off. For a given h, $g_f\left(m\right)$ decreases as m increases; however, $g_f(m-n)$ is non-decreasing with m, and $h^{🟉}\left(m\right)$ may become less than or equal to $(m-n)$.

11. Conclusions

This paper analyzed the effect of measurement faults on the estimation error, measurement residual, and the failure mode slope in linear over-determined systems, such as GNSS. Of particular interest is the case where m is large because the number of GNSS measurements is increasing. The singular value decomposition provides formulas for the fault-induced estimation error and residual vector (see Equations (10) and (20), respectively) in terms of the orthogonal basis (i.e., ${\mathbf{U}}_1$) for the range of $\mathbf{H}$ and the complement of the range space (i.e., ${\mathbf{U}}_2$). This decomposition provides convenient formulas for the failure mode slope and worst-case fault in Equations (36) and (40). Using this decomposition, Section 7 provides ${\mathbf{U}}_1$ as a basis for the n dimensional linear space of faults that are undetectable from the residual (or parity) vector and ${\mathbf{U}}_2$ as the basis for the linear space of faults that affect the residual vector without affecting the estimate. For each fault scenario affecting h measurements, finding the worst-case fault is known to be an eigenvalue problem. Herein, this eigenvalue problem is written in terms of these two linear spaces, as in Equations (51) and (53). The detectability of any fault in the h-fault scenario from the residual vector depends on whether the range space of the matrix ${\mathbf{D}}_{\pi }$, which is a permutation of h columns of the identity matrix, has a component in the linear space spanned by ${\mathbf{U}}_1$. Building on this idea, it is clear from Equation (40) that, as h increases, the numerator is bounded above by the squared reciprocal of the smallest singular value (i.e., $\frac{1}{{\alpha }_n^2}$), while the (worst-case) denominator can decrease toward zero, which causes the failure mode slope to increase toward infinity. Alternatively, for a fixed h, as the number of measurements m increases, the dimension of the vector $\parallel {\mathbf{r}}_{\mathbf{f}}\parallel$ in the denominator increases, which tends to decrease the failure mode slope for h. Simulated GNSS examples are included to demonstrate the results derived herein.