Kalman Filter Innovation-Based Optimal Integrity Monitoring Against Unknown Bias in Time Updates

Abstract

Kalman-filter-based multi-sensor navigation systems are promising for future mobility but vulnerable to sensor faults, necessitating integrity monitoring (IM). This paper proposes an optimal innovation-based IM approach that minimizes the integrity risk against dynamic sensor faults using current-time information and a single filter. While the solution separation (SS) concept provides an optimal solution yielding the snapshot-based SS IM, we approach this problem through innovation-domain optimization. We formulate a min-max problem, reparametrize it to apply the Neyman–Pearson lemma, and obtain the position-domain projection of the innovation as the optimal statistic. We prove that this solution is mathematically equivalent to the snapshot-based SS statistic. Performance evaluation of the innovation-based position-domain IM method under IMU faults demonstrates superior integrity performance compared with innovation-based range-domain IM. Comparison with the full SS IM approach using parallel filters shows the advantages of the proposed approach in improving computational efficiency while narrowing the integrity performance gap as measurement redundancy increases.

Keywords

1 INTRODUCTION

Multi-sensor navigation systems have become essential for future mobility systems, including autonomous cars and advanced air mobility. These systems integrate various navigation sensors with complementary characteristics through recursive filters, such as the Kalman filter (KF), to compensate for their individual limitations. The KF provides an unbiased and optimal estimate in terms of mean-squared error under nominal conditions by utilizing a knowledge of sensor measurement characteristics, which has led to the development of numerous KF-based multi-sensor navigation systems, such as those reported by Xu et al. (2022), Liao et al. (2021), Afia et al. (2015), Travis et al. (2005), and Miller et al. (1993). However, these systems are vulnerable to sensor failures, which introduce unknown biases into measurements. These faults can lead to significant position errors and pose a safety risk. Therefore, navigation integrity, which refers to a system's ability to provide timely warnings when position estimates should not be used for navigation purposes, is required in KF-based navigation systems when utilized for safety-critical applications.

The integrity monitoring (IM) procedure protects the user against faults by performing two primary functions: fault detection and integrity risk evaluation. In the fault detection step, a test statistic is formulated from the measurements. If the test statistic falls outside the detection region, as defined by Blanch et al. (2017), an alert is triggered to notify users not to use the navigation solution. The size of the detection region is set to maintain the false alert rate (or continuity risk) below the required level. Several prior works have proposed fault detection algorithms that, although not specifically developed for navigation integrity, form the foundation of the fault diagnosis research field. Patton and Chen (1991) and Odendaal and Jones (2014) utilized a parity space (or residual) vector for fault diagnosis. Hajiyev and Caliskan (1998) employed a sequence of KF innovation vectors to detect faults in dynamic systems. Although different fault detection algorithms have been proposed, they are exposed to the possibility of failing to detect faults due to nominal measurement errors, potentially leading to large state errors.

In this context, the second primary function for IM, referred to as the integrity risk evaluation, quantifies the impact of undetected faults on position estimates. Integrity risk is defined as the probability that an undetected position error exceeds a predefined acceptable limit, also known as the alert limit (AL). This risk is assessed against a predefined integrity risk requirement to determine whether to abort or continue a mission. The integrity risk evaluation involves both assessing the monitor's fault detection capability and quantifying the impact of undetected faults on position estimation errors. Several research efforts have proposed methods for evaluating the integrity risk of KF-based navigation systems (Joerger & Pervan, 2013; Gunning et al., 2018; Wang et al., 2022; Lee et al., 2023).

For unknown biases introduced in the measurement-update step of a KF, Joerger and Pervan (2013) proposed a residual-based (RB) IM algorithm that uses a weighted norm of the residual of a batch weighted least-square (WLS) estimator as the test statistic. The rationale behind this approach is that the current-time state estimates are identical for both a KF and a batch WLS estimator. This approach targets global navigation satellite system (GNSS) faults in KF-based GNSS/inertial measurement unit (IMU) navigation systems. The authors derived the worst-case fault profile, comprising all GNSS fault vectors that have occurred throughout the entire period of KF operation, that maximizes the integrity risk when using this RB test statistic and evaluated the corresponding integrity risk. Because the worst-case fault profile must be recomputed at each epoch to account for the recursive nature of the KF, this method is computationally intensive.

In contrast, Gunning et al. (2018) introduced a solution separation (SS) IM algorithm for a KF-based precise point positioning system by extending the SS receiver autonomous IM (RAIM) concept originally developed for GNSS-based navigation systems. For each fault hypothesis, the SS algorithm constructs a test statistic as the difference between the position estimate obtained with all available measurements and the estimate obtained after removing the measurement corresponding to the fault hypothesis (i.e., the fault-free solution). Accordingly, this approach requires maintaining a parallel KF for each fault hypothesis. Gunning et al. (2018) also compared the performance of their SS IM method with that of a batch RB IM method and demonstrated that the SS approach outperforms the RB method. The authors also showed that the RB statistic can be formulated using innovation vectors. Wang et al. (2022) subsequently proposed an SS IM for a centralized KF-based GNSS/multiple-IMU navigation system, addressing both GNSS and IMU fault hypotheses.

Although these IM algorithms are well established for ensuring navigation integrity, they suffer from high computational complexity. RB (or innovation-based) IM requires recomputing the worst-case fault profile at every epoch, involving matrix inversions whose cost grows over time. Moreover, the detection capability degrades as measurement noise accumulates. SS IM demands a separate filter for every fault hypothesis, which becomes burdensome as the number of hypotheses increases. Such complexity makes these methods impractical for applications with strict computational constraints. These limitations motivated us to develop an IM algorithm for dynamic sensor faults that relies solely on a single filter and current-time information, thereby reducing the computational load while ensuring navigation integrity.

Our prior work by Lee et al. (2023) proposed IM with a single filter, utilizing the magnitudes of the elements of the KF innovation vector as test statistics, to detect the faults of IMU sensors used in time-update steps. This monitor issues an alert if the magnitude of any innovation element falls outside the detection region. Although this monitor is designed for IMU failures, it is applicable to other sensors used in the time update because it is based on the fundamental characteristics of the KF, which do not depend on sensor type. This approach is simple to implement because it utilizes the current innovation vector, which is inherently available in the KF. Furthermore, Lee et al. (2023) showed that the current innovation vector fully captures the impact of the undetected fault profile on the position estimate. While it is straightforward to implement the approach of Lee et al. (2023) and evaluate the integrity risk, this approach is not optimal in terms of the availability, defined as the fraction of time during which the system is available (i.e., the integrity and continuity risks are below their corresponding requirements). In response, this paper aims to design a test statistic, as a function of the current-time innovation vector, that minimizes the integrity risk under the dynamic sensor failure condition.

The SS concept, which compares a fault-free solution computed without faulty measurements to a solution computed using all available measurements, provides the optimal solution for designing the optimal innovation-based statistic under dynamic sensor faults. The current-time innovation vector can be decomposed into two components: one component is a function of an a priori estimate directly affected by dynamic sensor faults, and the other component is a function of current fault-free measurements. Therefore, a fault-free solution can be constructed from the innovation vector using only the measurement component (i.e., a snapshot solution). Following the SS concept, the optimal innovation-based statistic is defined as the difference between this snapshot solution and an a posteriori estimate, which we refer to as the snapshot-based SS statistic in this paper.

Although the equivalence between the innovation-based optimal statistic and the snapshot-based SS statistic may be evident based on prior work (Blanch et al., 2017) demonstrating that the SS concept provides the optimal statistic against a single multi-dimensional fault, establishing this equivalence requires rigorously solving the underlying optimization problem. To this end, we formulate a min-max optimization problem for deriving the optimal innovation-based statistic. This problem involves minimization over the detection region defined solely in the innovation domain and maximization over dynamic sensor faults to account for the worst-case fault scenario. We then reparametrize this problem into a form compatible with the procedure originally developed for RAIM by Blanch et al. (2017), which was used to obtain the exact solution for a similar min-max problem. This reparameterization enables us to derive the optimal innovation-based statistic, and we verify that the procedure of Blanch et al. (2017) guarantees the optimal solution for our problem. Notably, the direct solution of this min-max problem takes the form of a position-domain projection of the innovation, rather than an explicit SS formulation.

To establish the rigorous equivalence between these two seemingly different formulations, we detail the construction of the snapshot-based SS statistic within the KF architecture and provide a formal mathematical proof demonstrating that it is exactly equivalent to the position-domain projection of the innovation. This formal proof confirms that the same result can be reached through innovation-domain optimization.

The remainder of this paper is organized as follows. Section 2 provides background on the KF and the prior method proposed by Lee et al. (2023). Section 3 formulates the optimization problem for deriving the innovation-based optimal test statistic that minimizes the integrity risk and derives the corresponding optimal solution, which is a position-domain projection of the innovation vector. Additionally, an SS statistic is introduced using a current-time snapshot solution as the fault-free solution within the KF architecture, and the theoretical equivalence between the position-domain projection of the innovation and this snapshot-based SS statistic is rigorously demonstrated. Section 4 evaluates the integrity performance of the proposed monitor in terms of integrity risk under the IMU fault condition for a KF-based Global Positioning System (GPS)/IMU navigation system. Section 5 compares the integrity performance and computational load of the snapshot-based SS and full SS IM methods to identify when the proposed approach offers practical advantages. Section 6 summarizes the paper.

2 BACKGROUND

This section provides background for this paper. Sections 2.1 and 2.2 detail the KF estimates under fault-free conditions and under sensor fault conditions with an unknown bias in the time update, respectively. Section 2.3 reviews the IM algorithm proposed by Lee et al. (2023), which utilizes the KF innovation-based range-domain monitor.

2.1 KF Estimate Error Under Fault-Free Conditions

This subsection describes errors of extended KF (EKF) position estimates under the nominal condition. The EKF sequentially estimates the state in two steps: a time-update step utilizing dynamic measurements and a measurement-update step utilizing observation measurements.

In the time-update step, the discretized dynamic model at epoch k is given as follows:

$x_{k} = F_{k - 1} x_{k - 1} + G_{k - 1} u_{k - 1} + w_{k - 1}$ 1

where $x_{k} \in R^{n}$ is the state vector, $F_{k - 1} \in R^{n \times n}$ is the state transition matrix, $G_{k - 1} \in R^{n \times l}$ is the input coefficient matrix, $u_{k - 1} \in R^{l}$ is the known control input, and $w_{k - 1} \in R^{n}$ is the process noise vector. The vector $w_{k - 1}$ is assumed to be normally distributed with zero mean and covariance matrix $Q_{k - 1} \in R^{n \times n}$ under the nominal condition. The following notation is used:

$w_{k - 1} \sim N (0_{n \times 1}, Q_{k - 1})$ 2

