Abstract
Prior work established a model for uncertain Gauss-Markov (GM) noise that is guaranteed to produce a Kalman filter (KF) covariance matrix that overbounds the estimate error distribution. The derivation was conducted for the continuous-time KF when the GM time constants are only known to reside within specified intervals. This paper first provides a more accessible derivation of the continuous-time result and determines the minimum initial variance of the model. This leads to a new, non-stationary model for uncertain GM noise that we prove yields an overbounding estimate error covariance matrix for both sampled-data and discrete-time systems. The new model is evaluated using covariance analysis for a one-dimensional estimation problem and for an example application in Advanced Receiver Autonomous Integrity Monitoring (ARAIM).
SUMMARY OF MAIN RESULT
Let πΊ be the Kalman filter (KF) estimate error covariance matrix and π be the true estimate error covariance matrix. The inequality πΊ β₯ π means that the predicted variance πΆππΊπΆ is greater than or equal to the true variance πΆπππΆ for any real vector πΆ.
Suppose that the measurement and process noise components are known to be first-order Gauss-Markov random processes (GMPs). These processes are completely specified by a time constant π, steady-state variance π2 and initial variance 
, and propagate according to the difference equation
where π is an arbitrary time index, Ξπ‘ = π‘π+ 1 βπ‘π is the discrete-time sampling interval, WGN(0, 1) indicates zero-mean white Gaussian noise with unit variance and π(0, 
) denotes a zero-mean normal random variable with variance 
.
When π and π2 are only known to reside in the intervals [πmin,πmax] and 
, a state-augmented KF that models ππ with
is guaranteed to produce a covariance matrix πΊ β₯ π.
1 INTRODUCTION
Safety-critical GNSS applications must ensure that the probability of the estimate error exceeding predefined bounds is acceptably small. This probability is often referred to as integrity risk. In aviation navigation applications, cumulative distribution function (CDF) overbounding has been used extensively over the past 10-20 years to derive upper bounds on integrity risk for snapshot least-squares, GPS-based navigation (DeCleene, 2000; Rife et al., 2004, Rife et al., 2006; Blanch et al., 2019). However, state estimation using a recursive filter can be more accurate and enables integration with additional sensors like an inertial measurement unit (IMU). One challenge with sequential estimation is how to conservatively account for time-correlated noise like GPS multipath. Precise multipath models are very difficult to obtain, as illustrated in Pervan et al. (2017), for an ensemble of empirically derived multipath autocorrelation functions (ACFs). Thus, even using the best models available, there is uncertainty when dealing with time-correlated noise.
Existing CDF overbounding methods do not account for model uncertainty over time, which has prompted the development of new techniques. Theoretical approaches were developed in (Rife & Gebre-Egziabher, 2007; Pulford, 2008) to bound the integrity risk of linear systems driven by a spherically symmetric random process. The authors in Langel et al. (2014) derive an integrity risk bound for the KF when measurement and process noise are Gaussian with ACFs that are unknown but can be upper- and lower-bounded. The methods in (Rife & Gebre-Egziabher, 2007; Pulford, 2008; Langel et al., 2014) are attractive because they do not require any knowledge about the mathematical form of the input noise ACFs. However, the price to pay is that (Rife & Gebre-Egziabher, 2007; Pulford, 2008; Langel et al., 2014) all provide a batch solution in the sense that they require storing matrices whose dimensions grow without bound at a rate dictated by sensor sampling rates. For real-time applications, these techniques are only suitable for short duration applications. If high-rate sensors like an IMU are being used, the maximum allowable duration could be on the order of minutes before computational and memory constraints begin to prohibit real-time operation. One way to derive a practical algorithm is to stipulate that the measurement and process noise ACFs are known up to a finite set of uncertain parameters.
State estimation with uncertain parameters has been studied extensively in the robust estimation literature. Guaranteed cost filtering emerged in the 1990s as a technique for defining a robust filter whose estimate error variance is guaranteed to be smaller than a given bound. These filters often require determination of one or more scale parameters using an optimization algorithm and can display unstable behavior. Dependence on scale parameters and questionable stability characteristics means that significant effort is required up front to design the estimator. Real-time capability is also a challenge, given the need to numerically solve an optimization problem. For these reasons, robust estimators have not been widely adopted so far by the navigation community to upper bound KF integrity risk when the input noise ACFs are uncertain.
This paper builds on the work in Tupysev et al. (2009), where a guaranteed upper bound on the KF estimate error variance is derived when process and measurement noise components are first-order GMPs with unknown, but bounded, time constants. Reference Tupysev et al. (2009) shows that in order to guarantee that the KF covariance matrix overbounds the estimate error, both the minimum and maximum time constants must be used in the KF matrices. The derivation is conducted for the continuous-time KF, and a covariance analysis is performed for an IMU-aided navigation problem to show that the estimate error variance is bounded. The authors in Tupysev et al. (2009) indicate that the initial variance of their model must be inflated relative to the true GM variance but do not specify how much inflation is necessary. They also did not discuss how their results apply for the more common case of discrete-time systems. We will address both deficiencies.
The paper begins with the problem statement and motivation for its study in Sections 2 and 3. Starting from a sampled-data system model, Sections 4 through 5.1 provide a more accessible derivation of Equation (23) in Tupysev et al. (2009), following many of the authorsβ key developments. In Section 5.2, we derive the minimum initial variance for the GM model (Equation (25)) that yields an overbounding covariance matrix for the mixed continuous-discrete KF. This expression is a previously unknown result that leads to a new, non-stationary model for GM noise that we prove in Section 5.3 results in an overbounding estimate error covariance matrix for the discrete-time KF. The proof is a new development that rigorously justifies using our model to account for uncertain GM noise in a discrete-time setting. An illustrative example shows that the non-stationary model provides a significantly tighter bound on the estimate error variance compared to a stationary model.
One application where time-sequential estimation is being explored for high integrity navigation is Vertical Advanced Receiver Autonomous Integrity Monitoring, or V-ARAIM, which seeks to provide worldwide vertical guidance of aircraft using dual-frequency measurements from GPS and Galileo. It was shown recently that V-ARAIM performance can be significantly improved by processing time-sequences of measurements (Joerger & Pervan, 2015, 2016). In Section 6, we carry out a performance analysis for time-sequential V-ARAIM. Multipath and residual tropospheric errors are first-order GMPs with time constants that reside within their own specified intervals. We use our new GM model to establish global availability maps that show sequential V-ARAIM outperforms traditional snapshot ARAIM even though the multipath and residual tropospheric error time constants are uncertain. Concluding remarks and suggestions for future study are given in Section 7.
2 PROBLEM SETUP
Consider the state estimation problem for the sampled-data system
1
where π β βπ is the state vector, ππ β βπ is the measurement vector and π(π‘) β βπ, ππ β βπ are the process and measurement noise vectors, respectively. The matrices π, π, ππ and ππ are known and of appropriate dimension.
A sampled-data system model is used for most of the theoretical development in this paper because many systems encountered in practice are described by Equation (1). State dynamics are often derived from physical laws (e.g., Newtonβs laws of motion) and take the form of continuous-time differential equations, whereas external measurements of a system are usually taken with digital sensors. Recognizing that Equation (1) is typically converted entirely to discrete-time prior to developing an estimator, results for sampled-data systems are extended to discrete-time systems in Section 5.3.
The vectors π and ππ are known to be linear combinations of a first-order GMP vector π β βπ and zero-mean WGN vectors π(π‘) β βπ and ππ β βπ . That is
2
such that πw(π‘) β βπΓπ and πv,π β βπ Γπ are known matrices, π(π‘) β βπΓπ is the power spectral density matrix of π(π‘), ππ β βπ Γπ is the covariance matrix of ππ and πΏ(π‘), πΏkl are the Dirac delta function and Kronecker delta, respectively. π(π‘) and ππ are unknown, but we assume that 
and 
can be determined.
For the GMP vector π, π and π are (Gelb et al., 1974)
3
where 
and ππ,π = 1,β¦,π are only known to reside in specified intervals. The goal of the paper is to develop a KF whose predicted estimate error covariance matrix accounts for uncertainty in 
and ππ. In subsequent sections, β(π‘)β will be omitted to lighten the notation.
State augmentation (Anderson & Moore, 1974) is commonly used to account for correlated noise and involves appending π to π to form the new linear system
4
With ππ = [ππ ππ], Equation (4) can also be written as
5
Equation (5) is a linear system driven by zero-mean WGN. Thus it is permissible to estimate π using a KF. Since 
and ππ are unknown, it is not clear what values should be used to define the KF. One criterion is to use parameter values so the predicted estimate error distribution overbounds the true error distribution over the range of admissible values.
To define an overbounding distribution, let 
be an estimate of a specified state π¦ and define 
as the associated estimate error. Then for a maximum tolerable error ππ¦, an overbounding distribution is one that produces an upper bound on the integrity risk π(|ππ¦| β₯ ππ¦). Given that Equation (5) is a linear system driven by zero-mean WGN and the fact that the KF is a linear unbiased estimator, it follows that 
. In this case, an upper bound on integrity risk is equivalent to an upper bound on the estimate error variance 
. Therefore, subsequent sections of the paper will focus on developing an upper bound on the estimate error variance for any state of interest.
3 MOTIVATIONAL EXAMPLE
Before developing a rigorous approach for defining an overbounding distribution in the presence of uncertainty, it is instructive to address the fallacies in current thinking. The conventional wisdom is that a KF using upper bound parameter values will produce an estimate error covariance matrix that overbounds the true error distribution. We use a simple one-dimensional estimation problem to demonstrate that this heuristic is not always accurate.
Consider the vehicle in Figure 1 moving at constant speed π’ along the π₯-axis from an initial point π0. Noisy measurements of the current position are available, described by the model
6
One-dimensional position and velocity estimation
where π£π is the sum of a first-order GMP and zero-mean WGN. That is, π£π = ππ +ππ such that 
and
7
The variances are known and given by 
and 
. However, the time constant π is only known to lie in the interval [10,100] s. Noting that π0 and π’ are constants, state augmentation results in the estimation problem
8
such that the initial estimate error covariance matrix π0 is diagonal with nonzero elements: [10 m2,1m2βs2,1m2].
The KF is constructed assuming that π = 100 seconds. However, suppose that π is actually equal to 50 seconds. We want to compare the KFβs predicted estimate error variance to the true variance in this case. The true estimate error variance is obtained through Monte Carlo simulation. Measurements of current position are generated according to the true measurement model and passed into the KF. The variance is then estimated from the sample estimate error distribution. We saw that the change in estimated variance between 800,000 trials to one million trials was indistinguishable. Thus, one million trials is sufficient to accurately determine the true variance.
Figure 2 shows the estimate error standard deviation predicted by the KF minus the true standard deviation for the initial position and speed states. One might expect to see that the estimated-minus-true difference is positive at all times when assuming that π = πmax. This is clearly not the case. For the initial position state, the assumption that π = πmax is eventually conservative after 60 seconds have elapsed. For the speed state, assuming that π = πmax always leads to optimistic predictions of the estimate error standard deviation. This simple example demonstrates that a more rigorous approach is needed to establish an overbounding estimate error distribution when correlated measurement error models are uncertain.
Difference between predicted and true estimate error standard deviations for a KF that assumes π = πmax
4 TRUE COVARIANCE MATRIX
This section derives propagation equations for the true estimate error covariance matrix when the KF uses models that do not correspond to the actual process and measurement noise statistics. In what follows, the notation (β )π|π will be used to denote a quantity at time index π based on all measurements up to and including time index π. The KF for the system in Equation (5) is defined by the recursion
9
The filter designer must specify 
(which is contained in 
) and 
. These matrices generally will not describe the true measurement and process noise statistics because the parameters defining them are uncertain. The KF covariance matrix 
therefore does not accurately capture the true estimate error distribution, and a new covariance matrix π must be defined that accounts for model uncertainty.
When 
, Appendix A shows that the KF estimate error vector 
propagates according to 
, where 
. Thus, πΊ and π must be considered simultaneously when determining the true covariance matrix for πΊ. The correlation between πΊ and π is properly accounted for by forming the augmented vector ππ = [πΊπ ππ]. This method of accounting for model uncertainty is similar to the approach in Gelb et al. (1974). The only difference is that we focus on uncertainty in π, whereas Gelb et al. (1974) considers uncertainty in the entire transition matrix π
. Appendix A shows that π = πΈ[πππ] evolves according to
10
such that
11
Equations (10) and (11) contain both true and estimator matrices and thus reflect that π, π, π and π
are uncertain. For example, 
and 
go into determining ππ, but π and ππ also appear explicitly. However, π cannot be computed in practice because of the unknown matrices involved. We therefore seek to define a covariance matrix πΊ free of dependence on unknown matrices such that (πΊ β π) β₯ 0, that is, πΊβπ is positive semi-definite. This criterion ensures that the predicted error variance for any state is an upper bound on the true variance.
5 MODELING UNCERTAIN GM NOISE
The goal of this section is to define a propagation structure for πΊ and use the criteria (πΊβπ) β₯ 0 to derive a model for uncertain GM noise. Following the approach in Tupysev et al. (2009), we stipulate that πΊ propagates in a manner similar to Equations (10) and (11), that is,
12a
12b
This structure was chosen because, as we will see, it enables 
, 
and the initial covariance matrix πΊ(0) to be defined explicitly from the requirement (πΊ β π) β₯ 0. The covariance error matrix π« = πΊβπ propagates as
13a
13b
Equation (13a) is a matrix differential equation with initial condition π«πβ1|πβ1 and whose solution at π‘ = π‘π is π«π|πβ1. The matrix π«π|πβ1 is positive semi-definite if π«πβ1|πβ1 β₯ 0 and 
(see Dieci & Eirola, 1994). The quadratic structure in Equation (13b) ensures that π«π|π β₯ 0 when π«π|πβ1 β₯ 0 and 
Therefore, we must have an initial matrix π«(0) β₯ 0 and 
to guarantee that π«(π‘) β₯ 0 for all π‘.
5.1 Continuous-time model
We will now provide a more accessible derivation of the continuous-time model for uncertain GM noise given in Tupysev et al. (2009), following many of the authorsβ key developments. In preparation for the first step, defining 
readers should re-familiarize themselves with the definitions in Equation (3).
Let 
formed with the upper bound values for 
and the lower bound values for ππ. Letβs also define 
as the matrix π populated with the true time constants ππ and as yet undetermined variances 
. Then with 
defined as
14
it is clear that 
for all admissible values of 
and ππ. However, 
could never be constructed because it depends on the unknown matrix 
. Therefore, we must use a mathematical approach to make πΊ insensitive to 
. To see how this can be accomplished, form the partitioned matrix 
so that Equation (12a) can be written in expanded form as
15a
15b
15c
where
16
We are only interested in propagating πΊπ₯, the block of πΊ corresponding to πΊ. In its current form, Equation (15a) cannot be propagated because it depends on the unknown matrix Ξπ
. We can eliminate Ξπ
by forcing πΊπ₯π(π‘) = π for all π‘. To achieve this, let 
and impose the constraint
17
which simplifies Equation (15b) to the unforced linear differential equation 
. Next, set the initial condition πΊπ₯π(0) = π. The only solution to an unforced linear differential equation with zero initial condition is the trivial solution. Therefore, it must be that πΊπ₯π(π‘) = 0 for all π‘, as desired.
Equation (17) has additional ramifications when we consider the alternative form πΊπ = β(Ξπ)β1π. The matrices Ξπ and π are time-invariant, which implies that 
. Making the substitutions 
into Equation (15c) results in the additional constraint
18
Noting that π,
, π and 
are diagonal matrices, Equation (18) is equivalent to the set of scalar conditions
19
It is worth reiterating that ππ and 
are unknown, 
is a quantity we need to determine, and 
must be greater than or equal to 
. To determine 
, first solve Equation (19) for 
, resulting in
20
We know that variances are never less than zero. The denominator in Equation (20) indicates that the only way to guarantee that 
for any ππ β [ππ,min,ππ,max] is to set 
equal to ππ,max.
With 
, the inequality in Equation (20) simplifies to 1 β₯ ππβππ,max, which is satisfied for all admissible values of 
. Thus, we satisfy the requirement that 
. Having determined the diagonal elements of 
and 
, we can conclude that πΊπ₯ β₯ ππ₯ when first-order GMPs are modeled according to
21
Equation (21) is identical to Equation (23) in Tupysev et al. (2009). This completes the proof of the main result in Tupysev et al. (2009). At this point, the model in Equation (21) is incomplete without specifying the initial variance. The reader also may be wondering about the role of 
. It is comforting that 
is not required to propagate πΊπ₯ because its dependence on the unknown time constant ππ means that 
can never be constructed. However, a previously unknown fact is that 
places constraints on the initial variance for the GM model. This new result is covered next.
5.2 Initialization
Reference Tupysev et al. (2009) indicates that the initial variance of the GM model must be inflated relative to the true variance but does not provide any insight as to how much inflation is required. In this section, we derive the tightest bound on the initial variance for the model in Equation (21) that still guarantees πΊπ₯ β₯ ππ₯.
Section 5.1 showed that the initial cross-covariance πΊπ₯π(0) must be zero to eliminate any dependence of πΊπ₯ on the unknown matrix Ξπ . Assuming that πΊπ,0 and πΊπ,0 are initially uncorrelated with covariance matrices πΊπ,0 and πΊπ,0, the requirement πΊπ₯π(0) = π leads to the initial covariance matrix
22
where πΊπ,0 is assumed to be known (e.g., from an initialization algorithm), πΊπ,0 is what we seek to define in this section, and 
is a diagonal matrix populated with the 
in Equation (20).
Appendix B shows that the true initial covariance matrix is
23
such that ππ,0 is a diagonal matrix comprised of the unknown variances 
. We showed earlier that πΊ(0) β₯ π(0) is a necessary condition to ensure that πΊ(π‘) β₯ π(π‘) for all π‘. We assume that πΊπ,0 can be specified so that πΊπ,0 β₯ ππ,0. Then using the fact that a block diagonal matrix is positive semi-definite if and only if each diagonal block is positive semi-definite, πΊπ,0 must be specified so that
24
Denoting an arbitrary diagonal element of πΊπ,0 as 
, Appendix B shows that Equation (24) is satisfied if
25
It is worth noting that the model in Equation (21) is stationary when initialized with a variance 
that is greater than the lower bound in Equation (25). A stationary model for uncertain GMPs is appealing because GMPs are in fact stationary. However, by initializing with the minimum variance in Equation (25), it is possible to obtain a tighter bound on the estimate error variance using a non-stationary model. This will be demonstrated in Subsection 5.4 for the one-dimensional example considered in Section 3.
5.3 Discrete-time model
Even though Equation (1) may be the native description of a system, it is often converted entirely to discrete-time when designing and implementing a KF. We will prove in this section that the discretized form of Equation (21) yields an overbounding estimate error covariance matrix for the discrete-time KF.
Consider a discrete-time system analogous to Equation (4)
26
where 
and 
. As before, we assume that 
and πΊπ,0 β₯ ππ,0 can be specified. Suppose that 
is estimated with a KF that models ππ using the discretized form of Equation (21). From the continuous-to-discrete mapping (Rogers, 2003)
27
with π = πβΞπ‘βπ, the discretized form of Equation (21) is seen to be
28
Following a similar approach to Subsection 5.1, Appendix C shows that the covariance error matrix π« = πΊβπ for the discrete-time KF propagates as
29a
29b
Given that our analysis of sampled-data systems also considered discrete-time measurements, it is no surprise that Equation (29b) is identical to Equation (13b). However, Equation (29a) differs from the continuous version in Equation (13a) in that there is an additional matrix dependent on Ξπ π that must now be considered.
To show that π«π|π β₯ 0, first note that whether the original system is a sampled-data system or a discrete-time system, initialization is the same. Thus, the results from Subsection 5.2 also apply here, that is, initializing with πΊπ,0 β₯ ππ,0 and the variances in Equation (28) ensures that πΊ0|0 β₯ π0|0. Second, notice that the quadratic structure in Equation (29b) guarantees that π«π|π β₯ 0 if π«π|πβ1 β₯ 0. Therefore, if we can establish that π«π+ 1|π β₯ 0, this will prove that the discrete-time KF covariance matrix overbounds the estimate error distribution.
Determining whether π«π+ 1|π β₯ 0 comes down to showing that
30
which we prove is true in Appendix D. Hence, we can conclude that for a discrete-time system, the model given in Equation (28) yields a KF covariance matrix that overbounds the estimate error distribution when the GM time constants are uncertain.
5.4 Motivational example revisited
A visual representation of the paperβs main result is obtained by revisiting the motivational example from Section 3. The KF is run twice, once with the non-stationary model in Equation (28) and once with the stationary model where 
. Just as before, the true estimate error variance is obtained through Monte Carlo simulation.
Figure 3 shows the estimate error standard deviation predicted by the KF and the true standard deviation using each model. The predicted standard deviation is always greater than the true standard deviation, indicating that both models provide an upper bound on the estimate error variance. For the speed state, there is almost no visible benefit of using the non-stationary model. However, the benefit is clear for the initial position state. Not only is the variance bound smaller using the non-stationary model compared to the stationary model, but it is also less conservative. That is, the gap between the predicted and true standard deviations is smaller when using the non-stationary model for GM noise.
KF standard deviation compared to truth when using stationary and non-stationary bounding GMP models
Remember that πΊ was derived so that πΊ β₯ π. This means that the ellipsoid corresponding to πΊ circumscribes the ellipsoid associated with π for any admissible π. Figure 4 shows the true covariance matrix obtained from Monte Carlo simulation in gray and the KF covariance matrix in black, at an elapsed time of 10 seconds, when using the stationary and non-stationary bounding models. The true covariance matrices are different for each bounding model because they are obtained using two separate estimators, derived from two different values for the initial GM state variance. We can see that the KF covariance ellipse does indeed circumscribe the true ellipse. The white space between π and πΊ along the π₯-axis is smaller for the non-stationary bound, and about the same for both bounds along the π¦-axis. This reiterates the conclusion from Figure 3 that the non-stationary bound is tighter than the stationary bound for the initial position state, but provides little benefit for the speed state.
True and bounding covariance ellipses for the motivational example in Section 3. Ellipses are shown at an elapsed time of 10 seconds
The ability to overbound the estimate error distribution for any linear combination of state vector components is critical for certain applications. In aircraft precision navigation, for example, horizontal integrity risk assessment requires overbounding a combination of the north and east components of positioning error. For the first time, the methods developed in this paper provide a practical means to determine such an overbound in the presence of uncertain GM noise. The next section will implement the non-stationary model in an example application of aircraft navigation using ARAIM (Working Group C - ARAIM Technical Subgroup, 2012, 2015, & 2016).
6 APPLICATION TO SEQUENTIAL ARAIM PERFORMANCE ANALYSIS
The baseline ARAIM algorithm is a βsnapshotβ implementation that uses dual-frequency multi-constellation measurements at one time-instant to achieve LPV-200 requirements, that is, requirements for localizer performance with vertical guidance down to 200 feet above the runway (Working Group C - ARAIM Technical Subgroup, 2015 & 2016).
However, when nominal GPS and Galileo constellations are depleted, LPV-200 can only be sparsely achieved using snapshot ARAIM (Working Group C - ARAIM Technical Subgroup, 2015). βSequential ARAIMβ addresses this limitation by processing measurements over time (Joerger & Pervan, 2015, 2016, & 2020). One major challenge in sequential ARAIM is to derive robust models of measurement errors over time. In prior work (Joerger & Pervan, 2016 & 2020), a βbias-plus-rampβ model was derived for satellite clock and orbit error dynamics using nine months of data. In parallel, assumptions were made on the time correlation of tropospheric and multipath errors. Because no method was available to account for uncertainty in correlation time constants, the largest measured value was employed for multipath, and sensitivity to residual tropospheric delay was analyzed (Joerger & Pervan, 2015, 2016, & 2020). But multipath time constants can take a wide range of values as reported in Figures 21-22 of Pervan et al. (2017).
In this paper, we implement our new non-stationary GM model to account for uncertainty in the correlation time constant of multipath and tropospheric errors in a sequential ARAIM performance analysis. This section briefly describes the main assumptions of the multipath and troposphere error models, and then focuses on quantifying the impact of error correlation model uncertainty on ARAIM performance. Readers may refer to Joerger and Pervan (2015, 2016, & 2020) for details on the sequential ARAIM implementation but should be aware that these references focus on the batch implementation. The algorithms in Joerger and Pervan (2015, 2016, & 2020) can equivalently be implemented using a KF because sequential ARAIM risk evaluation only requires current-time state estimates, which coincide for the KF and batch estimators. In KF-based sequential ARAIM, a bank of KFs is used to evaluate current-time full-set and subset solutions for multiple hypothesis solution separation (MHSS). Readers interested in the details of the KF versus batch implementations may refer to Joerger and Pervan (2013); Joerger (2009); and Crassidis and Junkins (2004).
6.1 Sequential ARAIM assumptions
Let π be the number of visible and healthy GPS and Galileo satellites at time π. The sequential ARAIM state propagation and measurement equations can be written in the form of Equations (4) and (5), where
π βΆ is the πΓ1 non-augmented state vector composed of three-dimensional antenna position coordinates and GPS/Galileo receiver clock biases, time-invariant, float-valued carrier phase cycle ambiguities, and time-invariant clock and orbit error bias and ramp parameters for each satellite. In this case, π = 5+3π .
π βΆ is the πΓ1 augmentation vector, with π = 3π , made of code and carrier multipath and troposphere error states for each satellite.
ππ βΆ is the πΓ1 measurement vector, with π = 2π , comprised of unfiltered ionosphere-free code and carrier ranging measurements for all satellites at time π.
For the measurement equation, ππ = ππ, and ππ and πv,π are
31a
31b
where ππ is an π Γ5 geometry matrix whose first three columns are made of unit line-of-sight row vectors, and whose next two columns have zeros and ones identifying whether the receiver clock bias is for a GPS or Galileo satellite; π‘0 and π‘π are the filter initialization time and current time, respectively; πM,π is a diagonal matrix of elevation-dependent multipath coefficients described in Equations (8) and (9) of Joerger and Pervan (2016); πT,π is a diagonal matrix of elevation-dependent troposphere error coefficients given in Equations (6) and (7) of Joerger and Pervan (2016). The measurement noise covariance matrix 
is diagonal, with elevation-dependent elements given in Equations (10) and (11) of Joerger and Pervan (2016).
The non-augmented process equation captures the fact that all elements of π are constant, except for the five position and receiver clock bias states, for which the time propagation is unknown. Thus, π = ππΓπ,π = ππ,πw = ππΓπ, and the process noise power spectral density matrix π has large values on the first five diagonal elements, and zeros otherwise.
For the GM state vector, π and π are
32
33
where 
, 
and 
are defined in Table 1. The values of ππ,ππ and ππ are unknown but bounded following the inequalities
34
Multipath and troposphere error models
Values for the bounds on ππ,ππ and ππ are provided in Table 1.
In order to evaluate sequential ARAIM performance, we assume nominal error model parameter values given in Joerger and Pervan (2020). A constant filtering period of 10 min is assumed throughout the section.
Additional ARAIM error model parameters, such as the bounded bias accounting for non-Gaussian ranging errors (Working Group C - ARAIM Technical Subgroup, 2015 & 2016), are included in the simulation as explained in Joerger and Pervan (2020) but are not described here because they are not directly relevant to the problem of modeling uncertain, time-correlated noise.
6.2 Sequential ARAIM integrity monitoring performance analysis
In this subsection, we compare sequential ARAIM performance using the non-stationary GM model to the performance of snapshot ARAIM. Nominal simulation parameters are defined in Joerger and Pervan (2020), and include:
an elevation mask of 5 degrees
an integrity risk requirement of 10β7 per approach
a continuity risk requirement of 4Γ10β6 per approach
a vertical alert limit π of 35 m (unless otherwise stated)
a prior probability of satellite fault of 10β5
a prior probability of constellation fault πconst of 10β4 (unless otherwise stated)
depleted constellations of β24-1β GPS satellites and β24-1β Galileo satellites (Stanford University GPS Lab)
In Figure 5, the integrity risk is evaluated using sequential MHSS (Working Group C - ARAIM Technical Subgroup, 2016; Joerger & Pervan, 2020) over 24 hours, at an example Blacksburg, VA, location (37.2β¦ N, β80.4β¦ E), assuming dual-frequency measurements from GPS and Galileo. In this figure, in order to accentuate differences between implementations, we use a vertical alert limit of 10 m (instead of 35 m) and a prior probability of constellation faults of 10β8.
Integrity risk bound for snapshot versus sequential ARAIM with the non-stationary GM model. Bounds are computed using depleted constellations, πconst = 10β8 and π = 10 m
Figure 5 shows the significant integrity risk reduction obtained using sequential ARAIM (thick black curve) as compared to snapshot ARAIM (thin black curve with diamond markers). The fraction of time where the integrity risk curves are below the horizontal dotted line is the availability. In this case, availability is 60% for sequential ARAIM compared to 26% for snapshot ARAIM. The main conclusion from Figure 5 is not only that sequential ARAIM outperforms snapshot ARAIM as in Joerger and Pervan (2020), but also that we are still able to achieve substantial availability improvement even after accounting for our lack of knowledge in ππ, ππ and ππ.
Figures 6 and 7 display global availability maps for a 10β¦ Γ10β¦ latitude-longitude grid of locations, for satellite geometries simulated at regular 10-minute intervals over a 24-hour period. Availability is computed at each location as the fraction of time where the integrity risk bound is lower than the 10β7 requirement. The larger presence of white areas in Figure 7 compared to Figure 6 clearly indicates that snapshot ARAIM is outperformed by sequential ARAIM, even after accounting for measurement error time correlation uncertainty. The worldwide availability metric given in the figure captions is the weighted coverage of 99.5% availability. Coverage is defined as the percentage of grid point locations exceeding 99.5% availability. The coverage computation is weighted at each location by the cosine of the locationβs latitude, because grid point locations near the equator represent larger areas than near the poles.
Availability map for snapshot ARAIM using depleted constellations, πconst = 10β4,π = 35 m (coverage of 99.5% availability is 84.6%)
Availability map for sequential ARAIM with the non-stationary GM model using depleted constellations, πconst = 10β4,π = 35 m (coverage of 99.5% availability is 97.8%)
Table 2 lists worldwide coverage of 99.5% availability, and of 95% availability (given in parentheses). It shows that the coverage of 99.5% availability increases from 84.6% for snapshot ARAIM to 97.8% for sequential ARAIM with the new non-stationary GM model, a 13% improvement for the case considered here. It is also worth pointing out that in Langel et al. (2019), we showed that if we knew π = πmax, 99.5% availability coverage for sequential ARAIM would be 99.5%. Thus, we sacrifice 1.7% of availability coverage due to time correlation uncertainty. These results would further improve with better knowledge of the measurement error sources, that is, if the difference between πmax and πmin could be reduced.
Coverage of 99.5% availability and coverage of 95% availability (in parentheses)
7 CONCLUSIONS
A new model for GM noise was developed that produces a KF covariance matrix that overbounds the estimate error distribution when the GM variances and time constants are uncertain. We proved that our model is applicable to both sampled-data and discrete-time systems, and demonstrated its performance for a one-dimensional estimation problem. For an aircraft navigation application using sequential ARAIM with our new model, it was shown that under specific circumstances, a larger than 10% improvement in 99.5% availability coverage is achievable compared to conventional snapshot ARAIM.
Future work will apply the results of this paper to other filtering applications like GNSS/INS integrated navigation. Research is also underway extending the ideas presented in the paper beyond first-order Gauss-Markov noise, with the goal of developing practical overbounding methods that do not rely on precise mathematical knowledge of the underlying time correlation.
HOW TO CITE THIS ARTICLE
Langel S, GarcΓa Crespillo O, Joerger M. Overbounding the effect of uncertain Gauss-Markov noise in Kalman filtering. NAVIGATION. 2021;68(2):259β276. https://doi.org/10.1002/navi.419
ACKNOWLEDGMENTS
We would like to thank Dr. Juan Blanch of Stanford University for his comments and suggestions. His review was instrumental in ensuring that the technical development of the paper was conveyed as clearly as possible.
APPENDIX A: TRUE COVARIANCE MATRIX
This appendix derives the true estimate error covariance matrix for a state-augmented KF when there are uncertain parameters in the augmented state dynamic models. With 
, the transition matrix in Equation (4) can be written as
A1
Then the dynamic model 
in Equation (5) becomes 
. Noting that the KF estimate vector propagates as 
, the estimate error vector 
satisfies the differential equation
A2
Defining Ξπ π = [π Ξππ] and incorporating the dynamic model for the Gauss-Markov states results in the following propagation equation
A3
The estimate vector after the measurement update is 
. It is straightforward to show that the estimate error vector is πΊπ|π = (πβ ππππ)πΊπ|πβ1 βππππππ. When combined with ππ = ππ, the update equation is
A4
The covariance matrix π = πΈ[πππ] propagates between measurements according to Simon (2006)
A5
where
A6
During a measurement update, the covariance matrix is updated through
A7
Equations (A5) and (A7) are identical to Equations (10) and (11) stated in Section 4.
APPENDIX B: INITIAL VARIANCE BOUND
This appendix derives the minimum allowable initial variance for the process in Equation (21). The derivation is based on the requirement that π«(0) β₯ 0 and πΊπ₯π(0) = π. If no prior information exists to initialize π (the case considered here), the initial estimate for the augmented states is 
. We also assume that 
and Γ’0 are initially uncorrelated, so that 
.
This leads to the following expected values
B8
where we have used the fact that because Γ’0 = 0, 
, 
and 
are all equal to π. Using these relations, the true initial covariance matrix of 
is
B9
The requirement πΊπ₯π = π implies that
B10
The goal is to determine πΊπ,0 so that πΊ(0) β₯ π(0). Subtracting π(0) from πΊ(0) results in
B11
Using Schur complements (Gallier, 2019), πΊ(0) β₯ π(0) if and only if
B12
It is assumed that πΊπ,0 can be determined so that (πΊπ,0 β ππ,0) β₯ 0. Then it suffices to show that
B13
All matrices in Equation (B13) are diagonal. Denoting an arbitrary diagonal element of πΊπ,0, 
and ππ,0 as 
, 
and π2, respectively, we can therefore consider the scalar condition
B14
It can be shown from Equation (20) that
B15
Substituting this result into Equation (B14) and simplifying yields
B16
With 
and π β [πmin,πmax], the initial variance 
must satisfy the inequality
B17
APPENDIX C: DISCRETE-TIME COVARIANCE ERROR MATRIX
This appendix derives propagation equations for the covariance error matrix Ξ associated with the discrete-time system
C18
where 
and ππ is a vector of first-order GMPs with unknown but bounded time constants. We assume 
and 
can be spec-ified, and that 
is estimated with a KF that models uncertain GMPs using Equation (28), restated below
C19
With π = πβΞπ‘βπ and 
, a discrete-time GMP is governed by the difference equation Rogers (2003) ππ+ 1 = πππ +πΎπ€π, such that π€π βΌ WGN(0, 1). Therefore,
C20
For the process in Equation (C19), let us also define 
and 
so that 
is ππ populated with 
and 
is ππ populated with 
. After defining
C21
the KF covariance matrix propagates according to Gelb et al. (1974)
C22a
C22b
where ππ is the Kalman gain matrix and πΊπ,0 is a diagonal matrix populated with the 
βs from Equation (C19).
Following the same approach as in Section 4, it can be shown that the true covariance matrix propagates as
C23a
C23b
C23c
such that 
, 
and ππ₯,π has the same form as 
, except that 
and 
are replaced with ππ and ππ. A comparison between the KF and true covariance matrices is made using the same approach as in Section 5. That is, by specifying a propagation structure for the larger covariance matrix πΊ that facilitates forming the error matrix π«π|π = πΊπ|π βππ|π. To this aim, we define propagation equations for πΊπ|π as
C24a
C24b
C24c
where 
is defined in Section 5.2 and 
is diagonal with πth diagonal element 
. The reader should be aware that these are not the only options for 
and 
. Other matrices could be specified. We will focus on the particular realizations above because they are a logical choice and facilitate subsequent analysis.
The structure in Equations (C24b) and (C24c) ensures that 1) Equations (C22a) and (C22b) remain intact, 2) πΊ0|0 β₯ π0|0, which was shown in Section 5.2, and 3) the lower right block of πΊπ|π is time-invariant and always equal to its initial value 
. The next step is to form the differences πΊπ+ 1|π βππ+ 1|π and πΊπ|π βππ|π to define propagation equations for the covariance error matrix π«. However, Equation (C24b) is not ideal because the matrix bracketing πΊπ|π on the right-hand side does not match the matrix bracketing ππ|π in Equation (C23b). In response, let us rewrite Equation (C24b) as
C25
where we utilize the fact that πΊ is a block diagonal matrix and that the lower right block is time invariant. Now we can subtract Equations (C23b) and (C23c) from Equations (C24b) and (C24c), resulting in the covariance error propagation equations
C26
C27
APPENDIX D: PROOF OF POSITIVE SEMI-DEFINITENESS
This appendix shows that the following matrix
D28
is positive semi-definite. Recalling that 
and 
, and using the definition for 
given in Equation (C21), π can be written in the 3Γ3 block form
D29
where we have used the fact that ππ has the same structure as 
with 
, 
replaced by ππ and ππ.
Making the definitions
D30a
D30b
allows π to be written as 
. To prove that π β₯ 0, we use the fact that with π > 0 and invertible, π β₯ 0 if and only if (Gallier, 2019)
D31
To show that π > 0 and invertible, first note from Appendix C that 
is a diagonal matrix with πth diagonal element
D32
Since 0 β€ ππ β€ 1 and 
. Thus, the necessary conditions that π > 0 and invertible are satisfied.
Substituting Equations (D30a) and (D30b) into Equation (D31) results in
D33
One way to prove that π β₯ 0 is to show that 
and π² are both positive semi-definite. Since we assumed that 
, clearly 
. It is more difficult to assess whether π² β₯ 0. To proceed, first notice that π² is a diagonal matrix because it is a linear combination of products of diagonal matrices. Thus, we can prove that π² β₯ 0 by showing that the diagonal elements are always nonnegative.
D.1 Defining an arbitrary diagonal element of Ξ
This subsection develops an expression for an arbitrary diagonal element of π². From Appendices B and C, the πth diagonal elements of each matrix in Equation (D33) are
D34
Substituting the expressions from Equation (D34) into Equation (D33), the πth diagonal element of π² is
D35
As 
approaches 
gets smaller. Thus, we will focus on proving that ππ β₯ 0 for the limiting case when 
. From this point forward, the π index will be dropped to simplify notation. Let π = πmaxβπmin and π½ = π2 πmaxβ(π β πmax). Then Equation (D34) shows that 
can be written as π2 = β2π½, and Equation (D35) simplifies to
D36
Defining the expression inside the square brackets in Equation (36) as βπ, we obtain the final form for π used in this appendix
D37
D.2 Monotonicity of π
We will now show that π is a monotonically increasing function of Ξπ‘ for any π, πmin and πmax satisfying π β [πmin,πmax] and πmin β€ πmax. With this fact established, the range of π(Ξπ‘) is entirely defined by its values at Ξπ‘ = 0 and Ξπ‘ = β. The proof is based on the fact that a differentiable function is monotonically increasing if its derivative is nonnegative at every point in the functionβs domain (Bermant, 1963). It is safe to use this criteria for monotonicity because π(Ξπ‘) is continuous and therefore differentiable.
From Equation (D37), we have the following derivatives
D38
which, after simplification, lead to the following expression for ππβπΞπ‘
D39
where πΎ = π2 +2π½ = π2(π + πmax)β(π β πmax). It will be useful to define the quantity in square brackets as
D40
so that ππβπΞπ‘ = β(2πΎπβπ)π.
The factor β(2πΎπβπ) is nonnegative, which follows from the fact that πΎ β€ 0. Therefore, ππβπΞπ‘ will be nonnegative as long as π β₯ 0. We can prove that π is nonnegative by examining its behavior as a function of Ξπ‘.Agiven value Ξπ‘β uniquely determines the quantities 
, 
, 
and πβ that, when substituted into Equation (D40), produce the value πβ. In addition to giving us πβ, Equation (D40) also tells us that πβ must be a root of the equation 
. A necessary condition for the quadratic equation to have at least one real root is that the discriminant 
is greater than or equal to zero (Manfrino et al., 2015), that is, if
D41
From the definitions in Equation (D40), notice that
D42
The term in brackets is always less than or equal to zero, which can be shown by considering the inequality
D43
Substituting the definitions for π½, π and πΎ into Equation (D43) results in the new inequality
D44
Equation (D44) can also be written as
D45
As a function of π, the left-hand side of Equatiion (D45) is a parabola that starts at zero when π = πmin, rises to its maximum value (πmax βπmin)2β4 when π = (πmin + πmax)β2, and then falls back to zero when π = πmax. Thus, the inequality in Equation (D45) is always true. Tracing this result back to Equation (D42), we can conclude that
D46
Given that Equation (D46) is valid for any π, πmin,πmax and Ξπ‘, we can conclude from Equation (D41) that πβ β₯ 0. Since Ξπ‘β is arbitrary, it must be that π is nonnegative for all values of π, πmin,πmax and Ξπ‘. This completes the proof that ππβπΞπ‘ β₯ 0 and therefore that π is a monotonically increasing function of Ξπ‘.
D.3 Establishing the range of π
Since π is a monotonically increasing function of Ξπ‘, the range of π is entirely defined by its value at Ξπ‘ = 0 and Ξπ‘ = β. Substituting these values for Ξπ‘ into Equation (D37) leads to the conclusion that
D47
The sign of the right endpoint is dictated by the sign of the quantity πβ2πmin +πmax. The minimum value of this quantity is πmax βπmin when π = πmin. Since πmin β€ πmax, the right endpoint is greater than or equal to zero.
We have therefore established that π β₯ 0 for all admissible values of π2,π,πmin,πmax and Ξπ‘. This completes the proof that the matrix in Equation (D28) is positive semidefinite.
Footnotes
Disclosure
The authorβs affiliation with The MITRE Corporation is provided for identification purposes only and is not intended to convey or imply MITREβs concurrence with, or support for, the positions, opinions, or viewpoints expressed by the author.
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
- β΅Anderson, B. D. O., & Moore, J. B. (1974). Optimal Filtering (pp. 296-303). Dover Publications.
- β΅Bermant, A. F. (1963). A Course in Mathematical Analysis: Part I, International Series of Monographs on Pure and Applied Mathematics, Vol. 44. Oxford, United Kingdom: Pergamon.
- β΅Blanch, J., Walter, T., & Enge, P. (2019). Gaussian bounds of sample distributions for integrity analysis. IEEE Transactions on Aerosp. and Electronic Systems, 55(4), 1806β1815. https://doi.org/10.1109/TAES.2018.2876583
- Brown, R. G., & Hwang, P. Y. C. (2012). Introduction to Random Signals and Applied Kalman Filtering. 4th ed. Hoboken, NJ: John Wiley & Sons.
- β΅Crassidis, J., & Junkins, J. (2004). Optimal Estimation of Dynamic Systems (pp. 123β174). Chapman & Hall/CRC.
- β΅DeCleene, B. (2000). Defining pseudorange integrity - overbounding. Proc. of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation, Salt Lake City, UT, 1916β1924.
- β΅Dieci, L., & Eirola, T. (1994). Positive definiteness in the solution of Riccati differential equations. Numerische Mathematik, 67(3), 303-313. https://doi.org/10.1007/s002110050030
- β΅Gallier, J. (2019). The Schur complement and symmetric positive semidefinite (and definite) matrices. https://www.researchgate.net/publication/303969420_Notes_on_the_Schur_Complement
- β΅Gelb, A., Kasper, J. F. Jr, Nash, R. A. Jr, Price, C. F., & Sutherland, A. A. Jr. (1974). Applied Optimal Estimation. The M.I.T. Press.
- β΅Joerger, M. (2009). Carrier phase GPS augmentation using laser scanners and using low Earth orbiting satellites. (Ph.D. dissertation). Illinois Institute of Technology.
- β΅Joerger, M., & Pervan, B. (2013). Kalman filter-based integrity monitoring against sensor faults. AIAA J Guid Control Dynam, 36(2), 349-361. https://doi.org/10.2514/1.59480
- β΅Joerger, M., & Pervan, B. (2015). Multi-constellation ARAIM exploiting satellite geometry change. Proc. of the 28th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2015), Tampa, FL, 2688β2704.
- β΅Joerger, M., & Pervan, B. (2016). Exploiting satellite motion in ARAIM: measurement error model refinement using experimental data. Proc. of the 29th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2016), Portland, OR, 1696β1712. https://doi.org/10.33012/2016.14552
- β΅Joerger, M., & Pervan, B. (2020). Multi-constellation ARAIM exploiting satellite geometry change. NAVIGATION, 67, 235-253. https://doi.org/10.1002/navi.334
- β΅Langel, S., GarcΓa Crespillo, O., & Joerger, M. (2019). Bounding sequential estimation errors due to Gauss-Markov noise with uncertain parameters. Proc. of the 32nd International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2019), Miami, FL, 3079β3098. https://doi.org/10.33012/2019.17014
- β΅Langel, S., Khanafseh, S., & Pervan, B. (2014). Bounding integrity risk for sequential state estimators with stochastic modeling uncertainty. Journal of Guidance Control Dynamics, 37(1), 36β46. https://doi.org/10.2514/1.62056
- β΅Pervan, B., Khanafseh, S., & Patel, J. (2017). Test statistic auto- and cross-correlation effects on monitor false alert and missed detection probabilities. Proc. of the 2017 International Technical Meeting of the Institute of Navigation, Monterey, CA, 562β590. https://doi.org/10.33012/2017.14874
- β΅Pulford, G. W. (2008). A proof of the spherically symmetric overbounding theorem for linear systems. NAVIGATION, 55, 283β292. https://doi.org/10.1002/j.2161-4296.2008.tb00437.x
- β΅Rife, J., & Gebre-Egziabher, D. (2007). Symmetric overbounding of correlated errors. NAVIGATION, 54, 109β124. https://doi.org/10.1002/j.2161-4296.2007.tb00398.x
- β΅Rife, J., Pullen, S., Enge, P., & Pervan, B. (2006). Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Transactions on Aerosp. and Electronic Systems, 42(4), 1386β1395. https://doi.org/10.1109/TAES.2006.314579
- β΅Rife, J., Walter, T., & Blanch, J. (2004). Overbounding SBAS and GBAS error distributions with excess-mass functions. Proc. of the 2004 International Symposium on GNSS/GPS, Sydney, Australia.
- β΅Rogers, R. M. (2003). Applied Mathematics in Integrated Navigation Systems. 2nd ed. (pp. 106-109). Reston, VA: American Institute of Aeronautics & Astronautics, Inc.
- β΅Simon, D. (2006). Optimal State Estimation: Kalman, H Infinity and Nonlinear Approaches. John Wiley & Sons.
- β΅Tupysev, V. A., Stepanov, O. A., Loparev, A. V., & Litvinenko, Y. A. (2009). Guaranteed estimation in the problems of navigation information processing. Proc. of the 18th IEEE International Conference on Control Applications, St. Petersburg, Russia, 1672β1677. https://doi.org/10.1109/CCA.2009.5281081
- β΅Working Group C - ARAIM Technical Subgroup. Interim report issue 1.0. EU-U.S. Cooperation on Satellite Navigation; December 2012.
- β΅Working Group C - ARAIM Technical Subgroup. Milestone 2 report. EU-U.S. Cooperation on Satellite Navigation; February 2015.
- β΅Working Group C - ARAIM Technical Subgroup. Milestone 3 report. EU-U.S. Cooperation on Satellite Navigation; February 2016.
- Tanil, C., Khanafseh, S., Joerger, M., & Pervan, S. (2018). 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
- β΅Manfrino, R. B., Ortega, J. A. G., & Delgado, R. V. (2015). Topics in Algebra and Analysis: Preparing for the Mathematical Olympiad. New York, NY: Springer International Publishing.
