Measurement Error Time-Correlation Modeling for Safety-Critical Navigation

2.1 Background on Time-Correlation Modeling in the Time Domain

We introduce the time-domain modeling method reported by Langel (2014) using the following simplified illustrative example of an unbiased, linear, scalar state estimation process using two scalar measurements arranged in a batch observation vector. The scalar state estimation error ε_τ at time τ is given as follows:

${\epsilon }_{\tau }=\left(\begin{array}{cc}{a}_0& a_{\tau }\end{array}\right)\left(\begin{array}{c}{v}_0\\ v_{\tau }\end{array}\right)=a_0{v}_0+a_{\tau }{v}_{\tau }$ 1

where a₀ and a_τ are the estimator coefficients (derived in Appendix F for an example least-squares estimator) and v₀ and v_τ are the measurement errors at times 0 and τ, respectively. The measurement errors are outputs of a zero-mean wide-sense stationary (WSS) process. This process has a known variance $E\left(v_0^2\right)=E\left(v_{\tau }^2\right)={\sigma }_v^2$ and an uncertain ACF r(τ), with $r\left(\tau \right)=E\left(v_0{v}_{\tau }\right)$, where E[.] is the expectation operator. The estimation error variance ${\sigma }_{\epsilon }^2$ at time τ can be expressed as follows:

${\sigma }_{\epsilon }^2={\sigma }_v^2\left(a_0^2+a_{\tau }^2\right)+2a_0{a}_{\tau }r\left(\tau \right)$ 2

We can upper-bound ${\sigma }_{\epsilon }^2$ in Equation (2) if bounds on r(τ) are known. For example, a tight upper bound on ${\sigma }_{\epsilon }^2$ is found by using a knowledge of the sign of a₀a_τ and assuming upper and lower bounds on r(τ) (Langel, 2014; Langel et al., 2019). Upper- and lower-bounding ACFs can be derived from FOGMP models and expressed as follows:

$\begin{array}{cccc}{r}_{\max }\left(\tau \right)={\sigma }_{\max }^2\exp \left(-\frac{\tau }{T_{\max }}\right)& \text{and}& r_{{}_{\min }}\left(\tau \right)={\sigma }_{\min }^2\exp \left(-\frac{\tau }{T_{\min }}\right)& \hspace{0.17em}\text{for}\hspace{0.17em}\end{array}\tau \in \left(0,{\tau }_{op}\right)$3

where the notation exp(·) designates an exponential function and τ_op is an operational limit on lag times of interest. The upper- and lower-bounding FOGMP model ACFs in Equation (3) are decaying exponential functions with FOGMP time constants T_max and T_min and variance ${\sigma }_{\max }^2$ and ${\sigma }_{\min }^2$ as parameters, respectively. In practice, the measurement error ACF is estimated from data for lag times τ ranging from zero to τ_op. For example, GPS satellite clock and orbit ephemeris errors do not need to be modeled over a lag-time interval larger than τ_op = 7 h, which is the maximum duration of a satellite pass for a ground user (Gallon et al., 2019; Perea, 2019).

Figure 1 illustrates the paired ACF bounding procedure for 80 simulated sample ACFs drawn from the same FOGMP with unit variance and a time constant of T = 50 s. Lag times are displayed from 0 to 100 s. The 80 ACFs are generated by using 80 independent 100-s-long time series. The 100-s limit captures the fact that the actual error data may be sparse or may need to be partitioned to meet the WSS assumption (Gallon et al., 2022). Even though all 80 curves are drawn from the same FOGMP, variations between sample ACFs are significant because of data sparsity.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 1

Overview of the paired ACF bounding process for a unit-variance FOGMP with a time constant of 50 s using 100-s-long data time series

The dotted and solid black curves in Figure 1 indicate the upper and lower FOGMP ACF bounds obtained with the assumption of unit variance. The gray sample ACFs cannot be fully lower-bounded because they have negative values for lag times τ approaching 35 s. In this example, if the time-correlation model is to be used over τ_op intervals longer than 35 s, a different method is needed.

2.2 Time-Correlation Modeling Using Lagged Products

This section aims to address the limitations of FOGMP ACF bounding by finding a FOGMP model that “overbounds lagged products” for all values of τ belonging to $\left(0\hspace{0.17em}{\tau }_{op}\right)$ and for all quantiles of the lagged product CDF. First, we introduce and define lagged products. Second, we formulate a closed-form expression for the probability density function (PDF) of the FOGMP lagged products. Third, we relate lagged product overbounding to ACF upper- and lower-bounding.

In the first part of this section, we define a lagged product as follows:

$q\left(\tau \right)=v_0{v}_{\tau }$4

