Overbounding the effect of uncertain Gauss-Markov noise in Kalman filtering

Prior work established a model for uncertain Gauss-Markov (GM) noise that is guaranteed to produce a Kalman filter (KF) covariance matrix that over-bounds 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  2 0 , and propagate according to the difference equation ∼ WGN (0, 1) and  0 ∼ (0,  2 0 ), 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,  2 0 ) denotes a zero-mean normal random variable with variance  2 0 .
When  and  2 are only known to reside in the intervals [ min ,  max ] and [0,  2 max ], a state-augmented KF that models   with  =  max  2 =  2 max ( max ∕ min ) is guaranteed to produce a covariance matrix  ≥ .

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 timecorrelated 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 continuoustime 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 continuousdiscrete 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.

PROBLEM SETUP
Consider the state estimation problem for the sampleddata system ξ = () + ()() 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 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 N() ≥ () and R ≥   can be determined.
For the GMP vector ,  and  are (Gelb et al., 1974) where  2  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  2  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 With   = [    ], Equation (4) can also be written as Equation ( 5) is a linear system driven by zero-mean WGN.Thus it is permissible to estimate  using a KF.Since  2  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 F I G U R E 1 One-dimensional position and velocity estimation 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 zeromean WGN and the fact that the KF is a linear unbiased estimator, it follows that   ∼ (0,  2  ).In this case, an upper bound on integrity risk is equivalent to an upper bound on the estimate error variance  2  .Therefore, subsequent sections of the paper will focus on developing an upper bound on the estimate error variance for any state of interest.

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 where   is the sum of a first-order GMP and zero-mean WGN.That is,   =   +   such that   ∼ (0,  2  ) and The variances are known and given by  2  = 0.5m 2 and  2  = 1m 2 .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 F I G U R E 2 Difference between predicted and true estimate error standard deviations for a KF that assumes such that the initial estimate error covariance matrix  0 is diagonal with nonzero elements: [10 m 2 , 1 m 2 ∕s 2 , 1 m 2 ].
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.

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 = P|−1  The filter designer must specify L (which is contained in F) 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 P therefore does not accurately capture the true estimate error distribution, and a new covariance matrix  must be defined that accounts for model uncertainty.
When L ≠ , Appendix A shows that the KF estimate error vector  =  − x propagates according to ε = F + Δ + , where Δ  = [ ( − L)  ].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 Equations ( 10) and ( 11) contain both true and estimator matrices and thus reflect that , ,  and  are uncertain.For example, N and R 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.

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, This structure was chosen because, as we will see, it enables F, Q and the initial covariance matrix (0) to be defined explicitly from the requirement ( − ) ≥ 0. The covariance error matrix  =  −  propagates as Equation ( 13a) is a matrix differential equation with initial condition  −1|−1 and whose solution at Dieci & Eirola, 1994).The quadratic structure in Equation ( 13b) ensures that  | ≥ 0 when  |−1 ≥ 0 and R ≥   .Therefore, we must have an initial matrix (0) ≥ 0 and Q ≥  to guarantee that () ≥ 0 for all .

Continuous-time model
We it is clear that ( Q − ) ≥ 0 for all admissible values of  2  and   .However, Q could never be constructed because it depends on the unknown matrix .Therefore, we must use a mathematical approach to make  insensitive to .
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 Δ =  − L and impose the constraint which simplifies Equation (15b) to the unforced linear differential equation Σ = F  +     .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 Σ =  and   = −(Δ) −1  into Equation (15c) results in the additional constraint Noting that , L,  and  are diagonal matrices, Equation ( 18) is equivalent to the set of scalar conditions With τ =  ,max , the inequality in Equation ( 20) simplifies to 1 ≥   ∕ ,max , which is satisfied for all admissible values of τ .Thus, we satisfy the requirement that  2  ≥  2  .Having determined the diagonal elements of L and Û, we can conclude that   ≥   when first-order GMPs are modeled according to 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  2  .It is comforting that  2  is not required to propagate   because its dependence on the unknown time constant   means that  2  can never be constructed.However, a previously unknown fact is that  2  places constraints on the initial variance for the GM model.This new result is covered next.

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 where  ,0 is assumed to be known (e.g., from an initialization algorithm),  ,0 is what we seek to define in this section, and  ,0 is a diagonal matrix populated with the  2  in Equation ( 20).Appendix B shows that the true initial covariance matrix is such that  ,0 is a diagonal matrix comprised of the unknown variances  2  .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 Denoting an arbitrary diagonal element of  ,0 as  2 0 , Appendix B shows that Equation ( 24) is satisfied if It is worth noting that the model in Equation ( 21) is stationary when initialized with a variance  2 max ( max ∕ min ) 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.

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) where  [     ] =    kl and  [     ] =    kl .As before, we assume that N ≥   , R ≥   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) with  =  −Δ∕ , the discretized form of Equation ( 21) is seen to be 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 F I G U R E 3 KF standard deviation compared to truth when using stationary and non-stationary bounding GMP models 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.

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  2 0 =  2 max ( max ∕ min ).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 True and bounding covariance ellipses for the motivational example in Section 3. Ellipses are shown at an elapsed time of 10 seconds less conservative.That is, the gap between the predicted and true standard deviations is smaller when using the non-stationary model for GM noise.
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 nonstationary 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.
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 nonstationary model in an example application of aircraft nav-igation using ARAIM (Working Group C -ARAIM Technical Subgroup, 2012Subgroup, , 2015Subgroup, , & 2016)).

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, 2015Subgroup, & 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).

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, timeinvariant, 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 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 R 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 where  2  ,  2  and  2  are defined in Table 1.The values of   ,   and   are unknown but bounded following the inequalities 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, 2015Subgroup, & 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.

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: 1. an elevation mask of 5 degrees 2. an integrity risk requirement of 10 −7 per approach 3. a continuity risk requirement of 4 × 10 −6 per approach 4. a vertical alert limit  of 35 m (unless otherwise stated) 5. a prior probability of satellite fault of 10 −5 All other simulation parameters are given in Joerger and Pervan (2020).
F I G U R E 5 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 6. a prior probability of constellation fault  const of 10 −4 (unless otherwise stated) 7. 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 .
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 F I G U R E 6 Availability map for snapshot ARAIM using depleted constellations,  const = 10 −4 ,  = 35 m (coverage of 99.5% availability is 84.6%) 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 F I G U R E 7 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%) cosine of the location's latitude, because grid point locations near the equator represent larger areas than near the poles.
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.

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.

