High-Integrity Modeling of Nonstationary Noise Processes for GNSS/INS Integration

Elisa Gallon; Mathieu Joerger,; Boris Pervan

doi:10.33012/navi.729

Abstract

This paper describes a power spectral density (PSD) bounding method for deriving high-integrity models of stationary and nonstationary time-correlated measurement error processes. The method is intended for safety-critical commercial aircraft navigation where robust sensor models are required to predict error bounds on position and orientation estimates. These bounds are used both in navigation system design for integrity performance analyses and in operation to determine whether a pilot should proceed with an operation. In prior work, we used PSD upper-bounding to obtain high-integrity models of time-correlated global navigation satellite system (GNSS) measurement errors. However, the method was limited to stationary processes. In this paper, we derive an approach to expand the concept of PSD bounding to nonstationary error modeling for Kalman filter-based estimation using GNSSs and inertial navigation systems in aircraft navigation applications.

Keywords

1 INTRODUCTION

In safety-critical aviation applications where flight path deviations can quickly become hazardous to passengers and crew, navigation integrity is of the utmost importance. Integrity is a measure of trust in the accuracy of the information supplied by the navigation system (International Civil Aviation Organization (ICAO) 2001). The integrity risk is the probability that the positioning error exceeds a predefined limit of acceptability, or alert limit, without a timely warning. For aircraft navigation, integrity risk requirements are extremely stringent, equal to or lower than 10^-7. The Global Positioning System (GPS) alone does not meet such requirements with high availability and continuity (Pullen & Joerger, 2020). Therefore, augmentation systems were developed for aviation applications, including ground-based augmentation systems (GBASs) (Pervan, 2020) and space-based augmentation systems (SBASs) (Walter, 2020).

Navigation integrity is ensured by predicting an upper bound on the probability that an undetected state-estimate error exceeds the alert limit. At the source of the estimation error is the uncertainty of the sensor measurements. Therefore, high-integrity navigation demands careful prior evaluation and modeling of measurement errors.

Overbounding theory aims to find models of empirical measurement error distributions that guarantee upper bounds on estimation errors. These bounds are predictable and are proven to be conservative in the sense that they do not underestimate actual errors of linear (and linearized) estimators. Overbounding is applied in SBASs, GBASs, and advanced receiver autonomous integrity monitoring for “snapshot” positioning, i.e., positioning using measurements received at one time instant (Advanced Receiver Autonomous Integrity Monitoring (ARAIM) Working Group C (WGC) 3.0, 2016; Blanch et al., 2019).

Algorithms that filter measurements over time have been investigated for high-integrity air and ground navigation applications (Joerger & Pervan, 2020). Filtering could provide a certifiable means for guiding aircraft closer to the runway using satellite-based ranging measurements. However, significant challenges remain in deriving measurement error models that guarantee integrity risk bounds for recursive algorithms when the measurement error time correlation is unknown.

In this paper, we introduce a methodology for high-integrity modeling of time-correlated stochastic errors when the time correlation is uncertain and not stationary.

This methodology is intended for a Kalman filter (KF)-based implementation combining a global navigation satellite system (GNSS) with an inertial navigation system (INS) for use in aircraft navigation applications. GNSS measurements are impacted by time-correlated atmospheric delays, satellite orbit and clock errors, and multipath. INS data are impacted by time-correlated accelerometer and gyroscope errors. Each source of error must be carefully modeled and bounded. For KF implementations, this not only implies that instantaneous error distributions be addressed using overbounding, but it also requires that the error time correlation be accounted for.

Section 2 of this paper provides an overview of the state of the art of high-integrity error modeling and highlights current limitations for KF applications— in particular, the restriction of existing power spectral density (PSD) bounding methods to stationary processes. To extend the scope of PSD bounding methods, it is necessary to first divide “nonstationary” processes into subcategories that can be individually analyzed. Section 3 defines criteria for classifying error processes with respect to their unique characteristics. We identify three categories of nonstationarity and explain their attributes. Sections 4 and 5 then develop corresponding robust modeling approaches. For each category of error, a bounding methodology is derived. Experimental examples are then given to demonstrate the proposed methodology using real GNSS and INS data. Section 6 provides concluding remarks. Elements of this research contributed to a Ph.D. dissertation by Gallon (2023).

2 OVERBOUNDING, AUTOCORRELATION FUNCTION BOUNDING, AND PSD BOUNDING

Overbounding theory provides tools to provably bound the impacts of Gaussian and non-Gaussian measurement errors on estimation uncertainty. Prior research on overbounding for high-integrity error modeling has yielded significant results that are currently implemented in commercial aviation navigation safety. The concept of cumulative distribution function (CDF) overbounding was introduced by DeCleene (2000). The original CDF overbounding method models an empirical, non-Gaussian error distribution with an overbounding model distribution. This method is limited to ideal cases in which the actual error distribution is zero-mean (or non-zero-mean with a known mean value), symmetric, and unimodal. Realistic error distributions do not always satisfy these conditions, which motivated the development of paired-CDF overbounding (Rife et al., 2006). Paired overbounding uses two bounding functions to strictly and conservatively overbound empirical error distributions.

Gaussian functions are often employed as overbounding models for sensor errors because their combination through linear estimators results in a known and predictable Gaussian estimation error distribution. Blanch et al. (2019) combined single-CDF, paired-CDF, and excess-mass overbounding methods to achieve tighter bounds than paired-CDF overbounding.