Given a lag time τ, the lagged product is a random variable whose mean value is the ACF, r(τ). The distribution of q(τ) can be modeled by using overbounding theory. Widely used in GNSS-based safety critical aviation applications (RTCA Special Committee 159, 2001, 2004; Working Group C - ARAIM Technical Subgroup, 2016), overbounding provides a means to compare and model non-Gaussian probability distributions. In addition, if an error model’s CDF overbounds sample measurement errors, then this model function (e.g., a Gaussian) can be used to predict an overbound on the actual estimation error’s CDF, even if the latter is non-Gaussian (Blanch et al., 2017; DeCleene, 2000; Rife et al., 2006). To determine an overbounding model for q(τ), all quantiles of its distribution must be considered for all values of $\tau \in \left(0,{\tau }_{op}\right)$. Because error modeling is performed offline, the computational load of this two-dimensional approach is not prohibitive.

Figure 2(a) shows 2500 sample lagged product curves for a FOGMP with variance ${\sigma }_v^2=1$ and time constant T = 50 s over lag times ranging from 0 to 100 s. For each value of τ, the q(τ)-sample CDF is determined, with quantile curves shown in Figure 2(b). In contrast with the FOGMP ACF, quantiles of lagged products can have negative values.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 2

(a) 2500 lagged product curves for a FOGMP with unit variance and time constant T = 50 s; (b) sample quantiles of the lagged product distribution

The second part of this section aims at deriving a parametric, closed-form expression of the lagged product distribution of a FOGMP model, which will then be used to overbound the empirical distribution of q(τ) samples. A discrete-time FOGMP v_n at time step n with time constant T can be expressed as follows:

$\begin{array}{c}\begin{array}{ccc}{v}_n=\exp \left(-\frac{\Delta t}{T}\right)v_{n-1}+{\eta }_{n-1}& \text{where}& {\eta }_n\sim N\left(0,{\sigma }_{\eta }^2\right)\end{array}\\ \begin{array}{cc}\text{with}& {\sigma }_{\eta }^2\triangleq {\sigma }_v^2\left(1-\exp \left(-\frac{2\Delta t}{T}\right)\right)\end{array}\end{array}$5

We use the notation $v_n\sim N\left(0,{\sigma }_v^2\right)$ to designate the fact that v_n is normally distributed with zero mean and variance ${\sigma }_v^2$.

We use the formulas from Cui et al. (2016) to derive the following closed-form expression of the FOGMP lagged product PDF:

$\begin{array}{c}{f}_{q\left(\tau \right)}\left(x\right)=\frac{1}{\pi {\sigma }_v^2\sqrt{1-{\alpha }^2}}\text{exp}\left(\frac{x\alpha }{{\sigma }_v^2\sqrt{1-{\alpha }^2}}\right)K_0\left(\frac{\left(x\right)}{{\sigma }_v^2\sqrt{1-{\alpha }^2}}\right)\\ \begin{array}{cc}\hspace{0.17em}\text{with}\hspace{0.17em}& \alpha \triangleq \text{exp}\end{array}\left(-\frac{\tau }{T}\right)\end{array}$6

where $k_0\left(.\right)$ is the modified Bessel function of the second kind of order zero. When $\tau =0$, the lagged product is $v_0^2$, and as expected, Equation (6) simplifies to a one-degree-of-freedom chi-square distribution scaled by ${\sigma }_v^2$ (as shown in Appendix D). As the lag time approaches infinity, the correlation between v₀ and v_τ approaches zero, and the lagged product tends to be distributed according to a modified Bessel function of the second kind (symmetric around zero).

Equation (6) is used in Figure 3(a) to represent the theoretical PDF of a FOGMP with T = 50 s and ${\sigma }_v^2=1$. We numerically integrate this PDF to obtain the model lagged product CDF for each lag time of interest. Lines of constant probability are displayed in Figure 3(b) and represent the theoretical quantile lines corresponding to the sample quantiles in Figure 2.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 3

(a) PDF of FOGMP lagged products; (b) CDF quantiles of a FOGMP with unit variance and time constant T = 50 s

The third part of this section shows that overbounding the FOGMP lagged product CDF can guarantee upper and lower ACF bounds. We apply the following theorem derived in Appendix B and first used by Jada and Joerger (2020):

Theorem: For an arbitrary CDF $F_A\left(x\right)$ with mean ${\mu }_A$ and for any CDF upper-bounding and lower-bounding functions F_U(x) and F_L(x) with mean values μ_U and μ_L, respectively, such that:

$\begin{array}{ccc}{F}_U\left(x\right)\ge F_A\left(x\right)\hspace{0.17em}\forall x& \text{and}& F_L\left(x\right)\le F_A\left(x\right)\hspace{0.17em}\forall x\end{array}$7