where the notation of $N (μ, Σ)$ represents the normal Gaussian distribution with mean μ and covariance matrix Σ. The vector $w_{k - 1}$ represents white noise, implying that $E [w_{k} w_{k}^{T},] = 0$ for $k \neq k^{'}$ , where $E [■]$ denotes the expectation operator.

Based on Equation (1), the time-updated state estimate and its corresponding covariance matrix are presented in Equations (3) and (4), respectively:

${\hat{x}}_{k}^{-} = F_{k - 1} {\hat{x}}_{k - 1}^{+} + G_{k - 1} u_{k - 1}$ 3

$P_{k}^{-} = F_{k - 1} P_{k - 1}^{+} F_{k - 1}^{T} + Q_{k - 1}$ 4

where ${\hat{x}}_{k}^{-} \in R^{n}$ is an a priori estimate (i.e., time-updated estimate) of state $x_{k}$ at epoch $k, {\hat{x}}_{k - 1}^{+} \in R^{n}$ is an a posteriori estimate (i.e., measurement-updated estimate) at epoch $k - 1, P_{k}^{-} \in R^{n \times n}$ is the covariance matrix of ${\hat{x}}_{k}^{-}$ , and $P_{k - 1}^{+} \in R^{n \times n}$ is the covariance matrix of ${\hat{x}}_{k - 1}^{+}$ . The control input vector $u_{k - 1}$ is constructed from dynamic sensor measurements, such as accelerometer, gyroscope, and odometer measurements. The estimate error of ${\hat{x}}_{k}^{-}$ , denoted by ${\tilde{x}}_{k}^{-}$ , follows a zero-mean Gaussian distribution under the nominal condition:

${\tilde{x}}_{k}^{-} = x_{k} - {\hat{x}}_{k}^{-} \sim N (0_{n \times 1}, P_{k}^{-})$ 5

In the measurement-update step, the linearized measurement model at epoch k is given as follows:

$y_{k} = H_{k} x_{k} + v_{k}$ 6

where $y_{k} \in R^{m}$ is the normalized observation measurement vector, $H_{k} \in R^{m \times n}$ is the normalized and linearized observation matrix, and $v_{k} \in R^{m}$ is the normalized measurement noise vector composed of zero-mean, unit-variance independent and identically distributed Gaussian random variables. The vector $y_{k}$ consists of measurements obtained from sensors such as GNSS, lidar, and altimeter sensors.

The measurement-updated state and its covariance matrix are described in Equations (7) and (8), respectively:

${\hat{x}}_{k}^{+} = {\hat{x}}_{k}^{-} + K_{k} (y_{k} - H_{k} {\hat{x}}_{k}^{-})$ 7

$P_{k}^{+} = (I_{n} - K_{k} H_{k}) P_{k}^{-}$ 8

where $I_{n} \in R^{n \times n}$ is the identity matrix and $K_{k} \in R^{n \times m}$ is the Kalman gain matrix, defined as follows:

$K_{k} = P_{k}^{-} H_{k}^{T} {(H_{k} P_{k}^{-} H_{k}^{T} + I_{m})}^{- 1}$ 9

The estimate error of ${\hat{x}}_{k}^{+}$ , denoted by ${\tilde{x}}_{k}^{+}$ , also follows a zero-mean Gaussian distribution under the nominal condition, as follows:

${\tilde{x}}_{k}^{+} = x_{k} - {\hat{x}}_{k}^{+} \sim N (0_{n \times 1}, P_{k}^{+})$ 10

2.2 KF Estimate Error Under Sensor Fault Conditions

This subsection demonstrates the impact of unmodeled biases in the time-update step on the KF estimates. In the presence of a dynamic sensor fault, an unmodeled bias ( $Δ u_{k - 1}$ ) is introduced into the control input, modifying the time-update equation in Equation (3) as follows:

${\hat{x}}_{k}^{-} = F_{k - 1} {\hat{x}}_{k - 1}^{+} + G_{k - 1} (u_{k - 1} + Δ u_{k - 1})$ 11

In the remainder of this paper, $G_{k - 1} Δ u_{k - 1}$ is denoted as $f_{k - 1} \in R^{n}$ . This vector $f_{k - 1}$ introduces biases into both a priori and a posteriori estimates. Additionally, these estimates are affected by faults from previous epochs as well as the current epoch. Therefore, a priori and a posteriori estimate errors, under dynamic sensor fault conditions, are modeled as follows:

${\tilde{x}}_{k}^{-} \sim N (F_{k - 1} d_{k - 1}^{+} + f_{k - 1}, P_{k}^{-})$ 12

${\tilde{x}}_{k}^{+} \sim N (d_{k}^{+}, P_{k}^{+})$ 13

where $d_{k}^{+} \in R^{n}$ is the vector of $L_{k} (F_{k - 1} d_{k - 1}^{+} + f_{k - 1})$ or $L_{k} d_{k}^{-}, L_{k}$ is the matrix of $I_{n} - K_{k} H_{k}, d_{k}^{-}$ is the bias in ${\tilde{x}}_{k}^{-}$ , and $d_{k - 1}^{+}$ is the bias in ${\tilde{x}}_{k - 1}^{+}$ . These equations imply that these biases in ${\tilde{x}}_{k}^{-}$ and ${\tilde{x}}_{k}^{+}$ depend on the fault history vector, which comprises all fault vectors that have occurred throughout the entire period of KF operation.

2.3 KF Innovation-Based IM Algorithm

This subsection briefly reviews the work by Lee et al. (2023), which proposed an IM algorithm that utilizes a KF innovation-based range-domain monitor to ensure integrity against IMU faults for KF-based GNSS/IMU navigation systems.

The innovation vector at epoch $k (γ_{k} \in R^{m})$ is as follows:

$γ_{k} = y_{k} - H_{k} {\hat{x}}_{k}^{-}$ 14

Under the nominal condition, the innovation vector follows a zero-mean Gaussian distribution with covariance matrix $S_{k} \in R^{m \times m}$ :

$γ_{k} ∣ H_{0} \sim N (0_{m \times 1}, S_{k})$ 15

where $H_{0}$ represents the nominal condition in which all sensors are fault-free. $S_{k}$ is computed as follows:

$S_{k} = H_{k} P_{k}^{-} H_{k}^{T} + I_{m}$ 16

In the presence of an IMU fault, the innovation vector follows a Gaussian distribution with a mean vector of $b_{k} \in R^{m}$ but the same covariance matrix as in the nominal condition:

$γ_{k} ∣ H_{i (i > 0)} \sim N (b_{k}, S_{k})$ 17

Here, $H_{i (i > 0)}$ represents the i-th hypothesis of sensor faults. For instance, $H_{1}$ represents the hypothesis in which the IMU sensor is faulted, and $H_{2}$ is the hypothesis for GNSS sensor faults. In the remainder of this paper, the subscript i represents only the hypothesis in which the dynamic sensor is faulted, unless otherwise stated.

The impact of the fault on $γ_{k}$ is recursively propagated through the KF. Thus, under the IMU fault hypothesis, $b_{k}$ is expressed as a function of the mean of ${\tilde{x}}_{k - 1}^{+}$ (i.e., $d_{k - 1}^{+}$ ), caused by any fault occurring before the current epoch, and the bias of $f_{k - 1}$ , caused by the fault occurring at the current epoch, as derived by Lee et al. (2023):

$b_{k} = - H_{k} F_{k - 1} d_{k - 1}^{+} - H_{k} f_{k - 1}$ 18

Lee et al. (2023) employed $γ_{k}$ as test statistics to detect IMU sensor faults. The monitor issues an alert if the magnitude of any element in $γ_{k}$ exceeds the predefined monitoring thresholds, as follows:

$⋃_{j} | γ_{k, j} | > T_{k, j}$ 19

where $γ_{k, j}$ is the j-th element of $γ_{k}$ and $T_{k, j}$ is its corresponding monitor threshold.

This innovation-based range-domain monitor proposed by Lee et al. (2023) is not limited to IMU fault detection; rather, it can be applied to detect faults in any sensor used in the KF time-update step. This versatility stems from its reliance on fundamental KF characteristics, rather than specific sensor characteristics. Additionally, this method requires minimal additional computational resources for IM, as the innovation is inherently calculated within the KF. However, Lee et al. (2023) did not aim to determine a monitor test statistic that minimizes the integrity risk and thereby maximizes system availability. The following section develops an optimal KF innovation-based monitor that can minimize the integrity risk against dynamic sensor faults.

3 INNOVATION-BASED OPTIMAL DETECTOR AGAINST DYNAMIC SENSOR FAULTS

This section derives the optimal test statistic for innovation-based monitoring that minimizes the integrity risk against dynamic sensor faults. Although the optimal innovation-based statistic can be obtained from the SS concept by formulating a snapshot-based SS statistic, this paper approaches the problem differently from the perspective of optimizing the detection region in the innovation domain. Section 3.1 formulates a min-max optimization problem that seeks the detection region that minimizes the integrity risk under worst-case dynamic sensor faults, constrained to use only the innovation vector. Section 3.2 reparametrizes this min-max problem into a form compatible with the procedure originally developed for RAIM (Blanch et al., 2017). Section 3.3 derives the optimal innovation-based statistic, which takes the form of a position-domain projection of the innovation, rather than an explicit SS formulation. Section 3.4 then demonstrates that this rigorously derived solution is mathematically equivalent to a snapshot-based SS statistic, thereby confirming that the same result can be reached by optimizing the detection region in the innovation domain.

3.1 Optimization Problem for Detection Regions

This subsection defines the optimization problem for deriving the innovation-based optimal test statistic that minimizes the integrity risk. Designing the innovation-based test statistic can be equivalently cast as determining a detection region ( $Ω_{k, q} \subseteq R^{m}$ ). When the innovation vector ( $γ_{k}$ ) falls outside the detection region (i.e., $γ_{k} \notin Ω_{k, q}$ ), the monitor triggers an alert.

System availability is evaluated by comparing the integrity risk under the i-th hypothesis of dynamic sensor faults ( $H_{i}$ ), defined as the probability that the position error exceeds the AL without triggering an alert (Blanch et al., 2017), against its requirement. Specifically, the system is considered available when the following holds:

$P (| x_{k, q} - {\hat{x}}_{k, q}^{+} | \geq A L_{q}, γ_{k} \in Ω_{k, q} ∣ H_{i}) P_{H_{i}} \leq P_{H M I, H_{i, q}}$ 20

Here, $x_{k, q}$ is the q-th coordinate state of $x_{k}, {\hat{x}}_{k, q}^{+}$ is an a posteriori estimate of $x_{k, q}, A L_{q}$ is the alert limit, $P_{H_{i}}$ is the prior probability of $H_{i}$ , and $P_{H M I, H_{i, q}}$ is the pre-allocated integrity risk requirement under $H_{i}$ for the q-th coordinate state. While the system attempts to meet the required integrity, $Ω_{k, q}$ must be set to maintain the false alert rate, i.e., the continuity risk under the nominal condition ( $H_{0}$ ), below its requirement ( $P_{f a, i, q}$ ):

$P (γ_{k} \notin Ω_{k, q} ∣ H_{0}) \leq P_{f a, i, q}$ 21

Therefore, in the context of maximizing the system availability, determining $Ω_{k, q}$ becomes an optimization problem of minimizing the integrity risk subject to constraints on the false alert risk. This optimization problem can be formulated as follows:

$\begin{matrix} min_{Ω_{k, q}} P (| x_{k, q} - {\hat{x}}_{k, q}^{+} | \geq A L_{q}, γ_{k} \in Ω_{k, q} ∣ H_{i}) : \\ P (γ_{k} \notin Ω_{k, q} ∣ H_{0}) \leq P_{f a, i, q} \end{matrix}$ 22

In Equation (22), evaluating the integrity risk requires the fault profile vector, which comprises all fault vectors that have occurred throughout the entire period of KF operation. This stems from the fact that the probability distributions of both ${\hat{x}}_{k, q}^{+}$ and $γ_{k}$ depend on the fault profile vector. However, prior information on faults is limited owing to their rare occurrence, unlike nominal measurement errors whose statistical characteristics can be reliably modeled using a large amount of data. Therefore, in safety-critical applications, the integrity risk should be evaluated in the most conservative manner to avoid making an optimistic assumption about faults. This conservative approach leads us to use the following maximum value of the integrity risk in Equation (22):

$max_{f_{0 : k - 1}} P (| x_{k, q} - {\hat{x}}_{k, q}^{+} | \geq A L_{q}, γ_{k} \in Ω_{k, q} ∣ H_{i}, f_{0 : k - 1})$ 23

where $f_{0 : k - 1}$ is the fault profile vector under the $H_{i}$ fault hypothesis. The $f_{0 : k - 1}$ vector that maximizes the integrity risk is referred to as the worst-case fault profile vector.

By substituting Equation (23) into Equation (22), the optimization problem can be formulated as the following min-max problem:

$\begin{matrix} min_{Ω_{k, q}} max_{f_{0 : k - 1}} P (| x_{k, q} - {\hat{x}}_{k, q}^{+} | \geq A L_{q}, γ_{k} \in Ω_{k, q} ∣ H_{i}, f_{0 : k - 1}) : \\ P (γ_{k} \notin Ω_{k, q} ∣ H_{0}) \leq P_{f a, i, q} \end{matrix}$ 24

3.2 Reparameterization

This subsection reparametrizes Equation (24) into a form compatible with the procedure established by Blanch et al. (2017) to obtain the optimal statistic for RAIM. Blanch et al. (2017) reparametrized the test statistic defined in the residual (or parity) space into two independent random variables: one variable fully captures GNSS measurement faults (i.e., the fault estimate derived from the residual vector), and the other variable accounts solely for noise. The fault estimate was further reparametrized into a new random variable, whose first component corresponds to the position bias. Using the Neyman–Pearson lemma, Blanch et al. (2017) showed that this first component becomes the optimal statistic and demonstrated that the remaining components do not affect the optimality of the solution to the min-max optimization problem for RAIM.

To apply this procedure for deriving the optimal innovation-based statistic, we reparametrize the innovation vector into two independent random variables: one variable captures the bias in an a priori estimate (i.e., Equation (32)), and the other variable accounts solely for the innovation noise (i.e., Equation (28)). This reparameterization is performed based on the following expression of the innovation vector under the fault hypothesis $H_{i}$ :

$γ_{k} ∣ H_{i} = y_{k} - H_{k} {\hat{x}}_{k}^{-} = v_{k} - H_{k, p} η_{k, p}^{-} - H_{k, p} d_{k, p}^{-}$ 25

where $H_{k, p} \in R^{m \times n_{p}}$ is the partial matrix of $H_{k}$ corresponding to the states whose information is provided by measurements, with the subscript “p” indicating that it is a partial matrix. When only pseudoranges are utilized in the measurement-update step, $H_{k, p}$ is constructed using the column vectors of $H_{k}$ corresponding to the position and receiver clock bias. $d_{k, p}^{-} \in R^{n_{p}}$ is the partial vector of $d_{k}^{-}$ corresponding to the partial states, and $η_{k, p}^{-} \in R^{n_{p}}$ is a random vector defined as follows:

$η_{k}^{-} \sim N (0_{n_{p} \times 1}, P_{k, p}^{-})$ 26

where $n_{p}$ is the number of partial states and $P_{k, p}^{-}$ is the partial matrix of $P_{k}^{-}$ corresponding to the partial state vector.

From Equation (25), the bias of an a priori estimate error induced by dynamic sensor faults $(d_{k}^{-})$ can be partially obtained as follows:

${\hat{d}}_{k, p}^{-} = {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} γ_{k}$ 27

where ${\hat{d}}_{k, p}^{-}$ denotes the bias estimate in ${\tilde{x}}_{k, p}^{-}, d_{k, p}^{-}$ , computed from the innovation vector. Equation (27) is applicable only when a pseudo-inverse of $H_{k, p}$ exists. This paper assumes that a sufficient number of measurements is provided in the measurement-update step, resulting in $H_{k, p}$ being of full rank (i.e., $m > n_{p} - 1$ ).

The part of the innovation that is statistically independent of ${\hat{d}}_{k, p}^{-}$ is computed as follows:

$ε_{⊥, k, p} = (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) γ_{k} = P_{⊥, k, p} γ_{k}$ 28

$ε_{⊥, k, p}$ is also orthogonal to $H_{k, p} {\hat{d}}_{k, p}^{-}$ . The statistical independence between $ε_{⊥, k, p}$ and ${\hat{d}}_{k, p}^{-}$ is proven in Appendix A.

We introduce the random variable ${\hat{δ}}_{k, p, 1}^{+}$ , which is the first component of the random vector ${\hat{δ}}_{k, p}^{+}$ (defined in Equation (32)) as follows:

${\hat{δ}}_{k, p, 1}^{+} = \frac{- e_{q}^{T} L_{k} E_{k} {\hat{d}}_{k, p}^{-}}{\sqrt{Var [- e_{q}^{T} L_{k} E_{k} {\hat{d}}_{k, p}^{-}]}}$ 29

where $e_{q}^{T}$ represents the row vector extracting the q-th coordinate state and $E_{k} \in R^{n \times n_{p}}$ is the matrix constructing the full state vector from the partial state vector. The full state constructed using $E_{k}$ has zero elements for all states except for the position and receiver clock bias when only pseudoranges are used in the measurement-update step. $Var [■]$ denotes the variance operator. After ${\hat{δ}}_{k, p, 1}^{+}$ has been defined, the remaining components of ${\hat{δ}}_{k, p}^{+}$ can be determined, as shown in Equation (32). From this, we can obtain a form (i.e., Equation (35)) for which the result of Blanch et al. (2017), which states that the first component of the defined random variable becomes the optimal statistic, is applicable.

The position bias of the q-th coordinate state ( $d_{k, q}^{+}$ ) can be expressed in terms of the mean of $γ_{k}, b_{k}$ , under the dynamic sensor fault hypothesis as described by Lee et al. (2023):

$d_{k, q}^{+} = - e_{q}^{T} L_{k} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} b_{k}$ 30