Although such error models can be employed in snapshot positioning algorithms, they are insufficient for sequential implementations because they do not address the stochastic dynamics of the errors over time.

To address this limitation, Langel et al. (2014) derived an analytical bound on sequential linear estimation errors using auto-correlation function (ACF) bounding. Unfortunately, the ACF-based approach requires cumulative storage of all state propagation and observation matrix coefficients over time following filter initiation. This approach is unsuitable for KF implementations over long filtering durations, especially when high-rate sensors such as INSs are used.

More recently, an approach to upper-bound the KF estimate error variance in the presence of uncertain time-correlated noise was developed by Langel et al. (2020). This approach relies on upper-bounding a modified PSD derived from the Fourier transform of a windowed ACF for stationary errors.

The PSD bounding method, unlike ACF bounding, is not restricted to fixed-interval implementations and is compatible with KF approaches. However, the treatment of experimental data for the modeling of wide-sense stationary (WSS) processes has yet to be formalized, and the case of nonstationary error processes remains unaddressed. For instance, inertial measurements have nonstationary noise components, including flicker noise, that cannot be directly modeled using the approach of Langel et al. (2020) because the PSD cannot be obtained from the ACF. Therefore, new modeling methods must be derived. This paper builds upon the concept of PSD bounding and extends it to nonstationary errors.

As shown in the following section, the working definitions of ACFs and PSDs vary depending on whether the process under consideration is stationary. Therefore, nonstationary characteristics must be defined and categorized. Section 3 defines the stationarity criteria used to classify error processes.

3 ERROR CLASSIFICATION

This section aims to classify stationary and nonstationary processes (Gallon, 2023), serving as an outline for the remainder of the paper, where modeling methods are derived for each process category.

3.1 Stationary Error Processes

Mean values and ACFs are key descriptors of the statistical properties of a random process. The mean value of a random process x at time t can be expressed as follows:

μ_{x} (t) = E [x (t)]

where E[] is the expected value operator. In addition, the ACF, i.e., the correlation between the values of the random process x at two different times $t - τ / 2$ and $t + τ / 2$ is defined as follows:

R_{x x} (t, τ) = E [x (t + \frac{τ}{2}) x (t - \frac{τ}{2})]

If $μ_{x} (t)$ and $R_{x x} (t, τ)$ change over time t, then process x is said to be nonstationary. If $μ_{x} (t)$ and $R_{x x} (t, τ)$ do not vary with t (i.e., $μ_{x} (t) = μ_{x}$ and $R_{x x} (t, τ) = R_{x x} (τ)$ )’ then process x is said to be WSS. Note that we use different notations for a stationary ACF R_xx and a non-stationary (instantaneous) ACF R_xx. The following subsection aims to further refine the classification of nonstationary errors.

3.2 Nonstationary Error Processes

When process x is nonstationary, the ACF in Equation (2) becomes an instantaneous ACF (Bendat & Piersol, 2010), depending on both t and the lag τ. The ACF can be determined in one of two ways:

(a) analytically, if the process is sufficiently well understood to evaluate Equation (2), or
(b) empirically, from previously collected data.

In practice, option (a) is seldom possible. Option (b) is feasible only when variations in the ACF with respect to time t are much slower than variations with respect to τ, which can be expressed as follows:

| \frac{\partial R_{x x}}{\partial t} | ≪ | \frac{\partial R_{x x}}{\partial τ} |

Therefore, we distinguish two relevant categories of nonstationary error processes (Gallon, 2023):

Slowly varying ACFs: Equation (3) holds, permitting computation of stationary sample ACFs.
Time-invariant PSDs: Equation (3) does not hold, yet the PSD is independent of t; thus, the sample PSDs remain valid.

A third case—where Equation (3) fails and the PSD is time-varying—is illustrated in Appendix A by a rapid tropospheric-delay change during a severe storm. Such transients are best treated as fault conditions, handled via real-time detection and measurement exclusion rather than stochastic modeling, and thus lie outside the scope of this work.

The classification of error processes is summarized in Figure 1.

FIGURE 1

Error classification tree outlining the sections of this paper

4 FREQUENCY-DOMAIN BOUNDING OF STATIONARY INPUT ERROR PROCESSES

Langel et al. (2020) and Langel et al. (2024) derived a frequency-domain criterion for modeling stationary correlated noise through PSD bounding. This section gives a simpler, alternative proof of this criterion for time-invariant and time-varying KFs. KF-based INS/GNSS processing may be impacted by linearization errors, which are assumed to be small in this research. Bounding linearization errors is outside the scope of this paper.

Because the KF is a linear filter, we can separately account for the impacts of independent additive error sources and then sum up their individual contributions. This approach is valid both for multiple sensors contributing to a state’s estimation (e.g., multiple independent GNSS satellite measurements) and for multiple sources of uncertainty impacting the same measurement (e.g., a GNSS ranging signal’s atmospheric delays, satellite orbit and clock ephemeris errors, and multipath errors). This approach is also valid for any linear combination of KF outputs. Therefore, for clarity of exposition in the following derivation, we consider a scalar KF fed by zero-mean sensor and process noise contributions.

4.1 Time-Invariant KF with WSS Input