the following inequalities are always satisfied:

$\begin{array}{ccc}{\mu }_U\le {\mu }_A& \text{and}& {\mu }_L\ge {\mu }_A\end{array}$8

Applied to lagged product overbounding, the mean of the CDF upper-bounding FOGMP, ${\sigma }_{\min }^2\exp \left(-\tau /T_{\min }\right)$, provides a lower bound for the mean of the actual distribution, $E\left(q\left(\tau \right)\right);$ the mean of the CDF lower-bounding FOGMP, ${\sigma }_{\min }^2\exp \left(-\tau /T_{\min }\right)$, provides an upper bound for $E\left(q\left(\tau \right)\right)$. Thus, the time-correlation models in Equation (3) are determined by finding the pair of FOGMPs whose lagged product CDFs upper- and lower-bound the sample CDF for all τ values at all quantiles.

2.3 Step-by-Step Procedure for Time-Correlation Modeling in the Time Domain

Consider a measurement error data set with LN samples partitioned into L time series of N samples each and expressed in matrix form as follows:

$V\triangleq \begin{array}{ccc}\left(\begin{array}{ccc}{v}_1\left({\tau }_0\right)& \cdots & v_1\left({\tau }_{N-1}\right)\\ ⋮& \ddots & ⋮\\ v_L\left({\tau }_0\right)& \cdots & v_L\left({\tau }_{N-1}\right)\end{array}\right),& \text{with}& N=\frac{{\tau }_{op}}{\Delta t}\end{array}$9

where ${\tau }_n=n\Delta t$ for n = 0, …, N – 1, Δt is the sampling interval, and τ_op is the required model validity period. Each row of V is a time series over τ_op. The lag τ_op is greater than two time constants, ensuring that the time series include uncorrelated samples, i.e., separated by a time interval larger than or equal to τ_op.

We multiply each row by its first element to obtain the following lagged product matrix:

$\begin{array}{ccccc}Q\triangleq \left(\begin{array}{ccc}{q}_1\left({\tau }_0\right)& \cdots & q_1\left({\tau }_{N-1}\right)\\ ⋮& \ddots & ⋮\\ q_L\left({\tau }_0\right)& \cdots & q_L\left({\tau }_{N-1}\right)\end{array}\right),& \text{with}& q_1\left({\tau }_n\right)\triangleq v_{l,0}{v}_{l,n}& \text{for}& l=1,\hspace{0.17em}\dots \hspace{0.17em},L\end{array}$10

At any given lag time τ_n, we have a column vector of L sample lagged products, ${\left(q_1\left({\tau }_n\right),\hspace{0.17em}\dots \hspace{0.17em},q_L\left({\tau }_n\right)\right)}^T$. Figure 2 presents results for an example Q matrix with N = 100 and L = 2500.

Next, we sort each column of Q in ascending order of its elements and arrange the elements in the following matrix form:

$\begin{array}{ccc}{Q}_{\text{sorted}}\triangleq \left(\begin{array}{ccc}{q}_{\left(1\right)}\left({\tau }_0\right)& \cdots & q_{\left(1\right)}\left({\tau }_{N-1}\right)\\ ⋮& \ddots & ⋮\\ q_{\left(L\right)}\left({\tau }_0\right)& \cdots & q_{\left(L\right)}\left({\tau }_{N-1}\right)\end{array}\right),& \text{Where}& q_{\left(1\right)}\left({\tau }_n\right)\end{array}\le q_{\left(l+1\right)}\left({\tau }_n\right)$11

The subscript notation (·) designates sorted indices. Let ${\left(F_{1,n\dots }{F}_{L,n}\right)}^T$ be the L × 1 empirical CDF vector derived from ${[q_{\left(1\right)}\left({\tau }_n\right)\hspace{0.17em}\dots \hspace{0.17em}{q}_{\left(L\right)}\left({\tau }_n\right)]}^T$ for lag time τ_n. The empirical CDF vectors at all N lags are captured in the following matrix form:

$F\triangleq \left(\begin{array}{ccc}{F}_{1,0}& & F_{1,N-1}\\ ⋮& ⋮& ⋮\\ F_{L,0}& & F_{L,N-1}\end{array}\right)$12

We find the FOGMP model parameters T_min, σ_min, T_max, and σ_max by ensuring that the upper- and lower-bounding model CDFs $F_{l.n}\left(\begin{array}{cc}{T}_{\min },& {\sigma }_{\min }\end{array}\right)$ and $F_{l.n}\left(\begin{array}{cc}{T}_{\max ,}& {\sigma }_{\max }\end{array}\right)$ satisfy the following inequalities:

$\begin{array}{c}{F}_{l,n}\left(T_{\text{max}},{\sigma }_{\text{max}}\right)\le F_{l,n}\le F_{l,n}\left(T_{\text{min}},{\sigma }_{\text{min}}\right),\\ \begin{array}{ccc}\text{for}\hspace{0.17em}l=1,\mathrm{...},L& \text{and}& n=0,\mathrm{...},N-1\end{array}\end{array}$13

Figure 4(a) shows 2000 gray lagged product curves for 100-s-long time series of a simulated FOGMP with T = 50 s and σ = 1. The top left-hand-side chart displays $F_{l.n}\left(\begin{array}{cc}{T}_{\max ,}& {\sigma }_{\max }\end{array}\right)$ with blue dashed lines bounding four example sample quantile curves at 2%, 16%, 84%, and 98% to determine T_max = 90 s and σ_max = 1.05. The bottom left-hand-side chart displays $F_{l.n}\left(\begin{array}{cc}{T}_{\min },& {\sigma }_{\min }\end{array}\right)$ (green dashed-dotted line) lower-bounding quantile curves to determine T_min = 10 s and σ_min = 0.96. Figure 4(b) shows q(τ)-CDF upper and lower bounds at example lag times of τ = 1 s, τ = 20 s, and τ = 50 s.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 4

Simulated data with upper- and lower-bounding FOGMP models for (a) example quantiles over all lag times and (b) example lag times at all quantiles

For a more complete visualization, Figures 5(a) and 5(b) show the differences between the sample CDF and the upper- and lower-bounding model CDFs, respectively, for all lag times and quantiles. The color code designates the tightness of the bound, with a tight bound indicated by blue and a loose bound indicated by yellow. If the pair of model functions failed to bound the empirical CDF, the CDF difference would be shown in red. In this example, the CDF bounds are loosest at the core of the distributions (for sample quantiles ranging from 5% to 95%) for τ ≥ 30 s for $F_{l,n}\left(T_{\max },{\sigma }_{\max }\right)$ and for τ ≥ 30 s for $F_{l,n}\left(T_{\min },{\sigma }_{\min }\right)$. This finding is consistent with well-established snapshot overbounding methods reported by Blanch et al. (2017), DeCleene (2000), and Rife et al. (2006): bounding a fat-tailed distribution via a low-order parametric model captures the distribution tails at the cost of a loose bound of the core. Loose measurement error bounds for some quantiles can cause loose estimation error bounds (Crespillo et al., 2023; Joerger et al., 2023; Langel et al., 2020), but this trade-off is necessary to achieve practical, computationally efficient, and memory-efficient error modeling in high-integrity estimation.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 5

Contour plots showing the color-coded difference between the model CDF and sample CDF for the (a) lower-bounding FOGMP model CDF and (b) upper-bounding FOGMP model CDF

Positive values throughout indicate that the models are bounding at all lag times and CDF quantiles.

Model determination using a previously recorded data set that is representative of operating conditions is performed offline and only needs to be performed once for a given data set. In this context, we visually tuned the FOGMP models in Figures 4 and 5 by iteratively modifying the FOGMP parameters, starting from an initial guess and adjusting the model until its distribution bounds the sample distribution while reducing the yellow-shaded areas in Figure 5. We selected the values of T_min, σ_min or T_max, σ_max to achieve as tight a bound as possible. Further model parameter optimization is relevant and is an active area of research for PSD bounding (Crespillo et al., 2023; Gallon & Veiga, 2024; Joerger et al., 2023; Langel et al., 2024), but is beyond the scope of this paper.

3.1 Background on Time-Correlation Modeling in the Frequency Domain

With time-domain modeling methods, there is no generic approach to recursively compute a bounding estimation error variance using parametric models of non-FOGMP time-correlated measurement errors. The method reported by Langel (2014), which is illustrated in Equation (2), requires that all individual estimator coefficients over time (a₀ and a_τ in Equation (1)) be stored and processed at each time step. For recursive, infinite-horizon estimators such as KFs, the number of estimator coefficients over time increases without bound. Therefore, storing estimator coefficients over time, as done by Langel (2014), is both computationally expensive and memory-expensive and can be intractable in long-duration, high-sampling rate implementations (Crespillo et al., 2023).

Frequency-domain time-correlation modeling provides a practical solution for estimation error modeling in recursive implementations (Gallon et al., 2022; Langel et al., 2024). We use the following illustrative scalar example from Gallon et al. (2021) to introduce PSD upper-bounding. A scalar output deviation ε(t) is estimated using a steady-state KF with scalar zero-mean stationary noise input measurement deviations v(t).