From Equations (13), (27), (29), and (30), the distribution of the q-th coordinate state of ${\tilde{x}}_{k}^{+}$ can be expressed as follows:

${\tilde{x}}_{k, q}^{+} \sim N (\sqrt{B_{k, q} S_{k} B_{k, q}^{T}} δ_{k, p, 1}^{+}, σ_{k, q}^{2})$ 31

where $B_{k, q}$ is $- e_{q}^{T} L_{k} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}, S_{k}$ is the covariance matrix of $γ_{k}$ , and $σ_{k, q}$ is the standard deviation (SD) of ${\tilde{x}}_{k, q}^{+}$ .

Finally, ${\hat{δ}}_{k, p}^{+}$ is fully defined as follows:

${\hat{δ}}_{k, p}^{+} = {[\begin{array}{lll} v_{1} & \dots & v_{n_{p}} \end{array}]}^{T} \cdot {\hat{d}}_{k, p}^{-} = V_{k, p} \cdot {\hat{d}}_{k, p}^{-}$ 32

where $v_{1} \in R^{n_{p}}$ is $- E_{k}^{T} L_{k}^{T} e_{q} / \sqrt{B_{k, q} S_{k} B_{k, q}^{T}}$ and $v_{2}, \dots, v_{n_{p}}$ are column vectors computed using the Gram-Schmidt process (Golub & Van Loan, 1996) from $v_{1}$ . The vectors $v_{1}, \dots, v_{n_{p}}$ ensure the following conditions as described by Blanch et al. (2017):