We first consider the simplest case of a steady-state scalar KF fed by zero-mean stationary noise. Consider a scalar error input x(t)¹ with PSD S_xx(f). In general, the KF is designed using an a priori measurement error PSD model ${\bar{S}}_{x x} (f)$ such that ${\bar{S}}_{x x} (f) \neq S_{x x} (f)$ because S_xx(f) is uncertain. Let $\bar{H} (f)$ be the transfer function of the KF from x(t) to some scalar output of interest y(t).

The mean of the KF output error is zero because the input is zero-mean, and the actual, unknown output error variance is given by the following:

σ_{y}^{2} = \int_{- \infty}^{\infty} {| \bar{H} (f) |}^{2} S_{x x} (f) d f

The KF produces a known output error variance that can be expressed as follows:

{\bar{σ}}_{y}^{2} = \int_{- \infty}^{\infty} {| \bar{H} (f) |}^{2} {\bar{S}}_{x x} (f) d f

Clearly, if the following inequality is satisfied

\begin{matrix} {\bar{S}}_{x x} (f) > S_{x x} (f) & \forall f \in ℝ \end{matrix}

then we achieve ${\bar{σ}}_{y}^{2} > σ_{y}^{2}$ In other words, the PSD upper-bound ${\bar{S}}_{x x} (f)$ defined in Equation (6) can be used to find a time-correlated measurement error model that guarantees an upper bound on the KF estimation error variance (Gallon et al., 2021b; Langel et al., 2020). In the case of a KF with multiple output variables, the variance bound ${\bar{σ}}_{y}^{2}$ can be derived for any linear combination y(t) of the output states, as proven in Appendix B. In state space, this variance bound lays on an ellipsoid entirely enclosing the actual state-estimate covariance ellipsoid (the actual KF output covariance matrix is bounded in the positive-definite sense). Integrity risk bounds are derived from such variance bounds (Advanced Receiver Autonomous Integrity Monitoring (ARAIM) Working Group C (WGC) 3.0, 2016).

4.2 Time-Varying KF with WSS Input

Let us consider a time-varying filter with a transient response $\bar{h} (τ, t)$ to a unit impulse $δ (τ)$ . The filter’s output for an input x(t) is given by the following convolution integral:

y (t) = \int_{- \infty}^{\infty} \bar{h} (α, t) x (t - α) d α

The variance of the output can then be derived as follows:

\begin{array}{l} σ_{y}^{2} (t) & = & E [y {(t)}^{2}] \\ = & E {\int_{- \infty}^{\infty} \bar{h} (α, t) x (t - α) d α \int_{- \infty}^{\infty} \bar{h} (β, t) x (t - β) d β} \\ = & \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \bar{h} (α, t) \bar{h} (β, t) E [x (t - α) x (t - β)] d α d β \\ = & \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \bar{h} (α, t) \bar{h} (β, t) R_{x x} (α - β) d α d β \end{array}

Using the Wiener-Khinchin theorem, we can express the ACF R_xx in terms of the PSD S_xx as follows:

R_{x x} (α - β) = \int_{- \infty}^{\infty} S_{x x} (f) e^{- i 2 π f (α - β)} d f

We can then rewrite Equation (8) as follows:

\begin{array}{l} σ_{y}^{2} (t) & = & \int_{- \infty}^{\infty} S_{x x} (f) \int_{- \infty}^{\infty} \bar{h} (α, t) e^{- i 2 π f α} d α \times \int_{- \infty}^{\infty} \bar{h} (β, t) e^{i 2 π f β} d β d f \\ = & \int_{- \infty}^{\infty} S_{x x} (f) \bar{H} (f, t) \bar{H} {(f, t)}^{*} d f \\ = & \int_{- \infty}^{\infty} S_{x x} (f) | \bar{H} (f, t) |^{2} d f \end{array}

Similar to Equation (5) for the time-invariant case, Equation (10) shows that for a time-varying KF if ${\bar{S}}_{x x} (f) > S_{x x} (f), \forall f \in ℝ$ , then ${\bar{σ}}_{y}^{2} (t) > σ_{y}^{2} (t), \forall t \in ℝ^{+}$ .

4.3 Experimental Procedure

The PSD bounding methodology was first introduced by Langel et al. (2014, 2020). The finite nature of experimental data can be modeled as a rectangular window applied over a random process. Rectangular windows, with their sharp edges, introduce spectral leakage during the conversion from the time domain to the frequency domain and degrade the quality of the resulting PSD estimate, thereby loosening the PSD upper bound that determines the error model. Any model whose PSD upper-bounds the PSD of the error can be a good candidate with regard to integrity. However, with respect to accuracy or continuity, we want a tight PSD upper bound to achieve as small a bounding estimation error variance as possible.

To counter this effect, tapering windows are often used (e.g., Hamming window, Hann window). Additionally, it is important to note that the longer the data set, the more accurate the ACF estimate will be. However, for most applications, the mission duration, T₁, is finite. In practice, the estimator does not use observations over lag times longer than T₁ (e.g., in aviation, 18 h is the longest commercial flight duration on record). Therefore, the information contained in the ACF for lag times longer than T₁ is not needed, but it can be leveraged in the windowing process to reduce spectral leakage (Langel et al., 2014), tighten the PSD upper bound, and ultimately lower the final bounding value of the estimate error variance (Gallon, 2023).

The proposed methodology (Langel et al., 2014, 2020) relies on a four-step process:

Estimate the ACF R_xx of the input process x.
Multiply R_xx by a tapering window Λ.
Take the Fourier transform of the windowed ACF to obtain the PSD of x.
Derive a high-integrity, time-correlated error model by tightly upper-bounding the PSD (Gallon, 2023).