The real-valued measurement error PSD s(ω) is a periodic function of the circular frequency $\omega \in \left(0,\pi \right)$ . The circular frequency ω, in units of radians, is related to the frequency f in Hertz through $\omega =2\pi f\Delta t$, where Δ(t) is the sampling interval (Chatfield & Xing, 2019, Chapter 6). If the KF is designed assuming $\overline{s}\left(\omega \right)$ and if the KF transfer function from v(t) to ε(t) is H(ω), then the mean of the KF output error is zero, and the true and predicted output error variances, respectively, are as follows:

$\begin{array}{ccc}{\sigma }_{\epsilon }^2={\int }_0^{\pi }{\left(H\left(\omega \right)\right)}^{2}s\left(\omega \right)d\omega & and& {\overline{\sigma }}_{\epsilon }^2={\int }_0^{\pi }{\left(H\left(\omega \right)\right)}^2\overline{s}\left(\omega \right)d\omega \end{array}$14

If $\overline{s}\left(\omega \right)\underset{\_}{>}s\left(\omega \right)$ for all $\omega \in \left(0\pi \right)$, then ${\sigma }^{\begin{array}{c}2\\ \end{array}}{}_{\epsilon }\underset{\_}{>}{\sigma }^{\begin{array}{c}2\\ \end{array}}{}_{\epsilon }$. The PSD upper bound $s¯\left(\omega \right)$ can be used to define a time-correlated measurement error model that guarantees an upper bound on the KF estimation error variance (Gallon et al., 2021; Langel et al., 2020).

In practice, measurement error PSDs can be estimated from data. The PSD is defined as the discrete-time Fourier transform of the ACF, which can be expressed using the Wiener-Khinchin theorem as follows (Chatfield & Xing, 2019, Chapter 6):

$s\left(\omega \right)\triangleq \frac{1}{\pi }\sum _{n=-\infty }^{\infty }r\left({\tau }_n\right)e^{-i\omega n}$15

where i is the unit imaginary number and r(τ_n) is the ACF at lag time $\left({\tau }_n\right)=n\Delta t$. ACF estimation from actual data requires averaging over multiple time series and using a tapered windowing function as described by Langel et al. (2020). These additional processing steps can introduce additional uncertainty in the PSD estimation process. An alternative is to use periodograms. For a measurement error time series v_n for n = 0, …, N = 1, a periodogram can be defined as follows:

$\begin{array}{ccccc}{\rm P}\left(\omega \right)\triangleq f\left(\omega \right)\overline{f\left(\omega \right)}& \text{with}& f\left(\omega \right)\triangleq \sum _{n=0}^{N-1}{v}_n{e}^{-i\omega n}& \text{for}& \omega \in \left(0,\pi \right)\end{array}$16

where f(ω) is the discrete Fourier transform (DFT)$\left(f\right(\omega )\in ℂ)$. The $\overline{\left(\cdot \right)}$ notation designates complex conjugation. The PSD can then be defined as follows (Chatfield & Xing, 2019, Chapter 7):

$s\left(\omega \right)\equiv \frac{1}{\pi }\underset{N\to \infty }{\lim }\frac{E\left({\rm P}\left(\omega \right)\right)}{N}$17

In Figure 6, the top block diagram displays PSD estimation from the ACF, whereas the bottom diagram represents a periodogram-based PSD estimate.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 6

A block diagram depicting two different PSD estimation schemes

The top branch shows ACF-based PSD estimation, and the bottom shows periodogram-based PSD estimation. The $\left(\hat{\cdot }\right)$ notation is used for estimates, and ${\left(\cdot \right)}^2$ represents the squared magnitude of a complex number.

Gallon et al. (2020), Gallon et al. (2021, 2022), and Langel et al. (2020) generated multiple PSD estimates to account for different satellite clock types when analyzing orbit and clock ephemeris errors and to account for seasonal variations when analyzing tropospheric errors. FOGMP models were derived whose PSDs provide an upper bound for the worst-case sample PSD, which is conservative, but does not account for the probability of occurrence of outlier events.

The following section describes a refined method using CDF overbounding of sample scaled periodograms $\left(\omega \right)/\left(\pi N\right)$ to determine a PSD upper bound. This periodogram distribution is derived from measurement error time series using the bottom block diagram in Figure 6, but without the last block, in which the mean value is taken. A key step in achieving periodogram overbounds is the derivation of a closed-form, parametric expression of FOGMP scaled periodogram distributions.

3.2 Scaled Periodogram Approach to Time-Correlation Modeling

In Appendix A, we show that the scaled periodogram of a FOGMP can be expressed as a linear combination of two independent chi-square random variables $y_1^2$ and $y_2^2$ with a single degree of freedom. Here, we use the notation $y_1^2\sim {\chi }^2\left(1\right)$ and $y_2^2\sim {\chi }^2\left(1\right)$. A scaled periodogram can be expressed as follows:

$\frac{{\rm P}\left(\omega \right)}{\pi N}=\frac{{\lambda }_1^2}{\pi N}{y}_1^2+\frac{{\lambda }_2^2}{\pi N}{y}_2^2$18

where ${\lambda }_1^2$ and ${\lambda }_2^2$ are the two real eigenvalues of the following symmetric positive definite matrix:

$\begin{array}{cc}\Lambda {\Lambda }^T\triangleq \left(\begin{array}{c}\mathrm{Re}{\left(c\right)}^T\\ \mathrm{Im}{\left(c\right)}^T\end{array}\right)\left(\begin{array}{cccc}{\sigma }_v^2& 0& \cdots & 0\\ 0& {\sigma }_{\eta }^2& \cdots & 0\\ ⋮& ⋮& \ddots & ⋮\\ 0& 0& \cdots & {\sigma }_{\eta }^2\end{array}\right)\left(\begin{array}{cc}\mathrm{Re}\left(c\right)& \mathrm{Im}\left(c\right)\end{array}\right)& \in {ℝ}^{2\times 2}\end{array}$ 19

Here, ${\sigma }_v^2$ and T are the variance and time constant of the FOGMP, respectively, ${\sigma }_{{\eta }_{}^{}}^2$ is the variance of the FOGMP driving noise defined in Equation (5), and $c\in {ℂ}^{N\times 1}$ is a vector of complex numbers. We use the notations Re(.) and Im(.) to designate the element-wise real and imaginary parts of a complex vector. The n-th element of c is defined as follows:

$c_n\triangleq \exp \left(-in\omega \right)\frac{1-exp{\left(-\frac{\Delta t}{T}-i\omega \right)}^{N-n}}{1-\exp \left(-\frac{\Delta t}{T}-i\omega \right)},\hspace{0.17em}n=0,\hspace{0.17em}\dots \hspace{0.17em},N-1$20

Thus, ${\lambda }_1^2$ and ${\lambda }_1^2$ depend on T, ${\sigma }_v^2$ Δ(t), N, and the DFT frequency ω. Because ${\lambda }_1^2$ and $y_2^2$ are not necessarily equal, the scaled periodogram in Equation (18) is distributed according to a generalized chi-square distribution. We can use the numerical approach from Davies (1980) to compute the scaled periodogram PDF and CDF.

Figure 7 shows the PDF and CDF of an example FOGMP with T = 50 s and ${\Delta }_v^2=1$ for =t = 1 s and N = 1000. For a sufficiently large number of samples N, the mean of the scaled periodogram approaches the PSD (see Equation (17)). Therefore, using the theorem in Section 2.2, a CDF lower bound on the scaled periodogram guarantees a PSD upper bound. To model the measurement error time correlation, i.e., to obtain a tight upper bound on measurement error PSD, we find the FOGMP whose scaled periodogram CDF gives a lower bound for the sample CDF at all frequencies.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 7

(a) The PDF surface and (b) quantiles of a scaled periodogram of a FOGMP with time constant T = 50 s and variance ${\sigma }_v^2=1$ for a time series of N = 1000 samples with a time step of Δt = 1 s

3.3 Step-by-Step Procedure for Time-Correlation Modeling in the Frequency Domain

Consider the partitioned measurement error matrix V defined in Equation (9). We compute scaled periodograms following the bottom branch in Figure 6 (without the last block) at discrete circular frequency values $w_n\in \left(0,\pi \right)$, for example, for n = 0, ⋯, N – 1 at N logarithmically spaced intervals. The sorted scaled periodograms are expressed in matrix form as follows:

${{\rm P}}_{sorted}\triangleq \frac{1}{N\pi }\left(\begin{array}{ccc}{p}_{\left(1\right)}\left({\omega }_0\right)& \cdots & p_{\left(1\right)}\left({\omega }_{N-1}\right)\\ ⋮& \ddots & ⋮\\ p_{\left(L\right)}\left({\omega }_0\right)& \cdots & p_{\left(L\right)}\left({\omega }_{N-1}\right)\end{array}\right)$21

where the time series index (l) is sorted in ascending order of scaled periodogram value for each given frequency. Thus, the elements in a row of P_sorted are scaled periodograms for the same time series at N different frequencies, and the elements in a column are scaled periodogram values for different time series at the same frequency.

Let ${\left(F_{1,n},\hspace{0.17em}\dots \hspace{0.17em},F_{L,n}\right)}^T$ be the empirical CDF vector evaluated at $\left(p_{\left(1\right)}\left(w_n\right)\right),\hspace{0.17em}\dots \hspace{0.17em},$ ${\left(p_{\left(L\right)}\left(w_n\right)\right)}^T$ for frequency ω_n. The empirical CDFs at all N frequencies are arranged in the following matrix:

$F\triangleq \left(\begin{array}{ccc}{F}_{1,0}& & F_{1,N-1}\\ ⋮& ⋮& ⋮\\ F_{L,0}& & F_{L,N-1}\end{array}\right)$ 22

We find the FOGMP model parameters T_max and ${\sigma }_{\max }^2$ such that the model CDF $F_{1,n}(T_{\max },{\sigma }_{\max }^2)$ satisfies the following inequality:

$\begin{array}{cccc}{F}_{l,n}\left(T_{\max },{\sigma }_{\max }^2\right)\le F_{l,n}& \text{for}\hspace{0.17em}l=1,\hspace{0.17em}\dots \hspace{0.17em},L& \text{and}& n-0,\hspace{0.17em}\dots \hspace{0.17em},N-1\end{array}$ 23

The FOGMP model’s CDF is evaluated at the same ω_n and quantile values as the empirical CDF.

Figure 8(a) shows 2000 gray lagged product curves for 100-s-long time series of a simulated FOGMP with T = 50 s and σ = 1. The left-hand-side chart displays $F_{\mathrm{l},n}(T_{\max },\hspace{0.17em}{\sigma }_{\max })$ with blue dashed lines bounding four example sample quantile curves (solid black) at 2%, 16%, 84%, and 98% to determine T_max = 50 s and σ_max = 1.3. Figure 8(b) shows q(_τ)-CDF upper and lower bounds at example circular frequencies of ${\omega }_n={10}^{-2}\pi \hspace{0.17em}\text{rad}$, ${\omega }_n={10}^{-1}\pi \hspace{0.17em}\text{rad}$, and ${\omega }_n=10\pi \hspace{0.17em}\text{rad}$.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 8

Simulated data with a CDF lower-bounding FOGMP model for (a) example quantiles at all circular frequencies and (b) example circular frequencies at all quantiles

Figure 9 shows the differences between the sample CDF and the CDF- lower-bounding model CDF for all circular frequencies and quantiles. The color code represents the tightness of the bound, with a tight bound indicated in blue and a loose bound indicated in yellow. If the model failed to bound the empirical CDF, i.e., if Equation (23) was not satisfied, the CDF difference would be shown in red. In this example, the CDF bounds are loosest at the core of the distributions and for ${\omega }_n\ge 0.03$ rad, i.e., near and above the FOGMP corner frequency $\Delta t/T\pi =0.06\hspace{0.17em}\text{rad}$. This trend, which was also observed for the time-domain method, is consistent with snapshot overbounding where low-order parametric models of fat-tailed distributions are achievable at the cost of loose bounds of the distribution cores. Loose measurement error bounds can cause loose estimation error bounds, but they are needed to achieve high-integrity navigation using practical estimators such as KFs.

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 9

Contour plot showing the tightness of the CDF lower-bounding FOGMP model over the sample periodogram distribution

Positive values throughout indicate that the model is bounding at all frequencies and quantiles.

4.1 Time-Correlation Modeling Using a Sparse Data Set

We can first assess the impact of data sparsity on bound tightness in the measurement domain using the example L = 2000 simulated FOGMPs of N = 100 samples each. Based on this same finite data set and using the methods in Sections 2 and 3, we achieve high-integrity FOGMP models with T_min = 10 s, σ_min = 0.96, T_max = 90 s, and σ_max = 1.05 in the time domain and with T_max = 50 s and σ_max = 1.3 in the frequency domain. If we increase L, the time constants and variance values would approach the simulated FOGMP parameters of T = 50 s and ${\sigma }_v^2=1$.

For comparison with previous ACF bounding methods reported by Langel (2014), we reprocess the data by considering groups of 25 time series out of the original L = 2000, and we average their lagged products to generate 80 ACFs. The ACF upper bound has the following parameter values: T_max = 100 s, σ_max = 1.414, giving a looser bound than that obtained with the new lagged product method (T_max = 90 s, σ_max = 1.05); the ACF lower bounding parameters are undefined because sample ACFs take negative values. A major benefit of the new method is that one can provide T_min and σ_min values for any data set over any range of operational lag times $\tau \in \left(0,{\tau }_{op}\right)$. In parallel, for comparison with the PSD upper-bounding approach of Gallon et al. (2020) and Langel et al. (2020), we generate 80 PSDs by averaging five scaled periodograms over 500-s-long time series. The time series lengths had to be extended as compared to N = 100 to limit the looseness of the PSD bound. This data partitioning favors the baseline PSD approach, which would otherwise produce an unrealistically loose bound. The rationale for this process, which requires some background information on PSD estimation, is described in Appendix E. The resulting PSD bound gives T_max = 50 s and σ_max = 2.5, which is again looser than that obtained with the new periodogram-based approach (T_max = 50 s, σ_max = 1.3).