$E [v_{j}^{T} {\hat{d}}_{k, p}^{-} \cdot v_{k}^{T} {\hat{d}}_{k, p}^{-}] = {\begin{matrix} 0, & if j \neq k \\ 1, & if j = k \end{matrix} \forall j, k$ 33

Subsequently, the random variable $η_{k, q}^{+}$ is defined as follows:

$η_{k, q}^{+} \sim N (0, σ_{k, q}^{2})$ 34

From Equations (31), (32), and (34), the reparametrized optimization problem is obtained as follows:

$\begin{matrix} min_{Ω_{k, q}} max_{δ_{k, p}^{+}} P (| η_{k, q}^{+} + \sqrt{B_{k, q} S_{k} B_{k, q}^{T}} δ_{k, p, 1}^{+} | \geq A L_{q}, ({\hat{δ}}_{k, p}^{+}, ε_{⊥, k, p}) \in Ω_{k, q} ∣ H_{i}, δ_{k, p}^{+}) : \\ P (({\hat{δ}}_{k, p}^{+}, ε_{⊥, k, p}) \notin Ω_{k, q} ∣ H_{0}) \leq P_{f a, i, q} \end{matrix}$ 35

In this form, the objective function no longer depends on $f_{0 : k - 1}$ , but rather on $δ_{k, p}^{+}$ , and only the first component of ${\hat{δ}}_{k, p}^{+}$ affects the position bias of the q-th coordinate state, $d_{k, q}^{+}$ .

3.3 Optimal Detection Region for One Coordinate

Using the theoretical result for an optimization problem similar to Equation (35), derived by Blanch et al. (2017), the optimal innovation-based detection region can be determined. The result of Blanch et al. (2017) states that the first component of the reparametrized random variable capturing the fault (which corresponds to Equation (29) in this paper) becomes the optimal statistic. The solution for Equation (35) is derived as follows:

$Ω_{k, q}^{*} = {γ_{k} | ∣ B_{k, q} γ_{k} ∣ \leq K_{f a, i, q} \sqrt{B_{k, q} S_{k} B_{k, q}^{T}}}$ 36

where $K_{f a, i, q}$ is a multiplier satisfying the false alert rate requirement $P_{f a, i, q}$ . A detailed process for obtaining this solution is described in Appendix B.

This test statistic is easily computed as the magnitude of $B_{k, q} γ_{k}$ , representing the projection of the innovation vector $γ_{k}$ onto the position domain. The corresponding threshold is also simply determined as the product of the matrices of $B_{k, q}$ and $S_{k}$ .

3.4 Equivalence with the SS Statistic

The optimal innovation-based statistic derived in Section 3.3 is a position-domain projection of the innovation (i.e., $B_{k, q} γ_{k}$ ). This subsection details the construction of the snapshot-based SS statistic and establishes its mathematical equivalence to the innovation-based position-domain statistic.

The fault-free estimate at epoch k under the fault hypothesis $H_{i} ({\hat{x}}_{F F, k} \in R^{n_{p}})$ , computed using the current-time normalized measurements, $y_{k}$ , is given by the following:

${\hat{x}}_{F F, k} = {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} y_{k}$ 37

Under the fault hypothesis $H_{i}$ , the estimate error of ${\hat{x}}_{F F, k}$ , denoted as ${\tilde{x}}_{F F, k}$ , is expressed as follows:

${\tilde{x}}_{F F, k} = E_{k}^{T} x_{k} - {\hat{x}}_{F F, k} \sim N (0_{n_{p} \times 1}, {(H_{k, p}^{T} H_{k, p})}^{- 1})$ 38

where $E_{k}^{T} x_{k}$ represents the true partial state vector and $n_{p}$ is the number of partial states. From Equation (37), the test statistic for the q-th coordinate is defined as follows:

$Δ_{k, q} = e_{q}^{T} {\hat{x}}_{k}^{+} - e_{q}^{T} E_{k} {\hat{x}}_{F F, k}$ 39

In Equation (39), $E_{k}$ is introduced to represent $e_{q}^{T}$ as a fixed-size row vector.

The innovation vector decomposes into two components: a function of an a priori estimate and a function of current measurements. Under dynamic sensor faults, as modeled in Equation (11), this fault mode induces a multi-dimensional bias in an a priori estimate through the time-update step. Therefore, a fault-free solution can be constructed from the innovation vector using only the measurement component. Based on the results of Blanch et al. (2017), which established that the SS statistic is optimal for monitoring a single fault hypothesis from an integrity perspective, the snapshot-based SS statistic defined using this fault-free solution is optimal against such faults. We demonstrate that this snapshot-based SS statistic is mathematically equivalent to the position-domain projection of the innovation. A proof of their equivalence is provided in Appendix C, confirming that the same result can be reached by optimizing the detection region in the innovation domain.

4 PERFORMANCE EVALUATION

This section presents a performance evaluation of the proposed test statistic in terms of integrity risk under the dynamic sensor fault scenario. The proposed test statistic is compared with the integrity risk computed using the range-domain monitor proposed by Lee et al. (2023) for a tightly coupled KF-based GPS/IMU navigation system. Section 4.1 details the computation of integrity risk for the IMU sensor fault scenario, and Section 4.2 presents the evaluation results. Because the performance of the proposed monitor depends on the allocation of false alert rate requirements across all directions, Section 4.3 presents a sensitivity analysis of this allocation for the direction of interest.

4.1 Integrity Risk Under IMU Faults

4.1.1 Derivation of the Worst-Case Innovation Bias

This subsection derives the worst-case innovation bias of the innovation-based optimal statistic. This worst-case result will be used to evaluate the integrity risks of both the proposed monitor and the range-domain monitor proposed by Lee et al. (2023) in the following subsection.

From Equations (29), (30), (35), and (36), the integrity risk for the proposed monitor under the $H_{i}$ fault hypothesis can be expressed as follows:

$P (| η_{k, q}^{+} + B_{k, q} b_{k} | \geq A L_{q}, | B_{k, q} γ_{k} | \leq K_{f a, i, q} \sqrt{B_{k, q} S_{k} B_{k, q}^{T}} ∣ H_{i}, b_{k})$ 40

In this subsection, the integrity risk is expressed as a function of the innovation bias $(b_{k})$ , rather than $δ_{k, p}^{+}$ defined in the reparameterization process. This approach facilitates evaluating the performance of the range-domain monitor along with the proposed optimal monitor, because the range-domain monitor requires determining the $b_{k}$ that corresponds to the worst-case $δ_{k, p}^{+}$ for the optimal monitor. If we derive the worst-case $δ_{k, p}^{+}$ rather than the worst-case $b_{k}$ , it is not straightforward to determine the specific $b_{k}$ that corresponds to the worst-case $δ_{k, p}^{+}$ because these terms do not have a one-to-one relationship (see Equations (30) and (32)).

By exploiting the independence between the KF's current innovation and an a posteriori estimate error, as proven by Tanil et al. (2017), Equation (40) is expressed as follows:

$\int_{| η_{k, q}^{+}' | \geq A L_{q}} \int_{| t_{k, q} | \leq K_{f a, i, q} \sqrt{B_{k, q} S_{k} B_{k, q}^{T}}} N_{η_{k, q}^{+}'} (B_{k, q} b_{k}, σ_{k, q}^{2}) \cdot N_{t_{k, q}} (B_{k, q} b_{k}, B_{k, q} S_{k} B_{k, q}^{T}) d t_{k, q} d η_{k, q}^{+}'$ 41

where $η_{k, q}^{+}'$ is $η_{k, q}^{+} + B_{k, q} b_{k}$ and $t_{k, q}$ is $B_{k, q} γ_{k}$ . Based on the approximation made in Equation (B8), the contribution from the opposite side of the mean in the integral is negligible; thus, Equation (41) can be simplified. Without a loss of generality, assuming that $B_{k, q} b_{k}$ is nonnegative, Equation (41) can be approximated as follows:

$\int_{A L_{q}}^{\infty} N_{η_{k, q}^{+}}' (B_{k, q} b_{k}, σ_{k, q}^{2}) d η_{k, q}^{+}' . \int_{- \infty}^{K_{f a, i, q \sqrt{B_{k, q} S_{k} B_{k, q}^{T}}}} N_{t_{k, q}} (B_{k, q} b_{k}, B_{k, q} S_{k} B_{k, q}^{T}) d t_{k, q}$ 42

By normalizing the integrands in Equation (42) with respect to their SD, this equation can be rewritten as follows:

$\int_{\frac{A L_{q}}{σ_{k, q}}}^{\infty} N_{η_{k, q}^{+}}'' (\frac{B_{k, q} b_{k}}{σ_{k, q}}, 1) d η_{k, q}^{+}'' . \int_{- \infty}^{K_{f a, i, q}} N_{t_{k, q}'} (\frac{B_{k, q} b_{k}}{\sqrt{B_{k, q} S_{k} B_{k, q}^{T}}}, 1) d t_{k, q}'$ 43

In Equation (43), the magnitude of $b_{k}$ , which is the mean of the current innovation vector, affects the two integrals in opposite directions. In this context, the ratio $E {[η_{k, q}^{+} ″]}^{2} / E {[t_{k, q}']}^{2}$ , referred to as the failure mode slope $(ρ_{k, q}^{2})$ , is defined and utilized to determine the worst-case innovation bias.

The value of Equation (43) is maximized when $ρ_{k, q}^{2}$ is maximized (Joerger et al., 2014). The worst-case innovation bias $({\bar{b}}_{k})$ that maximizes $ρ_{k, q}^{2}$ is expressed as follows:

$\begin{matrix} {\bar{b}}_{k} & = & \underset{b_{k}}{argmax} ρ_{k, q}^{2} \\ B_{k, q} b_{k} \geq 0 \end{matrix}$ 44

By substituting the means of $η_{k, q}^{+} ″$ and $t_{k, q}'$ given in Equation (43) into $E {[η_{k, q}^{+} ″]}^{2} / E {[t_{k, q}']}^{2}, ρ_{k, q}^{2}$ is expressed as follows:

$ρ_{k, q}^{2} = \frac{b_{k}^{T} B_{k, q}^{T} B_{k, q} b_{k} / σ_{k, q}^{2}}{b_{k}^{T} B_{k, q}^{T} B_{k, q} b_{k} / B_{k, q} S_{k} B_{k, q}^{T}} = \frac{B_{k, q} S_{k} B_{k, q}^{T}}{σ_{k, q}^{2}}$ 45

This result indicates that the integrity risk for the proposed monitor is independent of the bias direction, consistent with the failure mode slope analysis of SS RAIM presented by Joerger et al. (2014).

To evaluate the integrity risk of the innovation-based range-domain monitor proposed by Lee et al. (2023), we assume the worst-case innovation bias along the following direction in Equation (46) (i.e., $S_{k} B_{k, q}^{T}$ ). If the worst-case fault occurs along this direction, the position bias is nonnegative, as assumed in Equation (42), because of the positive definiteness of the covariance matrix $S_{k}$ :

${\bar{b}}_{k} = α_{k} S_{k} B_{k, q}^{T}$ 46

The magnitude of ${\bar{b}}_{k}$ , denoted as $α_{k}$ in Equation (46), is a scalar that maximizes Equation (43) along the worst-case fault direction, and this value is always positive to ensure that ${\bar{b}}_{k}$ is positive. Note that in this paper, we search for the worst-case bias in the innovation domain rather than in the sensor measurement domain. By maximizing the integrity risk over the entire innovation space, which includes the subset of innovation biases reachable from sensor faults, the obtained worst-case bias conservatively bounds all possible cases that could arise from sensor faults.

The proposed monitor defines the test statistic for each coordinate direction (east, north, and up), and triggers an alert if any of the test statistics falls outside its detection region. Thus, $K_{f a, i, q}$ is determined by uniformly allocating the false alert rate requirement to each coordinate.

4.1.2 Integrity Risk Computation

This subsection explains the procedure for computing the integrity risk for the IMU sensor fault scenario. In this analysis, the worst-case IMU fault introducing the worst-case innovation bias $({\bar{b}}_{k})$ , as given in Equation (46), is assumed to occur. Under this fault scenario, the integrity risks are evaluated for both the monitor proposed in this study and the innovation-based range-domain monitor proposed by Lee et al. (2023).

The integrity risk for the q-th coordinate state at epoch k with the positiondomain monitor is evaluated using Equation (42) with the worst-case fault ${\bar{b}}_{k}$ .

The integrity risk for the q-th coordinate state at epoch k, when using the range-domain monitor proposed by Lee et al. (2023), is computed as follows:

$\int_{A L_{q}}^{\infty} N_{η_{k, q}^{+}'} (α_{k} B_{k, q} S_{k} B_{k, q}^{T}, σ_{k, q}^{2}) d η_{k, q}^{+}' \cdot \int_{- T_{k, 1}}^{T_{k, 1}} \dots \int_{- T_{k, mml:m}}^{T_{k, mml:m}} N_{γ_{k}} (α_{k} S_{k} {B_{k, q}}^{T}, S_{k}) d γ_{k}$ 47

where m is the number of measurements (i.e., the number of pseudorange and pseudorange rate measurements used for the KF-based GPS/IMU navigation system) and $T_{k, j}$ is the monitor threshold for the j-th innovation element at epoch k. Because the range-domain monitor defines individual test statistics for all innovation elements, the thresholds are computed by allocating the requirement to each element. For a fair comparison with the proposed monitor, which optimizes integrity performance in a specific direction of interest, we optimally allocate the false alert rate requirement across innovation elements to minimize the protection level (PL), defined as the position error bound satisfying the predefined integrity risk requirement, in the direction of interest, as described by Lee et al. (2023). The MATLAB® function mvncdf is used to integrate the multivariate Gaussian distribution in Equation (47).

It is important to note that ${\bar{b}}_{k}$ does not represent the worst-case fault for the range-domain monitor. Therefore, this evaluation method of using ${\bar{b}}_{k}$ returns a conservative result in terms of verifying the benefit gained by using the proposed monitor compared with the range-domain monitor.

4.2 Performance Analysis Results

The integrity risk for the vertical coordinate is evaluated for simulated autonomous vehicle motion involving two 90° turns, as shown in Figure 1. The assumed KF-based GPS/IMU navigation system utilizes pseudorange and pseudorange rate measurements from GPS signals for the measurement update and angular rate and acceleration measurements from a tactical-grade IMU for the time update. The GPS constellation is simulated based on the 27-expandable GPS constellation given in Section 3.2 from the Department of Defense (2020).

FIGURE 1

Ground truth vehicle trajectory

In the measurement-update step, the ionosphere-free (IFree) combination based on GPS L1 and L5 measurements is employed to eliminate a first-order ionospheric delay, and the remaining residual ionospheric delay is assumed to be negligible. Further specifications for GPS and IMU measurement errors are detailed in Tables 1 and 2, respectively. For the multipath and receiver noise, an inflation factor of 2.6 is applied to account for the increased level of multipath and receiver noise in the IFree combination. The SD of the IFree measurement error is calculated as the root sum square of the SDs of the individual error sources listed in Table 1. The prior probability of IMU faults in Table 2 is roughly estimated from the mean time between failures reported by EMCORE (n.d.).

View this table:

Table 1 GPS Measurement Specifications

View this table:

Table 2 IMU Measurement Specifications

The false alert rate requirement is set at $10^{- 6} / h$ , as used in civil aviation enroute operations (International Civil Aviation Organization, 1996). This requirement is equally allocated between GPS and IMU fault monitors.

To compare the integrity risks for the proposed optimal position-domain monitor and the range-domain monitor, the integrity risks are computed with respect to a vertical AL (VAL) ranging from 10 to 70 m . In this evaluation, the IMU fault is assumed to occur after the filter has sufficiently converged.

Figure 2 shows the integrity risk results obtained when using the proposed innovation-based optimal monitor and the innovation-based range-domain monitor proposed by Lee et al. (2023) with optimal continuity allocation minimizing the vertical PL (VPL). The integrity risks were computed 60 s after the simulation initiation. Eight GPS satellites were available during this period. In Figure 2, the solid curve represents the integrity risk for the optimal monitor, and the dotted curve represents the integrity risk for the range-domain monitor. For the entire range of VAL, the integrity risk of the optimal monitor is lower than that of the innovation-based range-domain monitor, indicating that the proposed monitor can ensure higher availability than the monitor proposed by Lee et al. (2023). If the range-domain monitor is evaluated for the worst-case fault of itself instead of the optimal monitor's worst-case fault, its integrity risk for a given VAL would be even higher than this result. Conversely, the integrity risk of the optimal monitor would be lower than the result depicted in Figure 2. Thus, the results shown in Figure 2 represent a conservative benefit gained by using the optimal monitor. The integrity risks for both monitors decrease as the VAL increases. The position-domain monitor is optimal (i.e., minimizing integrity risk) against the faults assumed in the simulation compared with other innovation-based monitors. Consequently, as the worst-case fault magnitude increases along with increasing VAL, the integrity risk of the proposed monitor shows a more rapid decrease than that of the innovation-based range-domain monitor. For an integrity risk requirement of $10^{- 7}$ , service with a VAL of 28.9 m is available with the optimal position-domain monitor, whereas a VAL of 48.3 m is required for the range-domain monitor to ensure availability.

FIGURE 2

Simulated integrity risk results for a VAL ranging from 10 to 70 m under the worst-case IMU fault, with a zoomed-in view for a VAL range of 10 to 20 m shown in the inset plot

Figure 3 shows the ratio of the integrity risk of the range-domain monitor to that of the optimal monitor. For a VAL below 15 m, the ratio approaches unity, as the integrity risks in both cases approach the prior probability of an IMU fault (i.e., the conditional integrity risks become one) under the assumed simulation conditions. However, for a VAL of 35 m, which is the requirement for the localizer performance with vertical guidance (LPV)-200, the integrity risk of the proposed monitor is lower by approximately three orders of magnitude or more than that of the innovation-based range-domain monitor.

FIGURE 3

Ratio of the integrity risk of the innovation-based range-domain monitor to that of the innovation-based position-domain (optimal) monitor, with a zoomed-in view for a VAL range of 10 to 20 m shown in the inset plot

4.3 Sensitivity Analysis of False Alert Rate Requirement Allocation

In the previous subsection, the integrity risk of the proposed monitor was evaluated under a uniform false alert rate (or continuity risk) requirement allocation among all three directions. However, the performance of the proposed monitor (and consequently its benefit over the range-domain monitor) depends strongly on this allocation. Therefore, this subsection conducts a sensitivity analysis by varying the continuity risk budget allocated to the vertical direction from 1% to 99% of the total continuity risk requirement, covering operationally extreme allocation scenarios. The performance is compared against the range-domain monitor with optimal continuity allocation minimizing the VPL. All other simulation conditions are identical to those in Section 4.2.

Figure 4 shows the results of the sensitivity analysis of integrity risk (log-scale) with respect to continuity allocation to the vertical direction. Subplot (a) shows results for the position-domain monitor, and subplot (b) shows results for the range-domain monitor. Black solid curves represent integrity risk contours at $10^{- 7}$ and $10^{- 9}$ . As the allocated continuity risk budget increases, the monitor threshold decreases, which reduces the worst-case fault magnitude of the proposed monitor and thereby the integrity risk. The multiplier corresponding to the allocated continuity budget (i.e., $K_{f a, i, q}$ ) exhibits a sharper increase as the allocated budget decreases. Consequently, in Figure 4(a), the integrity risk contours show more noticeable changes in low-continuity-allocation regions. The proposed position-domain monitor is optimal in the sense that it minimizes the integrity risk in the presence of its worst-case fault, making it the most sensitive monitor against the assumed fault in the simulation. In contrast, the range-domain monitor is not as sensitive as the proposed monitor. Therefore, there are no notable changes in its integrity risk results across different continuity allocations, as shown in Figure 4(b).

FIGURE 4

Sensitivity analysis of integrity risk (log-scale) with respect to continuity risk allocation

Subplots (a) and (b) show results for the position-domain monitor and range-domain monitor, respectively. Black solid curves represent integrity risk contours at 10^–7 and 10^–9.

Figure 5 shows the integrity risk ratio of the range-domain monitor to the optimal position-domain monitor for different continuity allocations. Black solid curves represent ratio contours at $10^{1}, 10^{2}, 10^{5}, 10^{10}$ , and $10^{20}$ . Because the integrity risk of the range-domain monitor remains nearly constant across all continuity allocations, the ratio results exhibit the same trend as the optimal monitor's integrity risk shown in Figure 4(a). Even when the continuity budget allocated to the vertical direction is as low as 1%, the integrity risk of the proposed monitor remains lower than that of the range-domain monitor by two orders of magnitude or more for a VAL of 35 m.

FIGURE 5

Ratio of integrity risk (log-scale) for the range-domain monitor to the position-domain (optimal) monitor with respect to continuity risk allocation

Black solid curves represent ratio contours at 10¹, 10², 10⁵, 10¹⁰, and 10²⁰.

5 PERFORMANCE AND COMPUTATIONAL LOAD COMPARISON WITH THE FULL SS MONITOR

This section compares the integrity performance and computational load of the snapshot-based and full SS IM methods to identify when the proposed approach offers practical advantages. The performance difference between the two IM methods varies with the accuracy of the fault-free solution, which depends on the navigation system configuration. Thus, a comparison is performed across three sensor configurations: GPS/IMU, dual-constellation GNSS (GPS and Galileo)/IMU, and GNSS/altimeter/IMU. The GPS/IMU configuration represents the baseline with a single GNSS constellation. The dual-constellation GNSS/IMU configuration adds a Galileo constellation to examine whether increased measurement redundancy narrows the performance gap. However, because GNSS satellite geometry fundamentally limits vertical observability regardless of the number of constellations, direct altitude measurements are required to provide observability across all coordinate directions. Therefore, the GNSS/altimeter/IMU configuration is included to assess the performance difference when sufficient observability is available in all directions. The integrity performance of each IM method is compared in terms of the PL of each method against IMU faults.

5.1 PL Computation

Under the IMU fault hypothesis, the following upper bound of the integrity risk for the q-th coordinate proposed by Blanch et al. (2015) can be used to compute the PL:

$P (| {\tilde{x}}_{F F, k, q} | > A L_{q} - T_{k, q} ∣ H_{I M U}) P_{H_{I M U}}$ 48

where ${\tilde{x}}_{F F, k, q}$ denotes the error of the fault-free q-th coordinate position estimate at epoch k. For the full SS monitor, ${\tilde{x}}_{F F, k, q}$ is an a posteriori estimate error of the KF without IMU measurements. For the proposed snapshot-based SS monitor, ${\tilde{x}}_{F F, k, q}$ is the snapshot solution error computed using only the current-time measurements. $T_{k, q}$ is the monitoring threshold of the SS monitor and is computed using the key property of the SS statistic as follows (Blanch et al., 2015):

$T_{k, q} = K_{f a, I M U, q} \cdot \sqrt{σ_{F F, k, q}^{2} - σ_{k, q}^{2}}$ 49

Here, $K_{f a, I M U, q}$ is a multiplier satisfying the false alert rate requirement allocated to the IMU fault monitor and the q-th coordinate direction, and $σ_{F F, k, q}$ is the SD of ${\tilde{x}}_{F F, k, q}$ . This property was demonstrated for the full SS monitor implemented in KF-based navigation systems by Gunning et al. (2018). This property also holds for the snapshot-based SS monitor because the snapshot-based SS monitor is formulated using stacked measurements. When an a priori estimate and current measurements are formulated as stacked measurements, their WLS solution becomes mathematically identical to an a posteriori estimate. Therefore, the snapshot-based SS statistic, i.e., the difference between the snapshot solution (from current measurements only) and an a posteriori estimate, is the direct result of applying the SS concept to this stacked measurement formulation against dynamic sensor faults that induce bias in an a priori estimate.

Because ${\tilde{x}}_{F F, k, q}$ follows a zero-mean Gaussian distribution under $H_{I M U}$ , the PL for the q-th coordinate satisfying the integrity risk requirement $P_{H M I, H_{I M U, q}}$ is computed as follows:

$P L_{k, q} = T_{k, q} + K_{m d, I M U, q} \cdot σ_{F F, k, q}$ 50

where $K_{m d, I M U, q}$ is a multiplier satisfying the conditional integrity requirement under the IMU fault condition, computed as $- Q^{- 1} (P_{HMI, H_{I M U, q}} / P_{H_{I M U}}) . Q^{- 1} (■)$ denotes the inverse cumulative distribution function (CDF) of the standard normal distribution, and $P_{H_{IMU}}$ is the prior probability of IMU faults.

5.2 Comparison with the Full SS Across Sensor Configurations

5.2.1 Simulation Conditions

The dual-constellation GNSS/IMU configuration employs measurements from both the 27-satellite GPS constellation used in Section 4 and an additional 27-satellite Galileo constellation transmitting E1 and E5a signals. For the GNSS/ altimeter/IMU configuration, barometric altitude measurements with an SD of 2.08 m (Lee et al., 2020) are further integrated with the GNSS/IMU system.

The integrity risk requirement is set to $10^{- 7} / h$ , as suggested by Whitty and Walport (2018) for autonomous vehicles. One-third of the total integrity risk requirement is allocated to the IMU fault hypothesis. The allocated integrity risk budget is further uniformly distributed among the three coordinate directions (east, north, and up). Other simulation conditions are identical to those used in Section 4.

The performance and computational load are compared using MATLAB on a system with an Intel Core i7-14700K processor and 32 GB RAM. Computational execution times are derived from 1,000 independent 60-s (or 120 measurement-update epoch) simulations of autonomous vehicle motion identical to that of Section 4 after 10 warm-up runs to account for just-in-time compilation effects, which consume time to translate interpreted MATLAB code into machine code at runtime. For each simulation, the cumulative IM execution time across all epochs is measured. The median of these 1,000 cumulative times is then divided by 120 epochs to obtain the reported per-epoch execution time. Each IM execution includes test statistic computation, fault detection, and PL computation, timed in isolation from the main navigation filter execution.

5.2.2 Simulation Results

Figures 6 and 7 present the PLs across three sensor configurations: GPS/IMU, dual-constellation GNSS/IMU, and GNSS/altimeter/IMU. Figure 6 shows the PLs over 60 s (120 measurement-update epochs) for each coordinate direction, while Figure 7 compares the converged PLs (i.e., values at 60 s). The PLs in Figure 6 initially increase and then gradually converge. This trend arises from changes in the variance difference between the fault-free solution and an a posteriori estimate using all available information, which is used to compute the threshold term in the PL formulation (Equation (50)). This term is dominant because $K_{fa,IMU, q}$ is larger than $K_{m d, I M U, q}$ owing to the prior probability of IMU faults.

FIGURE 6

Results of PLs across sensor configurations

The top, middle, and bottom plots show the PLs in the east, north, and vertical (up) directions, respectively. Red solid curves with circle markers represent the GPS/IMU configuration with the snapshot-based SS IM (S-SS), yellow solid curves with diamond markers represent the GPS/IMU configuration with the full SS IM (F-SS), magenta solid curves with square markers represent the GNSS/IMU configuration with the S-SS, cyan solid curves with asterisk markers represent the GNSS/IMU configuration with the F-SS, green solid curves with star markers represent the GNSS/altimeter/IMU configuration with the S-SS, and blue solid curves with hexagram markers represent the GNSS/altimeter/IMU configuration with the F-SS.

FIGURE 7

Converged PLs across sensor configurations

Subplots (a)–(c) show results for the GPS/IMU configuration, GNSS/IMU configuration, and GNSS/altimeter/IMU configuration, respectively. In each subplot, the left three bars present results for the snapshot-based SS IM (S-SS), and the right three bars present results for the full SS IM (F-SS). Red, green, and blue bars represent PLs in the east, north, and vertical (up) directions, respectively.

The performance comparison reveals two distinct trends depending on the coordinate direction. In the horizontal directions (east and north), the performance gap between the snapshot-based and full SS IM methods diminishes rapidly as measurement redundancy increases. As shown in Figure 7, the gap starts at 10.9% (2.2 m) and 2.5% (0.5 ms) in the east and north directions when using GPS measurements only, becomes negligible with dual-constellation measurements (0.7%–2.0%, 0.1–0.3 m), and approaches zero with the GNSS/altimeter/IMU configuration (≤1.5%, ≤0.2 m). This result demonstrates that measurement redundancy from multiple GNSS constellations enables the snapshot solution to achieve an accuracy comparable to that of the filter solution without IMU measurements in the horizontal direction.

In contrast, the vertical direction exhibits different behavior, as shown in the bottom subplots of Figures 6 and 7. The performance gap persists when dual-constellation measurements are used. This result occurs because, although dual-constellation measurements improve horizontal geometry, vertical observability remains fundamentally limited by GNSS satellite geometry. However, the gap reduces significantly (to 24.5% or 3.7 m) when direct altitude measurements from the barometric sensor are introduced, providing direct vertical observability.

In terms of computational efficiency, the snapshot-based monitor is superior to the full SS monitor. In the GPS/IMU configuration, the full SS monitor requires approximately 9.6 times the execution time of the snapshot-based monitor (7.3 μs for the snapshot-based SS monitor and 70.4 μs for the full SS monitor). This substantial difference reflects the additional computational burden of maintaining and processing the fault-free filter solution required for the full SS monitor. The computational advantage of the snapshot-based approach exhibits consistent growth across sensor configurations as system complexity increases. The execution time ratio of the full SS IM to the snapshot-based SS IM increases from 9.6 (GPS/IMU) to 11.5 (GNSS/IMU) to 11.7 (GNSS/Altimeter/IMU), and the absolute time difference grows from 63.1 to 108.3 to 112.7 μs.

5.2.3 Discussion

The results presented above illustrate several observations for monitoring dynamic sensor faults. As measurement redundancy and diversity increase through multi-sensor integration, the accuracy of the snapshot solution improves, narrowing the performance gap with the full SS approach. The horizontal performance gap becomes negligible with dual-constellation measurements, demonstrating that measurement redundancy from multiple GNSS constellations enables the snapshot solution to achieve an accuracy comparable to that of the filter solution without dynamic sensor measurements in the horizontal direction. The vertical performance gap persists with dual-constellation GNSS owing to the limited vertical observability of satellite geometry, but decreases significantly when direct altitude measurements are introduced. These observations suggest that for a multi-sensor configuration that provides sufficient observability in all coordinate directions, the snapshot-based SS IM can achieve integrity performance comparable to that of the full SS IM.

Regarding computational efficiency, the execution time ratio of the full SS approach to the snapshot-based SS approach increases as the complexity of the navigation system configuration increases, reflecting the fundamental difference between single-filter and two-filter architectures. As the measurement vector size grows, both the absolute time difference and the execution time ratio increase, indicating that the computational advantage of the snapshot-based approach grows with system complexity, as the computational cost of the full SS IM method scales more significantly.

These observations suggest that the snapshot-based SS approach may be more beneficial in richer multi-sensor configurations. When measurement redundancy provides sufficient observability in all coordinate directions, the snapshot solution can achieve an accuracy comparable to that of the filter solution without dynamic sensor measurements. Consequently, the approach offers comparable integrity performance while providing a computational advantage.

6 CONCLUSION

This paper proposed an optimal test statistic for an innovation-based monitor that relies solely on a single filter to detect sensor faults occurring in the KF time-update step, thereby ensuring low computational complexity. Although the SS concept provides the optimal solution for designing the optimal innovation-based statistic under dynamic sensor faults, corresponding to the snapshot-based SS statistic, this paper approaches the problem from the perspective of innovation-domain optimization rather than applying the SS concept. We formulated the test statistic design as a min-max optimization problem and reparametrized it into a form in which the test statistic is expressed as two independent random variables: one variable capturing the bias in an a priori estimate and another variable accounting solely for innovation noise. We then applied the Neyman–Pearson lemma to derive the optimal statistic, defined as the projection of the innovation vector onto the position domain. We demon-strated that this independently derived solution is mathematically equivalent to the snapshot-based SS statistic, thereby confirming that the same result can be reached through innovation-domain optimization.

The proposed monitor was evaluated using a tightly coupled KF-based GPS/ IMU navigation system under worst-case IMU fault scenarios. Sensitivity analysis of integrity risk across VALs (10–70 m) demonstrated that the proposed position-domain monitor achieves significantly lower integrity risk than the innovation-based range-domain monitor proposed by Lee et al. (2023).

Additionally, to identify when the proposed IM method offers practical advantages over the full SS IM using parallel filters, we compared the integrity performance and computational load of both approaches across three sensor configurations (GPS/IMU, dual-constellation GNSS/IMU, and GNSS/altimeter/IMU). As measurement redundancy and diversity increase through multi-sensor integration, the accuracy of the snapshot solution improves, narrowing the performance gap with the full SS approach. The horizontal performance gap becomes negligible with dual-constellation measurements, whereas the vertical gap decreases significantly with direct altitude measurements. Regarding computational efficiency, the snapshot-based approach consistently provides lower execution times, with the absolute difference and ratio improving as system complexity grows. These observations suggest that the snapshot-based approach (i.e., the innovation-based position-domain IM) is particularly well-suited for multi-sensor configurations in which measurement redundancy is sufficient to enable the snapshot solution to achieve an accuracy comparable to that of the filter solution without dynamic sensor measurements.

Future work will focus on extending the proposed monitor to a decentralized filter-based multi-sensor navigation system, in which the state is estimated using multiple parallel filters, each utilizing distinct sensor subgroups. The integrity performance of multi-sensor systems can be improved with this decentralized framework, when the proposed monitor is applied individually to each subgroup of multiple dynamic sensors.

HOW TO CITE THIS ARTICLE

Nam, G., Min, D., Kim, N.M., & Lee, J. (2026). Kalman filter innovation-based optimal integrity monitoring against unknown bias in time updates. NAVIGATION, 73. https://doi.org/10.33012/navi.769

APPENDIX | A

This appendix proves the statistical independence of the bias estimate in an a priori estimate computed from the innovation vector, denoted as ${\hat{d}}_{k, p}^{-}$ , with its residual vector $ε_{⊥, k, p}$ . Because the innovation vector follows a Gaussian distribution, a lack of correlation between ${\hat{d}}_{k, p}^{-}$ and $ε_{⊥, k, p}$ implies their statistical independence.

From Equations (27) and (28), the correlation between ${\hat{d}}_{k, p}^{-}$ and $ε_{⊥, k, p}$ is given by the following:

$\begin{matrix} E [{(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} γ_{k} γ_{k}^{T} (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T})] \\ = {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} E [γ_{k} γ_{k}^{T}] (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) \end{matrix}$ A1

By substituting the definition of the covariance matrix of the innovation vector (i.e., Equation (16)) into Equation (A1), the correlation is given as follows:

${(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} (H_{k} P_{k}^{-} H_{k}^{T} + I_{m}) (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T})$ A2

To simplify Equation (A2), we re-order the KF state vector as follows:

$x_{k} = [\begin{matrix} x_{k, p} \\ x_{k, 0} \end{matrix}]$ A3

Here, $x_{k, p} \in R^{n_{p}}$ denotes the partial state vector whose information is provided by sensors used in KF measurement-update steps, and $x_{k, 0} \in R^{n_{0}}$ indicates the remaining state components. Then, the partial measurement matrix can be expressed as follows:

$H_{k, p} = [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}]$ A4

where m denotes the number of measurements.

By substituting Equation (A4) into Equation (A2), the correlation is expressed as follows:

$\begin{matrix} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} (H_{k, p} P_{k, p}^{-} H_{k, p}^{T} + I_{m}) (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) \\ = (P_{k, p}^{-} H_{k}^{T} + {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) \\ = P_{k, p}^{-} H_{k}^{T} - P_{k, p}^{-} H_{k}^{T} + {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} - {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} = 0 \end{matrix}$ A5

where $P_{k, p}^{-}$ is the partial matrix of $P_{k}^{-}$ , corresponding to the partial state.

Thus, ${\hat{d}}_{k, p}^{-}$ is statistically independent of $ε_{⊥, k, p}$ , and the innovation is expressed as the sum of these two independent random vectors, as shown below:

$γ_{k} = H_{k, p} {\hat{d}}_{k, p}^{-} + ε_{⊥, k, p}$ A6

APPENDIX | B

In this section, we apply the results of Blanch et al. (2017) to our problem. One of the key results from Blanch et al. (2017) is that a solution to an optimization problem, such as Equation (35), can be obtained from the following min-max problem:

$\begin{matrix} min_{Ω} max_{δ_{p}^{+}} \frac{1}{2} (F (δ_{p}^{+}, Ω) + F (- δ_{p}^{+}, Ω)) : \\ P (({\hat{δ}}_{p}^{+}, ε_{⊥, p}) \notin Ω ∣ H_{0}) \leq P_{f a} \end{matrix}$ B1

where $F (δ_{p}^{+}, Ω)$ is the integrity risk for bias $δ_{p}^{+}$ and detection region Ω under the fault hypothesis. For simplicity, we drop the unnecessary subscripts q, i, and k.

To apply this result, the following condition must be satisfied in our problem (Blanch et al., 2017): the optimal detection region of the objective function in Equation (B1) for a given bias $δ_{p}^{+} (Ω_{δ_{p}^{+}})$ is symmetric. This condition is shown to be satisfied in the following subsection.

Neyman–Pearson Lemma

For a given bias $δ_{p}^{+}$ , the corresponding optimal detection region $Ω_{δ_{p}^{+}}$ is derived using the Neyman–Pearson lemma (Van Trees, 1968), which provides the solution of the following constrained optimization problem over the integral region Ω with a constant C:

$\begin{matrix} min_{Ω} \int_{Ω} f (y) d y : \\ \int_{Ω} g (y) d y \geq 1 - C \end{matrix}$ B2

From the Neyman–Pearson lemma, the optimal region $Ω_{Opt}$ is given as follows:

$Ω_{O p t} = {y | \frac{f (y)}{g (y)} \leq T}$ B3

The threshold T is determined from the constraint as shown below:

$\int_{Ω_{O p t}} g (y) d y = 1 - C$ B4

Optimal Detection Region for a Bias $δ_{p}^{+}$

The lack of correlation between the KF's current innovation and its a posteriori estimate error proven by Tanil et al. (2017), along with the Gaussian assumption, ensures the statistical independence of two events used to define the integrity risk. Thus, $F (δ_{p}^{+}, Ω)$ can be expressed as follows:

$P (| η^{+} + \sqrt{B S B^{T}} δ_{p, 1}^{+} | \geq A L ∣ H, δ_{p}^{+}) P ((ε_{⊥, p}, {\hat{δ}}_{p}^{+}) \in Ω ∣ H, δ_{p}^{+})$ B5

To apply the Neyman–Pearson lemma, we define the function f as follows, using a similar formulation introduced by Blanch et al. (2017):

$\begin{aligned} f (ε_{⊥, p}, {\hat{δ}}_{p}^{+}) = & \frac{1}{2} (P (| η^{+} + \sqrt{B S B^{T}} δ_{p, 1}^{+} | \geq A L ∣ H, δ_{p}^{+}) p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, δ_{p}^{+}) \\ + P (| η^{+} - \sqrt{B S B^{T}} δ_{p, 1}^{+} | \geq A L ∣ H, - δ_{p}^{+}) p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, - δ_{p}^{+})) \end{aligned}$ B6