This process is represented in Figure 2.

FIGURE 2

PSD bounding procedure for WSS input

The window used in this work was developed by Langel et al. (2020) and is determined over lag time τ as follows:

Λ (τ, T_{1}, T_{2}) = {\begin{matrix} 0, & for τ < - T_{2} \\ \frac{- 1}{e^{4 η / (1 - η^{2})} + 1} + 1, & for - T_{2} \leq τ < - T_{1} \\ 1, & for - T_{1} \leq τ \leq T_{1} \\ \frac{1}{e^{4 η / (1 - η^{2})} + 1}, & for T_{1} < τ \leq T_{2} \\ 0, & for τ > T_{2} \end{matrix}

where:

η = {\begin{matrix} 1 - 2 \frac{T_{1} + τ}{T_{1} - T_{2}} & for - T_{2} \leq τ < - T_{1} \\ 1 + 2 \frac{T_{2} - τ}{T_{1} - T_{2}} & for T_{1} < τ < T_{2} \end{matrix}

−T₂ and +T₂, respectively, are the lag times before and after which the window reaches a value of zero.

With this window, ACF values for lag times smaller than T_l remain unchanged, whereas ACF values for lag times larger than +T₂ or smaller than −T₂ are set to zero. When T_l and T₂ are further apart, less spectral leakage is found in the PSD estimate. An example of such a window is shown in Figure 2 for T₁ = 7 h and T₂ = 14 h.

The sample PSD estimate of process x can then be estimated by taking the Fourier transform of the windowed ACF. This step is expressed as follows:

\begin{matrix} S_{x} (f) = \int_{- T_{2}}^{T_{2}} Λ (τ, T_{1}, T_{2}) R_{x x} (τ) e^{- j 2 π f τ} d τ & \forall f \in ℝ \end{matrix}

5 FREQUENCY-DOMAIN BOUNDING OF NONSTATIONARY INPUT ERROR PROCESSES

In Section 3, we defined a condition in Equation (3) on the instantaneous ACF $R_{x x} (τ, t)$ for a nonstationary process x(t) at time t and lag time τ. This condition is used to categorize “slowly” versus “rapidly” changing nonstationarity. In this section, we present methodologies for error modeling in each of these two cases, which we then illustrate using experimental data (Gallon, 2023).

5.1 Slowly Varying ACF

5.1.1 Theoretical Development

Given an instantaneous ACF $R_{x x} (τ, t)$ , we can define the instantaneous PSD as follows (Bendat & Piersol, 2010):

S_{x x} (f, t) = \int_{- \infty}^{\infty} R_{x x} (τ, t) e^{- j 2 π f τ} d τ

Conversely, we have the following:

R_{x x} (τ, t) = \int_{- \infty}^{\infty} S_{x x} (f, t) e^{j 2 π f τ} d f

The KF is typically designed using an a priori parametric model ${\bar{S}}_{x x} (f, t)$ such that ${\bar{S}}_{x x} (f, t) \neq S_{x x} (f, t)$ .

Following a development similar to that of the time-varying KF with a stationary input process in Section 4.2, we can express the output variance as follows:

σ_{y}^{2} (t) = \int_{- \infty}^{\infty} \int_{- \infty}^{\infty} \bar{h} (α, t) \bar{h} (β, t) R_{x x} (α - β, t) d α d β

which leads to the following result:

σ_{y}^{2} (t) = \int_{- \infty}^{\infty} {| \bar{H} (f, t) |}^{2} S_{x x} (f, t) d f

In this case again, if $\begin{matrix} {\bar{S}}_{x x} (f, t) > S_{x x} (f, t) & \forall f \in ℝ \end{matrix}$ and $\forall t \in ℝ^{+}$ , then $\begin{matrix} {\bar{σ}}_{y}^{2} (t) > σ_{y}^{2} (t) & \forall t \in ℝ^{+} \end{matrix}$ .

Regardless of whether ${\bar{H}}_{x x} (f, t)$ is time-invariant, it is desirable for the input error model PSD ${\bar{S}}_{x x} (f, t)$ to be time-independent because it is difficult, if not impossible, to ensure the instantaneous validity of error models obtained from experimental data collected prior to the time at which the KF was actually implemented.

Instead, we follow a more robust and practical approach: we collect a long, non-stationary error data set that captures the time variability of S_xx(f, t). We then define a time-invariant upper-bounding PSD as follows:

\begin{matrix} {\bar{S}}_{x x} (f) > S_{x x} (f, t), & \forall t \in ℝ^{+} & and & \forall f \in ℝ \end{matrix}

We obtain this time-invariant upper-bounding PSD by dividing a long continuous data set into smaller quasi-stationary subsets that satisfy Equation (3), with center times t₁, t₂,…, t_N, to compute S_xx (f, t₁), S_xx (f, t₂),..., S_xx (f, t_N). We then choose ${\bar{S}}_{x x} (f)$ to satisfy the following inequality:

\begin{matrix} {\bar{S}}_{x x} (f) \geq \max_{1 \leq k \leq N} S_{x x} (f, t_{k}), & \forall f \in ℝ \end{matrix}

Based on our previous work (Gallon, 2023; Gallon et al., 2021b; Jada & Joerger, 2021; Joerger et al., 2023), we have found that the observation period required to obtain an accurate sample PSD estimate, S_xx (f, t_k), should be at least ten times the correlation time constant. We also assume that the sample PSDs already capture the worst-case measurement-error conditions expected in practice. The resulting error model can guarantee bounding only for these observed conditions; thus, the data collection plan must ensure that every relevant operational scenario is represented.

5.1.2 Experimental Results

In our prior work (Gallon et al., 2021b), high-integrity, time-correlated error models were derived for tropospheric errors by using the PSD bounding method in Equation (19). GNSS tropospheric delays were generated at 100 worldwide ground stations over the year 2018. Tropospheric residuals were evaluated after correction from two different models named MOPS and GPTw2 (as described by Gallon et al. (2021b)). Nonstationarities were regularly observed throughout the tropospheric error analysis.

For example, the GPT2w residuals observed at the GPS receiver and meteorological station “BADG” in Russia are shown in Figure 3. This time profile over 2018 shows a nonstationary process whose variance changes over time, with a clear increase between days of the year 117 and 270. Following the procedure described in Equation (19), this data set was split into three quasi-stationary subsets delimited by the red dashed lines shown in the figure.

FIGURE 3

Practical example of slowly varying nonstationarity in the tropospheric residuals of a station in Russia in 2018

This approach was performed for tropospheric residuals at all 100 stations. The process was partly automated using mathematical stationarity criteria detailed by Gallon et al. (2021b). The stationary subsets were then analyzed in the frequency domain, as shown in Figure 4, which contains all of the sample PSDs superimposed. The PSDs were then upper-bounded via a first-order Gauss-Markov random process (FOGMRP) model to find the upper-bounding PSD model ${\bar{S}}_{x x} (f)$ , as displayed by a red curve in Figure 4. We selected a FOGMRP because it is a two-parameter model (in this case, a time constant of 20 h and standard deviation of 9 cm) that can be easily incorporated into a KF by state augmentation.

FIGURE 4

Experimental results: Upper bounding of tropospheric residual PSDs of 100 stations worldwide in 2018

5.2 Rapidly Varying ACF with a Time-Invariant PSD

This section addresses the case in which the inequality in Equation (3) is not satisfied but the PSD is not a function of time. This category relates to processes that are nonstationary because they are unstable - e.g., 1/fⁿ noise (Brown & Hwang, 2012).

In these cases, it may be tempting to directly use Equations (17) and (19). However, these results were derived assuming that an instantaneous ACF was obtainable, i.e., assuming that Equation (3) was satisfied.

In the example of 1/f noise, the process is nonstationary, with an ACF that depends on both time t and lag time t. When t = τ, the ACF values go to infinity (Brown & Hwang, 2012), and the error variance becomes infinitely large. Therefore, the PSD estimation methodology that was introduced above, which relied on taking the Fourier transform of the ACF, cannot be applied.

Instead, we evaluate the PSD using periodograms rather than ACFs. We also assume that the KF error output process is stable (even if the input is not), which must be true for a useful KF.

5.2.1 Theoretical Development

First, we express the PSD of the output y(t) of a time-invariant KF as follows:

\begin{array}{ll} S_{y y} (f) & = \lim_{T \to \infty} \frac{1}{T} E \{Y (f) Y {(f)}^{*}\} \\ = \lim_{T \to \infty} \frac{1}{T} E \{[\bar{H} (f) X (f)] {[\bar{H} (f) X (f)]}^{*}\} \\ = {|\bar{H} (f)|}^{2} \lim_{T \to \infty} \frac{1}{T} E \{X (f) X {(f)}^{*}\} \\ = {|\bar{H} (f)|}^{2} S_{x x} (f) \end{array}

where X(f) and Y(f) are the Fourier transforms of x(t) and y(t), respectively. As long as the output error process y(t) is stable, the following integrals converge:

σ_{y}^{2} = \int_{- \infty}^{\infty} S_{y y} (f) d f = \int_{- \infty}^{\infty} {| \bar{H} (f) |}^{2} S_{x x} (f) d f

This form is the same as in Equation (5), and the same conclusion on variance bounding using a PSD bound ${\bar{S}}_{x x} (f)$ applies.

Next, we consider the PSD of the output y(t) of a time-varying KF, which can be written as follows:

\begin{array}{ll} S_{y y} (f, t) & = \lim_{T \to \infty} \frac{1}{T} E {Y (f, t) {Y (f, t)}^{*}} \\ = \lim_{T \to \infty} \frac{1}{T} E {[\bar{H} (f, t) X (f)] {[\bar{H} (f, t) X (f)]}^{*}} \\ = {| \bar{H} (f, t) |}^{2} \lim_{T \to \infty} \frac{1}{T} E {X (f) {X (f)}^{*}} \\ S_{y y} (f, t) & = {| \bar{H} (f, t) |}^{2} S_{x x} (f) \end{array}

Again assuming a stable output error process y(t), we can express the variance at time t as follows:

σ_{y}^{2} (t) = \int_{- \infty}^{\infty} S_{y y} (f, t) d f = \int_{- \infty}^{\infty} {| \bar{H} (f, t) |}^{2} S_{x x} (f) d f

This equation is of the same form as Equation (10), and the same conclusion on $σ_{y}^{2} (t)$ bounding using ${\bar{S}}_{x x} (f)$ applies.

5.2.2 Experimental Results

In this section, we apply the periodogram-based PSD estimation to the example INS accelerometer errors, which are impacted by “flicker” (1/f) noise, among other sources. (Gyroscope errors can be similarly treated.)

We use the following simple example model of accelerometer error:

a_{m} = a_{t} + b (t) + p (t) + v_{s} (t)

where a_m is the measured accelerometer output. The true value of the measured variable a_t is impacted by a time-varying “bias” b(t), an acceleration random walk (RW) p(t), and a velocity RW v_s(t).

The time-varying bias can be expressed as follows:

b (t) = b_{r} + b_{s} (t)

The initial bias b_r, known as bias repeatability (or turn-on bias stability), differs at every power-up because of signal processing initialization and physical properties (thermal, mechanical, and electrical variations). The bias instability b_s(t) is the time-varying component capturing the drift from the starting value. The bias instability can be modeled as flicker noise.

The error components modeled in Equations (24) and (25) are present in data collected over several 10-h-long periods in static conditions at a 125-Hz sampling frequency using a STIM-300 inertial measurement unit (IMU) from Sensonor (Gallon et al., 2021a; Sensonor, n.d.). The STIM-300 is a small, tactical-grade, low-weight, high-performance IMU that contains three highly accurate micro-electromechanical system gyroscopes, three high-stability accelerometers, and three inclinometers. Additionally, a temperature sensor was used to monitor environmental temperature variations. Five data sets meeting the IMU’s constant temperature requirements were selected. The PSDs derived from these sets using periodograms are represented in gray in Figure 5.

FIGURE 5

Experimental results: Upper-bounding sample accelerometer error PSDs of the STIM-300 IMU

To upper-bound these errors in the PSD domain, we define a bounding model ${\bar{S}}_{W G R}$ consisting of a white noise (WN) component ${\bar{S}}_{W N}$ , a FOGRMP component ${\bar{S}}_{G}$ , and an RW component ${\bar{S}}_{R W}$ :

{\bar{S}}_{W G R} (f) = {\bar{S}}_{W N} (f) + {\bar{S}}_{G} (f) + {\bar{S}}_{R W} (f)

The velocity RW v_s is modeled as white noise with PSD S_WN (defined by the parameter N₀), the acceleration RW p is modeled as a RW with PSD S_RW (defined by the parameter K), and the time-varying bias b is modeled as a FOGMRP with PSD S_G (defined by the parameters σ_G and τ_G ). These parameters are listed in Table 1. The time-varying bias b is actually flicker noise, but we bound its error contribution using a FOGMRP because flicker noise cannot be directly modeled in the state-space domain (Brown & Hwang, 2012).

View this table:

TABLE 1

Inertial Error Model Components

In Figure 5, the black curve represents the accelerometer specifications as provided by the manufacturer. Although the black curve appears to accurately model the average case, it does not upper-bound the sample PSDs. Instead, we can derive a PSD bounding model by selecting a set of values for the parameters in Table 1, i.e., K, σ_G, τ_G, and N₀, such that ${\bar{S}}_{W G R} (f)$ . This approach results in the red curve in Figure 5, which upper-bounds the sample PSD curves over the entire frequency range.

5.2.3 PSD Estimation Challenges for Unstable Errors

In practice, estimating the PSD S_xx (f) from experimental data to derive a model ${\bar{S}}_{x x} (f)$ is challenging. The first limitation stems from the range of frequencies achievable using observations, which is limited as follows:

Δ f \leq f \leq f_{N}

where f_N is the Nyquist frequency defined by the sampling interval At: $Δ t : f_{N} = 1 / (2 Δ t)$ . The second limitation is the frequency resolution Δf, which depends on the length of the data set, T_D: Δf = 1/T_D. Therefore, the estimate of S_xx(f) is determined by T_D and Δt.

In addition, PSD estimates derived from periodograms are impacted by spectral leakage, which can cause PSD bounds to be larger than necessary. Langel et al. (2014) and Gallon et al. (2021b) estimated PSDs using Fourier transforms of ACFs tapered with smooth-edged windows. Tapering the windows helps to limit spectral leakage (see Section 4.3) and tighten PSD bounds. Unfortunately, this approach cannot be applied to the case of unstable errors because PSDs derived from peri-odograms can be impacted by spectral leakage.

A loosely controlled data collection environment (e.g., an INS undergoing temperature variations or vibrations) can also contribute to degraded measurement models. For instance, INS manufacturers often provide specifications in the form of sensor error Allan variance (AV) curves rather than PSDs. The AV data are often of high quality, collected in highly controlled and stable environments using temperature and vibration isolation chambers. However, the parameters extracted from these data sets prioritize model accuracy over integrity, i.e., they are not designed to overbound measurement error distributions.

Therefore, for unstable errors, it would be beneficial for future work to complement PSD bounding with an AV-domain method to obtain tighter bounding error models.

6 CONCLUSION

Recent studies on high-integrity error time correlation modeling for aircraft navigation have assumed stationary error processes. This assumption is rarely applicable in practice.

This article classified error processes with respect to their stationarity properties, used PSD upper-bounding to model WSS processes, and addressed error bounding for two types of nonstationary processes. We derived a criterion to identify slowly varying nonstationary processes for which the error’s ACF decays quickly relative to the time required for the process itself to change. For such processes, we showed that the PSD upper-bounding approach used for WSS processes can be extended by splitting nonstationary time series into segments of stationary data. We then showed that for rapidly varying nonstationary processes whose PSDs are not a function of time, we can find high-integrity error models by bounding PSDs derived from periodograms.

Rapidly varying time-dependent processes such as large, rarely occurring, step-like errors may be treated as measurement faults requiring detection and exclusion and will be addressed in future research.

HOW TO CITE THIS ARTICLE:

Gallon, E., Joerger, M., & Pervan, B. (2026). High-Integrity modeling of nonstationary noise processes for GNSS/INS integration. NAVIGATION, 73. https://doi.org/10.33012/navi.729

ACKNOWLEDGMENT

The authors would like to thank the Federal Aviation Administration for their support of this research. However, the opinions in this paper are those of the authors and do not necessarily represent those of any other person or organization.

APPENDIX A | RAPIDLY VARYING NONSTATIONARITY AND TIME-VARYING PSD: TROPOSPHERIC DELAYS DURING STORMS

This appendix addresses the case in which the input process does not satisfy Equation (3) and its PSD is a function of time. This category of nonstationary error is the most challenging category because its behavior is highly unpredictable. Such cases can be observed, for example, in GNSS tropospheric errors during storm conditions.

Figure 6 shows tropospheric residuals at the GNSS receiver station “KOKV,” which is located in the southeast region of Hawaii. The residual generation procedure has been described by Gallon et al. (2021b). In August 2018, Hurricane Lane moved across the Hawaiian Islands, reaching Category 5 on August 22, which is day of year (DoY) 234. Hurricane Lane was the wettest tropical cyclone on record in Hawaii, with rainfall accumulations of up to 1,500 mm. It remained a Category 5 hurricane for 5 days before being downgraded to a tropical depression on August 28 (DoY 240) and dissipating the following day. The hurricane’s impact on tropospheric residuals can be clearly observed in Figure 6. On DoY 234, the first day of the hurricane, the residual mean shifted from -0.05 m to 0.05 m and remained near this value for the entire duration of the storm. When the storm ended on DoY 241, the residual mean returned to the initial value of -0.05 m. This type of nonstation-arity is rare. The impact of such events on estimation integrity would likely be best mitigated by outlier detection.

FIGURE 6

Practical example of rapidly varying nonstationarity in tropospheric residuals from 2018 at a reference station in Hawaii

APPENDIX B | PROOF OF VARIANCE BOUNDS FOR ARBITRARY LINEAR COMBINATIONS OF STATE-ESTIMATE ERRORS

Let us consider the following linear continuous-time system

\begin{matrix} \dot{x} = A x + B w & z = C x + v \end{matrix}

where x is the state vector, w is the process-noise vector, z is the measurement vector, and v is the measurement-noise vector. The Kalman-Bucy filter is given as follows:

\overset{\cdot}{\hat{x}} = A \hat{x} + K (z - C \hat{x})

Thus, the state-estimate error $\tilde{x} ≜ \hat{x} - x$ evolves as follows:

\overset{\cdot}{\tilde{x}} = (A - K C) \tilde{x} + B w + K v

Taking the Laplace transform yields the classical transfer-function representation:

\tilde{X} (s) = {\bar{H}}_{X / W} (s) W (s) + {\bar{H}}_{X / V} (s) V (s)

where:

\begin{matrix} {\bar{H}}_{X / V} (s) ≜ {(s I - A + K C)}^{- 1} B, & {\bar{H}}_{X / V} (s) ≜ {(s I - A + K C)}^{- 1} K \end{matrix}

are the transfer-function matrices from the input errors to the state-estimate errors.

Let the scalar output of interest be y = d^Tx, with $x, d \in ℝ^{m \times 1}$ . Then, we have the following:

\tilde{Y} (s) = {\bar{H}}_{Y / W} (s) W (s) + {\bar{H}}_{Y / V} (s) V (s)

where:

\begin{matrix} {\bar{H}}_{Y / W} (s) ≜ d^{T} {\bar{H}}_{X / W} (s), & {\bar{H}}_{Y / V} (s) ≜ d^{T} {\bar{H}}_{X / V} (s) \end{matrix}

In the proofs of Section 4, the “input error” can be any single element of the random vectors w or v. We assume that (i) the elements of w and v are mutually independent and (ii) the filters ${\bar{H}}_{X / W} (s)$ and ${\bar{H}}_{X / V} (s)$ are designed via upper bounds on the PSDs of w and v. Because the overall system is linear and d can be any vector, the proofs are sufficient to guarantee that the variance of $\tilde{y}$ —any scalar linear combination of the state-estimate errors—remains bounded for any linear combination of the elements of w and v.

Footnotes

↵1 The scalar function x(t) can represent any stationary sensor or process noise input to the KF. The function can also be the sum of multiple error components—e.g., $x (t) = \sum_{i} x_{i} (t)$ .

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

↵
Advanced Receiver Autonomous Integrity Monitoring Working Group C 3.0. (2016). ARAIM technical subgroup: Milestone 3.0 report [Working Group C: EU-US Cooperation on Satellite Navigation]. https://www.gps.gov/sites/default/files/2025-09/ARAIM-milestone-3-report.pdf
Google Scholar
↵
Bendat, J., & Piersol, A. (2010). Non stationary data analysis. In Random data analysis and measurement procedures (4th ed., 417–472). Wiley. https://doi.org/10.1002/9781118032428.ch12
Google Scholar
↵
Blanch, J., Walter, T., & Enge, P. (2019). Gaussian bounds of sample distributions for integrity analysis. IEEE Transactions on Aerospace and Electronic Systems, 55(4), 1806–1815. https://doi.org/10.1109/TAES.2018.2876583
Google Scholar
↵
Brown, R. G., & Hwang, P. Y. C. (2012). Mathematical description of random signals. Comparison of first and second order Gauss-Markov processes. In Introduction to random signals and applied Kalman filtering (4th ed., 326). Wiley.
Google Scholar
↵
DeCleene, B. (2000). Defining pseudorange integrity -overbounding. Proc. of the 13th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GPS 2000), Salt Lake City, UT, 1916–1924. https://www.ion.org/publications/abstract.cfm?articleID=1603
Google Scholar
↵
Gallon, E. (2023). High-integrity modeling of non-stationary Kalman filter input error processes and application to aircraft navigation [Doctoral dissertation, Illinois Institute of Technology]. https://www.proquest.com/openview/ea11612221576ce61bfe4c851833f45a/1?pq-origsite=gscholar&cbl=18750&diss=y
Google Scholar
↵
Gallon, E., Joerger, M., & Pervan, B. (2021a). Development of stochastic IMU error models for INS/GNSS integration. Proc. of the 34th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2021), St. Louis, MO, 2565–2577. https://doi.org/10.33012/2021.17962
Google Scholar
↵
Gallon, E., Joerger, M., & Pervan, B. (2021b). Robust modeling of GNSS tropospheric delay dynamics. IEEE Transactions on Aerospace and Electronic Systems, 57(5), 2992–3003. https://doi.org/10.1109/TAES.2021.3068441
Google Scholar
↵
International Civil Aviation Organization. (2001). International standards of recommended practices: Annex 10 to the convention on international civil aviation. https://www2023.icao.int/Meetings/anconf12/Document%20Archive/AN10_V2_cons%5B1%5D.pdf
Google Scholar
↵
Jada, S., & Joerger, M. (2021). High-integrity measurement error time correlation modeling using PSD-bounding of experimental data. Proc. of the Royal Institute of Navigation’s Navigation 2021 Conference, Edinburgh, Scotland, 238–246. https://doi.org/10.48448/kne6-z328
Google Scholar
↵
Joerger, M., Jada, S., Langel, S., Crespillo, O. G., Gallon, E., & Pervan, B. (2023). Practical considerations in PSD upper bounding of experimental data. Proc. of the 36th International Technical Meeting of the Satellite Division of the Institute of Navigation (ION GNSS+ 2023), Denver, CO, 441–448. https://doi.org/10.33012/2023.19196
Google Scholar
↵
Joerger, M., & Pervan, B. (2020). Multi-constellation ARAIM exploiting satellite motion. NAVIGATION, 67(2), 235–253. https://doi.org/10.1002/navi.334
Google Scholar
↵
Langel, S., Crespillo, O. G., & Joerger, M. (2020). A new approach for modeling correlated Gaussian errors using frequency domain overbounding. Proc. of the 2020 IEEE/ION Position, Location and Navigation Symposium (PLANS), Portland, OR, 868–876. https://doi.org/10.1109/PLANS46316.2020.9110192
Google Scholar
↵
Langel, S., Crespillo, O. G., & Joerger, M. (2024). Frequency-domain modeling of correlated gaussian noise in Kalman filtering. IEEE Transactions on Aerospace and Electronic Systems, 60(6), 8945–8959. https://doi.org/10.1109/TAES.2024.3442775
Google Scholar
↵
Langel, S., Khanafseh, S., & Pervan, B. (2014). Bounding integrity risk for sequential state estimators with stochastic modeling uncertainty. AIAA Journal of Guidance Control and Dynamics, 37(1), 34–46. https://doi.org/10.2514/L62056
Google Scholar
↵
Pervan, B. (2020). Ground-based augmentation system. In Y. Morton, F. Diggelen, J. Spilker, B. Parkinson, S. Lo, & G. Gao (Eds.), Position, navigation, and timing technologies in the 21st century. https://doi.org/10.1002/9781119458449.ch12
Google Scholar
↵
Pullen, S., & Joerger, M. (2020). GNSS integrity and receiver autonomous integrity monitoring (RAIM). In Y. Morton, F. Diggelen, J. Spilker, B. Parkinson, S. Lo, & G. Gao (Eds.), Position, navigation, and timing technologies in the 21st century. https://doi.org/10.1002/9781119458449.ch23
Google Scholar
↵
Rife, J., Pullen, S., Enge, P., & Pervan, B. (2006). Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Transactions on Aerospace and Electronic Systems, 42(4), 13861395. https://doi.org/10.1109/TAES.2006.314579
Google Scholar
↵
Sensonor. (n.d.). STIM-300 IMU product specifications Retrieved April 30, 2021 from https://www.sensonor.com/products/inertial-measurement-units/stim300/
Google Scholar
↵
Walter, T. (2020). Satellite-based augmentation systems (SBASs). In Y. Morton, F. Diggelen, J. Spilker, B. Parkinson, S. Lo, & G. Gao (Eds.), Position, navigation, and timing technologies in the 21st century. https://doi.org/10.1002/9781119458449.ch13
Google Scholar

PDF

Previous Next

High-Integrity Modeling of Nonstationary Noise Processes for GNSS/INS Integration

Abstract

1 INTRODUCTION

2 OVERBOUNDING, AUTOCORRELATION FUNCTION BOUNDING, AND PSD BOUNDING