A C K N O W L E D G E M E N T S
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.

R E F E R E N C E S
The covariance matrix  =  [  ] propagates between measurements according to Simon ( 2006) where During a measurement update, the covariance matrix is updated through 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 â0 = .We also assume that ξ0 and â0 are initially uncorrelated, so that  [ ,0   ,0 ] = .This leads to the following expected values where we have used the fact that because â0 = , [ ,0 â 0 ], [ 0 â 0 ] and  [ â0 â 0 ] are all equal to .Using these relations, the true initial covariance matrix of The requirement   =  implies that The goal is to determine  ,0 so that (0) ≥ (0).Subtracting (0) from (0) results in ] ≥ 0 (B12) It is assumed that  ,0 can be determined so that ( ,0 −  ,0 ) ≥ 0. Then it suffices to show that ( ,0 −  ,0 ) −  ,0 ( ,0 −  ,0 ) All matrices in Equation (B13) are diagonal.Denoting an arbitrary diagonal element of  ,0 ,  ,0 and  ,0 as  2 0 ,  2 and  2 , respectively, we can therefore consider the scalar condition It can be shown from Equation (20) that Substituting this result into Equation (B14) and simplifying yields With  2 ∈ [0,  ) ( With  =  −Δ∕ and  =  √ 1 −  2 , a discrete-time GMP is governed by the difference equation Rogers ( 2003)  +1 =   +   , such that   ∼ WGN(0, 1).Therefore, For the process in Equation (C19), let us also define φ =  where   is the Kalman gain matrix and  ,0 is a diagonal matrix populated with the  2 0 's from Equation (C19).Following the same approach as in Section 4, it can be shown that the true covariance matrix propagates as such that Δ   = [0 ΔL   ],    = [    ] and  , has the same form as Q, , except that N 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 where  ,0 is defined in Section 5.2 and   is diagonal with  th diagonal element  2  (1 −  2  ).The reader should be aware that these are not the only options for  ,0 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  ,0 .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 One way to prove that  ≥ 0 is to show that   ( N −   )   and  are both positive semi-definite.Since we assumed that N ≥   , clearly   ( N −   )   ≥ 0. 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 Substituting the expressions from Equation (D34) into Equation (D33), the  th diagonal element of  is As  2  approaches  2 ,max ,   gets smaller.Thus, we will focus on proving that   ≥ 0 for the limiting case when  2  =  2 ,max .From this point forward, the  index will be dropped to simplify notation.Let  =  max ∕ min and  =  2  max ∕( −  max ).Then Equation (D34) shows that

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 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 Δ.A given value Δ * uniquely determines the quantities  * 0 ,  * 1 ,  * 2 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  * 2  2 +  * 1  +  * 0 −  * = 0.A necessary condition for the quadratic equation to have at least one real root is that the discriminant  = ( * 1 ) 2 − 4 * 2 ( * 0 −  * ) is greater than or equal to zero (Manfrino et al., 2015), that is, if From the definitions in Equation (D40), notice that The term in brackets is always less than or equal to zero, which can be shown by considering the inequality 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 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 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.
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 Q, readers should re-familiarize themselves with the definitions in Equation (3).Let Û be  formed with the upper bound values for  2  and the lower bound values for   .Let's also define  as the matrix  populated with the true time constants   and as yet undetermined variances  2  ≥  2  .Then with Q defined as
28)Following a similar approach to Subsection 5.1, Appendix C shows that the covariance error matrix  =  −  for the discrete-time KF propagates as kl , [     ] =     and   is a vector of first-order GMPs with unknown but bounded time constants.We assume N ≥   and R ≥   can be spec- −Δ∕ max and γ =  ,| = ( −     ) ,|−1 ( −     ) ], and using the definition for Q given in Equation (C21),  can be written in the 3 × 3 block form used the fact that   has the same structure as Q with N , Û replaced by   and   .  ≤ 1 and   ≤  ,max , (  −   )  > 0.