where $p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, δ_{p}^{+})$ denotes the joint probability density function (PDF) of $ε_{⊥, p}$ and ${\hat{δ}}_{p}^{+}$ under the H fault hypothesis. Because $ε_{⊥, p}$ and ${\hat{δ}}_{p}^{+}$ are statistically independent, as proven in Appendix A, this joint PDF can be expressed as follows:

$p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, δ_{p}^{+}) = N_{ε_{⊥, p}} (0_{m - n_{p} \times 1}, P_{⊥, p} S P_{⊥, p}^{T}) \cdot N_{{\hat{δ}}_{p}^{+}} (δ_{p}^{+}, I_{n_{p}})$ B7

where the notation of $N_{x} (μ, Σ)$ represents the PDF of x that follows $N (μ, Σ)$ . The covariance matrix of ${\hat{δ}}_{p}^{+}$ is determined from the conditions described in Equation (33).

By exploiting the fact that the contribution from the opposite side of the position mean in the integral is negligible (i.e., the contribution of the region $η^{+} + \sqrt{B S B^{T}} δ_{p, 1}^{+} \leq - A L$ is negligible when $δ_{p, 1}^{+}$ is positive), the approximation used by Blanch et al. (2017) can be applied to simplify this function. Without a loss of generality, the case in which $δ_{p, 1}^{+}$ is positive is considered. Assuming that $δ_{p, 1}^{+}$ is positive, Equation (B6) can be simplified as follows:

$\begin{matrix} \frac{1}{2} {(1 - Q (\frac{A L}{σ} - \frac{\sqrt{B S B^{T}} δ_{p, 1}^{+}}{σ})) p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, δ_{p}^{+}) \\ + Q (- \frac{A L}{σ} + \frac{\sqrt{B S B^{T}} δ_{p, 1}^{+}}{σ}) p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, - δ_{p}^{+})} \end{matrix}$ B8

where $Q (■)$ represents the CDF of the standard normal distribution and σ is the SD of $η^{+}$ . Because Q is a point-symmetric CDF about (0,1/2), Equation (B8) can also be expressed as follows:

$\frac{1}{2} Q (- \frac{A L}{σ} + \frac{\sqrt{B S B^{T}} δ_{p, 1}^{+}}{σ}) (p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, δ_{p}^{+}) + p (ε_{⊥, p}, {\hat{δ}}_{p}^{+} ∣ H, - δ_{p}^{+}))$ B9

Based on Equations (B7) and (B9), the optimization problem for determining the optimal detection region corresponding to a given bias ${\hat{δ}}_{p}^{+}$ can be formulated as follows:

$\begin{array}{r} min_{Ω} \int_{(ε_{⊥, p}, {\hat{δ}}_{p}^{+}) \in Ω} N_{ε_{⊥, p}} (0_{m - n_{p} \times 1}, P_{⊥, p} S P_{⊥, p}^{T}) \cdot (N_{{\hat{δ}}_{p}^{+}} (δ_{p}^{+}, I_{n_{p}}) + N_{{\hat{δ}}_{p}^{+}} (- δ_{p}^{+}, I_{n_{p}})) d ε_{⊥, p} d {\hat{δ}}_{p}^{+} : \\ \int_{(ε_{⊥, p}, {\hat{δ}}_{p}^{+}) \in Ω} N_{ε_{⊥, p}} (0_{m - n_{p} \times 1}, P_{⊥, p} S P_{⊥, p}^{T}) \cdot N_{{\hat{δ}}_{p}^{+}} (0_{n_{p} \times 1}, I_{n_{p}}) d ε_{⊥, p} d {\hat{δ}}_{p}^{+} \geq 1 - P_{f a} \end{array}$ B10