4.2 Time-Correlation Modeling of Non-FOGMP Data

The methods described in this paper are applicable not only to samples drawn from a single FOGMP, but also to data with an unknown time-correlation structure. To analyze this more general case, in this section, we simulate non-FOGMP data consisting of L = 2000 time series with a length of 100 s and a sampling interval Δt = 1 s, i.e., N = 100. The data set is built as a composite of two distinct unit-variance FOGMPs including 1000 sample time series with time constant T_A = 50 s and 1000 time series with T_B = 15 s. Thus, this data set is a mixture of two unit-variance random processes v_A,n and v_B,n with respective PDFs p(v_A,n) and p(v_B,n), which can be expressed as follows:

$v_n\sim \frac{1}{2}p\left(v_{A,n}\right)+\frac{1}{2}p\left(v_{B,n}\right)$24

where:

$\begin{array}{ccc}{v}_{A,n}=\text{exp}\left(-\frac{\Delta t}{T_A}\right)v_{A,n-1}+{\eta }_{A,n-1},& \text{and}& v_{B,n}=\text{exp}\left(-\frac{\Delta t}{T_B}\right)v_{B,n-1}+{\eta }_{B,n-1}\end{array}$25

with:

$\begin{array}{ccc}{\eta }_{A,n}\sim N\left(0,1-\text{exp}\left(-\frac{2\Delta t}{T_A}\right)\right)& \text{and}& {\eta }_{B,n}\sim N\left(0,1-\text{exp}\left(-\frac{2\Delta t}{T_B}\right)\right)\end{array}$26

The statistics of this sparse non-FOGMP data set are modeled using FOGMPs in the time and frequency domains.

Using the time-domain procedure in Section 2.3, we found upper- and lower-bounding FOGMP models defined by T_max = 85 s, σ_max = 1.1 and T_min = 10 s, σ_max = 1, respectively. Using the frequency-domain procedure in Section 2.3, we found a PSD upper-bounding FOMGP model with T_max = 40 s and σ_max = 1.5. The tightness of the time-domain and frequency-domain bounding models is shown in Appendix C.

5.1 TIME-DOMAIN METHOD

The normalized IF CMC data for all 32 satellites are sampled at regular intervals of Δt = 1 s and partitioned into L = 904 time series of N = 600 samples each. The data are first processed via the time-domain method to derive a model based on a pair of FOGMPs. The differences between sample data and bounding model CDFs are shown in Figure 12 for lag times of τ = 1 s to τ = 600 s and over all quantiles. The FOGMP time-correlation model parameters are T_max = 1500 s, σ_max = 1.3 and T_min = 100 s, σ_min = 0.975. For comparison, the ACF upper-bounding approach reported by Langel (2014) would give parameter values of T_max = 1800 s, σ_max = 1.9; the ACF lower bounding parameters over τ_op = 10 min are undefined because the ACFs take negative values (shown in Figure 20(a) in Appendix E).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 12

Time-domain model: tightness of the (a) lower and (b) upper bounding FOGMP model CDF over the experimental data’s lagged product CDF

5.2 Frequency-Domain Method

The normalized IF CMC data are processed via the frequency-domain method with the same partitioning as in Section 5.1. The scaled periodograms are evaluated at circular frequencies ${\omega }_n\in \left(\pi \times {10}^{-4},\pi \right)$. The upper frequency limit is determined by the Nyquist frequency of 0.5 Hz for the 1-s sampling interval, which corresponds to a circular frequency of π rad. The lower limit should be significantly lower than $2\pi \Delta t(N\Delta t)$ to capture low-frequency content variations. For the 600-s-long time series data partition, the frequency limit should be lower than 10^-3 rad, where the periodograms start flattening. The parameters of the bounding FOGMP time-correlation model are T_max = 180 s and σ_max = 1.35. Figure 13 presents the differences between the sample IF CMC and the model FOGMP scaled periodogram CDF surfaces. For comparison, using the same partition, the PSD upper-bounding approach reported by Gallon et al. (2020) and Langel et al. (2020) yields T_max = 270 s and σ_max = 2.2, corresponding to a looser bounding model than that obtained by using the new periodogram overbounding approach (shown in Figure 20(b) in Appendix E).

• Download figure
• Open in new tab
• Download PowerPoint

FIGURE 13

Frequency-domain model: tightness of the CDF lower-bounding FOGMP model over the experimental data’s periodogram distribution