From the Neyman–Pearson lemma, the optimal detection region $Ω_{δ_{p}^{+}}$ is computed as follows:

$Ω_{δ_{p}^{+}} = {{\hat{δ}}_{p}^{+}, ε_{⊥, p} | \exp (- \frac{1}{2} δ_{p}^{+ T} δ_{p}^{+}) \cosh (δ_{p}^{+ T} {\hat{δ}}_{p}^{+}) \leq T}$ B11

where $\cosh (■)$ represents the hyperbolic cosine function and T is chosen to meet the false alert requirement. Because only the first component of $δ_{p}^{+}$ (and consequently ${\hat{δ}}_{p}^{+}$ ) affects the integrity risk, as proven by Blanch et al. (2017) for an optimization analogous to Equation (35), Equation (B11) can be simplified as follows:

$Ω_{δ_{p}^{+}} = {{\hat{δ}}_{p}^{+}, ε_{⊥, p} ∣ \cosh (δ_{p, 1}^{+} {\hat{δ}}_{p, 1}^{+}) \leq T^{'}}$ B12

Symmetry of $Ω_{δ_{p}^{+}}$

$Ω_{δ_{p}^{+}}$ is trivially symmetric because the hyperbolic cosine function is a symmetric function.

Thus, based on the result of Blanch et al. (2017), it can be concluded that Equation (B12) is the optimal solution of the original optimization problem (i.e., Equation (35)).

Further Simplification of $Ω_{δ_{p}^{+}}$

Given that the hyperbolic cosine function is symmetric and monotonically increasing for positive arguments, Equation (B12) can be further simplified as follows (Blanch et al., 2017):

$Ω_{δ_{p}^{+}} = {{\hat{δ}}_{p}^{+}, ε_{⊥, p} | | {\hat{δ}}_{p, 1}^{+} ∣ \leq K_{f a}}$ B13

Because ${\hat{δ}}_{p, 1}^{+}$ follows a standard normal distribution under the nominal condition, $K_{f a}$ is computed as $- Q^{- 1} (P_{f a} / 2)$ , where $Q^{- 1} (■)$ denotes the inverse function of $Q (■)$ . Consequently, by substituting the relationship between ${\hat{δ}}_{p, 1}^{+}$ and the innovation vector, Equation (36) is obtained.

APPENDIX | C

This appendix proves the equivalence between the snapshot-solution-based SS statistic and the innovation-based position-domain statistic. To show this equivalence, we first express each statistic in terms of error components. The two statistics are given by the following:

$\begin{matrix} B_{k, q} γ_{k} = & - e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} y_{k} \\ + e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} {\hat{x}}_{k}^{-} \end{matrix}$ C1

$Δ_{k, q} = e_{q}^{T} {\hat{x}}_{k}^{-} + e_{q}^{T} K_{k} (y_{k} - H_{k} {\hat{x}}_{k}^{-}) - e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} y_{k}$ C2

By adding and subtracting the term $e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} x_{k}$ in Equation (C1), $B_{k, q} γ_{k}$ can be written as follows:

$- e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k} + e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} {\tilde{x}}_{k}^{-}$ C3

By adding and subtracting the term $e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} x_{k} - e_{q}^{T} K_{k} H_{k} x_{k}$ in Equation (C2), $Δ_{k, q}$ can be expressed as follows:

$- e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} x_{k} + e_{q}^{T} {\hat{x}}_{k}^{-} + e_{q}^{T} K_{k} v_{k} - e_{q}^{T} K_{k} H_{k} {\tilde{x}}_{k}^{-} - e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k}$ C4

To express $Δ_{k, q}$ solely in terms of error components, we re-order the KF state vector as follows:

$x_{k} = [\begin{matrix} x_{k, p} \\ x_{k, 0} \end{matrix}]$ C5

where $x_{k, p} \in R^{n_{p}}$ denotes the partial state vector whose information is provided by sensors used in KF measurement-update steps and $x_{k, 0} \in R^{n_{0}}$ indicates the remaining state components. Then, the partial measurement matrix can be expressed as follows:

$H_{k, p} = [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}] [\begin{matrix} I_{n_{p}} \\ 0_{n_{0} \times n_{p}} \end{matrix}] = H_{k} E_{k}$ C6

where m denotes the number of measurements.

By substituting Equation (C6) into the first term of Equation (C4), Equation (C4) can be simplified as follows:

$\begin{matrix} - e_{q}^{T} [\begin{matrix} I_{n_{p}} \\ 0_{n_{0} \times n_{p}} \end{matrix}] {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}] x_{k} = - e_{q}^{T} [\begin{matrix} I_{n_{p}} \\ 0_{n_{0} \times n_{p}} \end{matrix}] [\begin{matrix} I_{n_{p}} & 0_{n_{p} \times n_{0}} \end{matrix}] x_{k} \\ = - e_{q}^{T} [\begin{array}{cc} I_{n_{p}} & 0_{n_{p} \times n_{0}} \\ 0_{n_{0} \times n_{p}} & 0_{n_{0} \times n_{0}} \end{array}] x_{k} = - e_{q}^{T} x_{k} \end{matrix}$ C7

Thus, Equation (C2) is given by the following:

$Δ_{k, q} = e_{q}^{T} (I_{n} - K_{k} H_{k}) {\tilde{x}}_{k}^{-} - e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k} + e_{q}^{T} K_{k} v_{k}$ C8

Now, both test statistics consist of two error components: (1) the contribution of the current-time measurement noise ( $v_{k}$ ) and (2) the contribution of an a priori estimate error $({\tilde{x}}_{k}^{-})$ . In the following subsections, a term-by-term comparison of the two test statistics is conducted to prove the equivalence of these two statistics.

Equivalence of Current-Time Measurement Noise Components

The measurement noise components of the two statistics are as follows:

${(B_{k, q} γ_{k})}_{v_{k}} = - e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k} + e_{q}^{T} K_{k} H_{k} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k}$ C9

${(Δ_{k, q})}_{v_{k}} = - e_{q}^{T} E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k} + e_{q}^{T} K_{k} v_{k}$ C10

Because the first terms of each test statistic are identical, we focus on the second terms of Equations (C9) and (C10). The difference between the second terms of Equations (C9) and (C10) is given by the following:

$e_{q}^{T} K_{k} H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} v_{k} - e_{q}^{T} K_{k} v_{k} = - e_{q}^{T} K_{k} (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) v_{k}$ C11

By applying the matrix inversion lemma (Guttman, 1946) to the Kalman gain $K_{k}$ , Equation (C11) is expressed as follows:

$- e_{q}^{T} P_{k}^{-} (I_{m} - H_{k} {(P_{k}^{- - 1} + H_{k}^{T} H_{k})}^{- 1}) H_{k}^{T} (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) v_{k}$ C12

From Equation (C6), the term $H_{k}^{T} (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T})$ in Equation (C12) can be written as follows:

$[\begin{matrix} H_{k, p}^{T} \\ 0_{n_{0} \times m} \end{matrix}] (I_{m} - H_{k, p} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T}) = [\begin{matrix} 0_{n_{p} \times m} \\ 0_{n_{0} \times m} \end{matrix}]$ C13

Thus, Equation (C11), which represents the difference between Equations (C9) and (C10), becomes zero.

Equivalence of a Priori Estimate Error Components

The a priori estimate error components of the two statistics are as follows:

${(B_{k, q} γ_{k})}_{{\tilde{x}}_{k}^{-}} = e_{q}^{T} (I_{n} - K_{k} H_{k}) E_{k} {(H_{k, p}^{T} H_{k, p})}^{- 1} H_{k, p}^{T} H_{k} {\tilde{x}}_{k}^{-}$ C14

${(Δ_{k, q})}_{{\tilde{x}}_{\bar{k}}} = e_{q}^{T} (I_{n} - K_{k} H_{k}) {\tilde{x}}_{k}^{-}$ C15

From Equation (C7), Equation (C14) is given by the following:

${(B_{k, q} γ_{k})}_{{\tilde{x}}_{\bar{k}}^{-}} = e_{q}^{T} (I_{n} - K_{k} H_{k}) [\begin{array}{cc} I_{n_{p}} & 0_{n_{p} \times n_{0}} \\ 0_{n_{0} \times n_{p}} & 0_{n_{0} \times n_{0}} \end{array}] {\tilde{x}}_{k}^{-}$ C16

The Kalman gain can also be expressed as follows:

$K_{k} = [\begin{matrix} K_{k, p} \\ K_{k, 0} \end{matrix}]$ C17

By substituting Equations (C6) and (C17) into Equation (C16), ${(B_{k, q} γ_{k})}_{{\tilde{x}}_{\bar{k}}^{-}}$ is expressed as follows:

$\begin{matrix} e_{q}^{T} (I_{n} - [\begin{matrix} K_{k, p} \\ K_{k, 0} \end{matrix}] [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}]) [\begin{array}{cc} I_{n_{p}} & 0_{n_{p} \times n_{0}} \\ 0_{n_{0} \times n_{p}} & 0_{n_{0} \times n_{0}} \end{array}] {\tilde{x}}_{k}^{-} \\ = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} \\ 0_{n_{0}} \end{array}] - e_{q}^{T} [\begin{matrix} K_{k, p} \\ K_{k, 0} \end{matrix}] [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}] {\tilde{x}}_{k}^{-} \\ = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} \\ 0_{n_{0}} \end{array}] - e_{q}^{T} [\begin{matrix} K_{k, p} H_{k, p} & 0_{n_{p} \times n_{0}} \\ K_{k, 0} H_{k, p} & 0_{n_{0} \times n_{0}} \end{matrix}] {\tilde{x}}_{k}^{-} = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} - K_{k, p} H_{k, p} {\tilde{x}}_{k, p}^{-} \\ - K_{k, 0} H_{k, p} {\tilde{x}}_{k, p}^{-} \end{array}] \end{matrix}$ C18

Similarly, ${(Δ_{k, q})}_{{\tilde{x}}_{k}^{-}}$ is expanded to the same expression as follows:

$\begin{matrix} e_{q}^{T} (I_{n} - [\begin{array}{c} K_{k, p} \\ K_{k, 0} \end{array}] [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}]) {\tilde{x}}_{k}^{-} = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} \\ 0_{n_{0}} \end{array}] - e_{q}^{T} [\begin{array}{c} K_{k, p} \\ K_{k, 0} \end{array}] [\begin{matrix} H_{k, p} & 0_{m \times n_{0}} \end{matrix}] {\tilde{x}}_{k}^{-} \\ = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} \\ 0_{n_{0}} \end{array}] - e_{q}^{T} [\begin{matrix} K_{k, p} H_{k, p} & 0_{n_{p} \times n_{0}} \\ K_{k, 0} H_{k, p} & 0_{n_{0} \times n_{0}} \end{matrix}] {\tilde{x}}_{k}^{-} = e_{q}^{T} [\begin{array}{c} {\tilde{x}}_{k, p}^{-} - K_{k, p} H_{k, p} {\tilde{x}}_{k, p}^{-} \\ - K_{k, 0} H_{k, p} {\tilde{x}}_{k, p}^{-} \end{array}] \end{matrix}$ C19

Equations (C18) and (C19) confirm that the innovation-based position-domain statistic is equivalent to the snapshot-solution-based SS statistic.

ACKNOWLEDGMENTS

This work was supported by National Research Foundation of Korea grants funded by the Korea government (Ministry of Science and ICT) (No. RS-2024-00354326 and No. RS-2025-02213804).

The authors would like to acknowledge Dr. Juan Blanch for valuable discussions and insightful comments on this work.

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

References

↵
Afia, A. B., Escher, A. C., Macabiau, C., & Roche, S. (2015, April). A GNSS/IMU/WSS/VSLAM hybridization using an extended Kalman filter. In Proceedings of the ION 2015 Pacific PNT Meeting (pp. 719–732). https://www.ion.org/publications/abstract.cfm?articleID=12761
Google Scholar
↵
Blanch, J., Walter, T., & Enge, P. (2017). Theoretical results on the optimal detection statistics for autonomous integrity monitoring. NAVIGATION, 64, 123–137. https://doi.org/10.1002/navi.175
Google Scholar OpenURL
↵
Blanch, J., Walter, T., Enge, P., Lee, Y., Pervan, B., & Rippl, M. (2015). Baseline advanced RAIM user algorithm and possible improvements. IEEE Transactions on Aerospace and Electronic Systems, 51(1), 713–732. https://doi.org/10.1109/TAES.2014.130739
Google Scholar OpenURL
↵
Department of Defense. (2020). Global Positioning System standard position service performance standard (5th ed.). https://www.gps.gov/technical/ps/2020-SPS-performance-standard.pdf
Google Scholar
↵
EMCORE. (n.d.). SDI500 tactical grade inertial measurement unit (IMU) datasheet. http://www.imajteknik.net/uploads/sdi500-sdi505.pdf
Google Scholar
↵
Golub, G. H., & Van Loan, C. F. (1996). Matrix computations (3rd ed.). Baltimore, MD, Johns Hopkins.
Google Scholar
↵
Groves, P. D. (2007). Principles of GNSS, inertial, and multisensor integrated navigation systems. Artech. IEEE. http://ieeexplore.ieee.org/document/9106151
Google Scholar
↵
Gunning, K., Blanch, J., Walter, T., de Groot, L., & Norman, L. (2018, September). Design and evaluation of integrity algorithms for PPP in kinematic applications. In Proceedings of the 31st International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2018) (pp. 1910–1939). https://doi.org/10.33012/2018.15972
Google Scholar
↵
Guttman, L. (1946). Enlargement methods for computing the inverse matrix. Annals of Mathematical Statistics, 17 (3), 336–343. https://doi.org/10.1214/aoms/1177730946
Google Scholar OpenURL
↵
Hajiyev, C. M., & Caliskan, F. (1998, September). Fault detection in flight control systems via innovation sequence of Kalman filter. In UKACC International Conference on Control ‘98 (Conf. Publ. No. 455, pp. 1528–1533). https://doi.org/10.1049/cp:19980456
Google Scholar
↵
International Civil Aviation Organization. (1996). International standards and recommended practices. Annex 10 to the Convention on International Civil Aviation (Vol. 1).
Google Scholar
↵
Joerger, M., Chan, F. C., & Pervan B. (2014). Solution separation versus residual-based RAIM. NAVIGATION, 61(4), 273–291. https://doi.org/10.1002/navi.71
Google Scholar OpenURL
↵
Joerger, M., & Pervan, B. (2013). Kalman filter-based integrity monitoring against sensor faults. Journal of Guidance, Control, and Dynamics, 36, 349–361. https://doi.org/10.2514/1.59480
Google Scholar OpenURL
↵
Lee, D., Nedelkov, F., Akos, D., & Park, B. Barometer based GNSS spoofing detection. (2020, September). In Proceedings of the 33rd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2020) (pp. 3268–3282). https://doi.org/10.33012/2020.17562
Google Scholar
↵
Lee, J., Kim, M., Min, D., Pullen, S., & Lee, J. (2023). Navigation safety assurance of a KF-based GNSS/IMU system: Protection levels against IMU failure. NAVIGATION, 70(4). https://doi.org/10.33012/navi.612
Google Scholar
↵
Liao, J., Li, X., Wang, X., Li, S., & Wang, H. (2021). Enhancing navigation performance through visual-inertial odometry in GNSS-degraded environment. GPS Solutions, 25, 50. https://doi.org/10.1007/s10291-020-01056-0
Google Scholar OpenURL
↵
Miller, B. L., Phillips, C. A., Evans, A. G., & Bibel, J. E. (1993, September). A Kalman filter implementation for a dual-antenna GPS receiver and an inertial navigation system. In Proceedings of the 6th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 1993) (pp. 593–602). https://www.ion.org/publications/pdf.cfm?articleID=4402
Google Scholar
↵
Misra, P., & Enge, P. (2006). Global Positioning System: Signals, measurements, and performance (2nd ed.). Ganga-Jamuna Press, Lincoln.
Google Scholar
↵
Odendaal, H. M., & Jones, T. (2014). Actuator fault detection and isolation: An optimized parity space approach. Control Engineering Practice, 26, 222–232. https://doi.org/10.1016/j.conengprac.2014.01.013
Google Scholar OpenURL CrossRef
↵
Patton, R. J., & Chen, J. (1991). A review of parity space approaches to fault diagnosis. IFAC Proceedings Volumes, 24, 65–81. https://doi.org/10.1016/S1474-6670(17)51124-6
Google Scholar OpenURL
↵
Psiaki, M. L. (2021). Navigation using carrier Doppler shift from a LEO constellation: TRANSIT on steroids. NAVIGATION, 68(3), 621–641. https://doi.org/10.1002/navi.438
Google Scholar OpenURL
↵
RTCA. (2001). Minimum operational performance standards for Global Positioning System/wide area augmentation system airborne equipment. RTCA/DO-229C.
Google Scholar
↵
Tanil, C., Khanafseh, S., Joerger, M., & Pervan, B. (2017). An INS monitor to detect GNSS spoofers capable of tracking vehicle position. IEEE Transactions on Aerospace and Electronic Systems, 54(1), 131–143. https://doi.org/10.1109/TAES.2017.2739924
Google Scholar OpenURL
↵
Travis, W., Simmons, A. T., & Bevly, D. M. (2005, June). Corridor navigation with a LiDAR/INS Kalman filter solution. In IEEE Proceedings. Intelligent Vehicles Symposium (pp. 343–348). https://doi.org/10.1109/IVS.2005.1505126
Google Scholar
↵
Van Trees, H. L. (1968). Detection, estimation, and modulation theory (Vol. 1). Wiley.
Google Scholar
↵
Wang, S., Zhan, X., Zhai Y., Zheng, L., & Liu, B. (2022). Enhancing navigation integrity for urban air mobility with redundant inertial sensors. Aerospace Science and Technology, 126, 107631. https://doi.org/10.1016/j.ast.2022.107631
Google Scholar OpenURL
↵
Whitty, C., & Walport, M. (2018). Satellite-derived time and position: A study of critical dependencies. Government Office for Science: London, UK.
Google Scholar
↵
Xu, Q., Lv, J., & Gao, Z. (2022, September). Evaluation on low-cost GNSS/IMU/vision integration system in GNSS-denied environments. In 2022 IEEE 12th International Conference on Indoor Positioning and Indoor Navigation (IPIN) (pp. 1–6). https://doi.org/10.1109/IPIN54987.2022.9918147
Google Scholar

PDF

Previous Next

Abstract

1 INTRODUCTION

2 BACKGROUND

2.1 KF Estimate Error Under Fault-Free Conditions

2.2 KF Estimate Error Under Sensor Fault Conditions

2.3 KF Innovation-Based IM Algorithm

3 INNOVATION-BASED OPTIMAL DETECTOR AGAINST DYNAMIC SENSOR FAULTS

3.1 Optimization Problem for Detection Regions

3.2 Reparameterization

3.3 Optimal Detection Region for One Coordinate

3.4 Equivalence with the SS Statistic

4 PERFORMANCE EVALUATION

4.1 Integrity Risk Under IMU Faults

4.1.1 Derivation of the Worst-Case Innovation Bias

4.1.2 Integrity Risk Computation

4.2 Performance Analysis Results

4.3 Sensitivity Analysis of False Alert Rate Requirement Allocation

5 PERFORMANCE AND COMPUTATIONAL LOAD COMPARISON WITH THE FULL SS MONITOR

5.1 PL Computation

5.2 Comparison with the Full SS Across Sensor Configurations

5.2.1 Simulation Conditions

5.2.2 Simulation Results

5.2.3 Discussion

6 CONCLUSION

HOW TO CITE THIS ARTICLE

APPENDIX | A

APPENDIX | B

Neyman–Pearson Lemma

Optimal Detection Region for a Bias δp+

Symmetry of Ωδp+

Further Simplification of Ωδp+

APPENDIX | C

Equivalence of Current-Time Measurement Noise Components

Equivalence of a Priori Estimate Error Components

ACKNOWLEDGMENTS

References

In this Issue

NAVIGATION: Journal of the Institute of Navigation

Alerts for this Article

Email This Article

Citation Tools

Kalman Filter Innovation-Based Optimal Integrity Monitoring Against Unknown Bias in Time Updates

Citation Manager Formats

Share This Article

Kalman Filter Innovation-Based Optimal Integrity Monitoring Against Unknown Bias in Time Updates

SAVE ITEM IN SELECTED FOLDER.

Jump to section

Related Articles

Cited By...

More in this TOC Section

Similar Articles

Keywords

Optimal Detection Region for a Bias $δ_{p}^{+}$

Symmetry of $Ω_{δ_{p}^{+}}$

Further Simplification of $Ω_{δ_{p}^{+}}$