Abstract
According to aviation minimum operational performance standards (MOPS), protection levels for satellite-based augmentation systems (SBASs) are to be computed as the product of the estimated standard deviation of errors and a scaling, or K-factor. MOPS recognized that K-factors were originally chosen to be consistent with certain assumptions that may not hold under all conditions. Considering the limited applicability of aviation-based K-factors, it will be important to identify a more rigorous method for deriving new SBAS applications for road, rail, or maritime use. Here we describe an innovative method applicable to any integrity risk (e.g., at 10−7 for aviation or 10−5 for maritime applications) and time interval T (e.g., 150 seconds to 1 hour for aviation or 3 hours for maritime use). This new method relies on rigorous probability justification. No restrictive assumptions are needed for the time correlation pattern of the errors. The method is easy-to-implement and applicable to any type of integrity risk or time interval T.
1 INTRODUCTION
Integrity as a theme is generally defined as “…a measure of trust which can be placed in correctness of the information supplied by a navigation system. Integrity includes the ability of the system to provide timely warnings to users when the system should not be used for navigation.” (US Federal Radionavigation Plan, 2008, p. B-4). This concept was initially developed and modeled within the scope of civil aviation for use in safety-of-life applications. Integrity is currently implemented in specific systems, including satellite-based augmentation systems (SBASs) or receiver autonomous integrity monitoring (RAIM), and is standardized by the International Civil Aviation Organization (ICAO; 1996) and minimum operation performance standards (MOPS) (2016). Today, there is a growing interest in the generalization of this concept beyond aviation as it may be critical for the safe navigation of maritime, rail, automotive, and unmanned aerial vehicles, among others.
The concept of aeronautic integrity relies on the definition of four specific features: the alert limit (AL); the time to alert (TTA); integrity risk (IR); and protection level (PL). The AL is the position error (PE) threshold that cannot be exceeded before triggering an alert. The TTA is the maximum allowable time elapsed from the loss of tolerance by the navigation system until the equipment issues the alert. The IR is the probability that the PE exceeds the AL beyond the TTA for at least one time interval T. The time interval T is adjusted to an appropriate value for a given targeted application. For example, at T = 150 seconds is used for precision aviation approaches. Finally, the PL is defined as the statistical bound such that the probability of the PE exceeding PL is smaller than or equal to the target IR. In practice, the PL is the product of the estimated errors standard deviation and a scaling factor known as the K-factor or K value. More precisely, horizontal PL (HPL) and vertical PL (VPL) are computed as indicated in Equation (1):
1
where dmajor corresponds to the error uncertainty along the semi-major axis of the error ellipse that depends on the geometry matrix and the error variance provided by SBAS, du is an over-bounding of the vertical error standard deviation, and KH, NPA, KH, PA and KV are the K-factors. K-factors can also be understood as a quantile of the distribution of the safety indexes (SFIs), where SFIs are positioning errors divided by their estimated standard deviations. Consequently, K-factors play crucial roles in maintaining SBAS service performance. Decreasing the K-factor will improve SBAS availability and continuity, but may lead to an unacceptable IR. Here, we focus on the computation of K-factors for the IR targets required by new SBAS applications, for example, 10−5/150 s for aviation approach Special Authorization Category I, 10−5/3 hours for maritime, and 10−9/1 hour for rail.
The MOPS (2016) has fixed K-factors for several aviation approaches, for example, KV = 5.33 for the category I (CAT I) vertical PL (IR = 10−7/150 s). Roturier et al. (2001) have explained the paradigms underlying these K-factors. However, MOPS (2016) also recognized that “…the values of KH,NPA, KH,PA, and KV were originally chosen to be consistent with certain assumptions on the distribution of position error and on error correlation time. It was then realized that these assumptions may not hold under all conditions, but that the choice of values is somewhat arbitrary” (p. J5).
Pervan et al. (2017) investigated the impact of the time correlation on the probability of false alerts and missed detection of ground based augmentation system (GBAS) monitoring tests. In particular, for ionospheric errors, The authors showed that, for ionospheric errors in particular, it is safer to consider more than one independent sample in a 15-second time window.
Bang and Milner (2021) recently evaluated the number of effectively independent samples over one hour for horizontal advanced RAIM using Monte-Carlo simulations of the Gauss-Markov processes. Their findings highlighted the need to challenge the MOPS assumption for many of the new integrity applications. Specifically, the authors reported that the number of effective samples per hour is greater than 10 when fixed by MOPS (2016); this leads to larger K factors for the Gauss-Markov process with correlation times equal to 100 or 500 seconds. We did not adopt the same methodology we would need to perform Monte-Carlo simulations for each correlation pattern and level of integrity; this would be prohibitive given the large number of integrity needs that will be encountered in maritime, rail, and road applications.
Here, we propose a rigorous method to compute K-factors suitable for any IR and time interval T. We refer to this method as Gauss-Markov K-factor determination (GMK). This method presents several distinct advantages:
GMK makes no restrictive assumptions on time correlation patterns;
GMK relies on rigorous formulations that apply to conditions associated with any integrity risk and time interval T.
The article is organized as follows. Section 2 presents the current aviation method (i.e., MOPS). Section 3 describes and justifies the new GMK method. Section 4 features the use of the GMK method in several different SBAS applications and compares it with MOPS for aviation users. Finally, Section 5 provides several conclusions and perspectives on this work.
2 STATE OF THE ART: MOPS K-FACTORS
Roturier et al. (2001) described the paradigms used to determine K-factors in MOPS (2016). Section 2.1 details the safety index (SFI) time correlation (or auto-correlation) pattern that is assumed implicitly by these paradigms. The MOPS K-factor formulae generated by Roturier et al. (2001) are explained in full in Section 2.2.
2.1 SFI Auto-Correlation Modeling
According to Roturier et al. (2001), the MOPS K-factors are based on crude modeling of the SFI correlation time. It is assumed that the SFI at time t and t + Δt are:
perfectly correlated if Δt is less than 360 seconds; and
independent of one another if Δt is greater than 360 seconds.
As explained by Roturier et al. (2001), the choice of 360 seconds was motivated by ionospheric corrections.
2.2 K-Factor Formulae
In this Section, we explain how to retrieve a K-factor formulation that justifies the MOPS K-factors described in Roturier et al. (2001). This formulation relies on two main assumptions:
the SFI auto-correlation pattern is compliant with the parameters described in Section 2.1; and
the SFI is a standard Gaussian vector.
The IR evaluated over time T is bounded by the probability of at least one misleading information (MI) over T, i.e., the probability that SFI exceeds the K-factor by at least one epoch. This bounding, while conservative, is quite stringent given a short TTA or high correlation time.
Normality and independence assumptions over 360 seconds lead to Equation (2):
2
where Xt denotes the SFI in one dimension at epoch t, Φ is the cumulative distribution function (CDF) of the standard normal and ⌈T / 360⌉ denotes a ceil of T / 360.
Because IR is close to zero, a Taylor series can be used to obtain the formula provided by Roturier et al. (2001) shown in Equation (3):
3
In the MOPS (2016), this formula has been applied to determine K-factors for vertical and horizontal cross-track positioning errors associated with precision approaches (PAs).
For non-precision approaches (NPAs), integrity must be ensured not only in the cross-track but in all horizontal directions. Consequently, IR is defined with respect to the Euclidian norm of the SFI in the horizontal plane and involves a two-dimensional SFI vector (e.g., one with both east and north components). The framework can be changed to include uncorrelated components in the SFIH vector using a Rayleigh distribution; the Euclidian norm of two centered normal variables with identical standard deviation follows a Rayleigh distribution. In practice, the standard deviations of the uncorrelated components of horizontal SFI may vary. However, using a conservative approach, it is possible to normalize horizontal SFI with the greatest standard deviation, i.e., the parameter dmajor described in annex J of MOPS (2016). These considerations and using similar reasoning for the vertical component will lead to Equation (4):
4
where ΦR is the CDF of the Rayleigh distribution.
In the MOPS (2016), Equation (4) has been applied to determine the K-factor for horizontal positioning errors in NPAs.
3 GMK METHOD: GAUSS-MARKOV FOR K-FACTOR DETERMINATION
Assuming stationarity, the auto-correlation of the SFI time series can generally be estimated with only a few hours of data. Unfortunately, these data are generally not available at the user level in real time because they involve real (or at least highly precise) positions to facilitate an evaluation of the errors. Alternatively, an SFI time series might be modeled beforehand to produce a K-factor based on the results. The Gauss-Markov (GM) process may be attractive because it is simple, familiar, and has been used to represent numerous real processes, for instance, in applications associated with finance and econometrics.
In practice, however, it may be impossible to guarantee a perfectly-fitted SFI time series as the result of a GM process (see Figure 1). Thus, our idea was to use the GM process to generate overbounding, rather than fitting, of the SFI auto-correlation pattern. We show here that overbounding guarantees integrity under comparatively mild assumptions. Section 3.1 details these results as they apply to SFI auto-correlation modeling.
Once it becomes clear that modeling is justified, K-factors must be computed for the chosen GM process. Section 3.2 describes the approach used to identify a convenient equation. First, we assessed the most straightforward formula, which was one derived from the GM auto-correlation function. The implementation of this formula revealed numerical instabilities for a time interval T greater than a few minutes and for an IR interval at 10−7. To the best of our knowledge, no alternative formula(s) hve been documented for the GM process. Interestingly, De Long (1981) provided both exact and approximate formulae for the tail probability of the Ornstein-Uhlenbeck process maximum over a time interval T. We show here that these results can be used to compute the K-factor associated with the GM process. We validated and compared both the exact and the approximate formulae and recommend an approximation, which was more robust in cases with a low IR or large time interval T that appears to be fully convenient for use in the SBAS context because it is explicit and easy to implement.
To summarize, the proposed GMK method includes two steps:
Overbounding of the SFI auto-correlation pattern via a GM process with correlation time τ, as will be explained in Section 3.1 and illustrated in Figure 2;
Evaluation of the K-factor calculated for the GM process with correlation time τ, as will be explained in Section 3.2.
3.1 SFI Auto-Correlation Modeling
The auto-correlation pattern of the SFI plays a critical role in the determination of K-factors and overall SBAS performances. Thus, it will be important to identify a correlation pattern that is realistic with respect to the SFI time series. The correlation pattern should also permit a tractable determination of K-factors.
The GM process is very attractive because it is very simple and familiar to most end-users as GM process or auto-regressive functions with order 1. These processes are defined below.
Definition 1: Gauss-Markov (GM) process
The standardized GM process (Yt) with correlation time τ can be defined by Equation (5):
5
where time , εt follows a centered normal law with variance (1 – e−2/τ).
This computation is straightforward given that the variance of Yt is 1 and that the GM process has an exponential auto-correlation function as shown in Equation (6):
6
As shown in Figure 1, the GM correlation pattern is more realistic than the MOPS model for the SFI time series because the former assumes perfect correlation followed by independence. However, the GM correlation pattern is not fully adequate. Typically, the SFI auto-correlation function will decrease more slowly than the exponential rate or will display some correlation peaks due to the periodicity of the SBAS corrections.
Theorem 1 shows that the assumptions associated with the SFI correlation pattern can be relaxed due to Power Spectral Density (PSD) overbounding as defined below.
Definition 2: PSD overbounding
For any pair Xt and Yt of stationary Gaussian processes, Xt is said PSD over-bounded by Yt if, for all frequencies (f), PSDX(f) ≤ PSDY(f).
Theorem 1
If Xt and Yt are two stationary Gaussian processes such that Yt PSD over-bounds Xt, then, as shown in Equation (7):
7
Demonstration
We note that rX and rY are the auto-correlation functions (ACFs) of Xt and Yt, respectively.
The vectors X = (Xt)t=0,⋯,T and Y = (Yt)t=0⋯T are Gaussian random vectors with mean 0T and correlation matrices defined by Equation (8):
8
We know from Langel et al. (2020a) that the PSD overbounding condition implies that ∑Y ≥ ΣX in the positive definite sense (i.e., that ∑Y – ΣX is positive semi-definite matrix). Based on Corollary 4 described by Anderson (1955), the implication is that P(maxt =0, ⋯, T |Yt| ≤ K) ≤ P (maxt =0, ⋯, T |Xt| ≤ K).
Theorem 1 implies that a K-factor computed according to a process Yt (i.e., that satisfies P(maxt =0, ⋯, T |Yt| > K) = IR) guarantees IR for the process Xt, if Yt is PSD overbounding Xt. Consequently, this theorem justifies, under the usual assumptions (i.e., normality and stationarity of the SFI time series), the use of the GMK approach based on two steps:
SFI PSD overbounding with a GM process,
K-factor computation for the overbounding GM process.
Theorem 1 assumes the normality of the SFI time series distribution. In practice, this assumption is not always satisfied. Typically, the SFI distribution might display bias or a non-Gaussian tail. Conditions under which the normality assumption in PSD overbounding might be relaxed is an area of active research because it also concerns positional integrity based on the Kalman algorithm. Recently, Langel et al., (2020b) proposed a method to relax the normality assumption based on a combination of symmetric and PSD overbounding. Some stringent assumptions are still required by the overbounded process, including that the functions must be zero-mean, symmetric, unimodal, and expressed as a linear transformation of a spherically-symmetric random vector. A similar methodology that could permit us to relax the normality assumption for GMK factors will be considered and developed in future work.
Note that the normality assumption is also applicable to the MOPS K-factors, even if the MOPS correlation pattern (i.e., independence every 10 minutes) is satisfied. More precisely, the CDF overbounding concept, which was introduced by DeCleene (2000) to justify Gaussian positioning errors when range errors are not Gaussian (but overbounded by a Gaussian distribution in the CDF sense), is not directly applicable to K-factor computation. Indeed, while the CDF overbounding theorem ensures the stability of any linear combination, the maximum absolute value would need to remain stable for the computation of a K factor.
The next section focuses on the second step of the approach, i.e., the determination of K-factors for the GM process.
3.2 K-Factor Formulae
This section addresses the determination of K-factors for a GM process Yt with correlation time τ.
A full numerical evaluation based on Monte-Carlo simulations of GM processes would be too computationally burdensome for an IR of 10−7 over one hour. This would require the simulation of a GM process for at least 107 hours, which is more than a thousand years.
When the time interval T is a few seconds, the K-factor can be computed based on the Gaussian vector Y = (Yt)t=1…T as shown in Equation (9):
9
where ΦY is the CDF of the Gaussian vector Y with dimension T, mean 0T, and covariance matrix V such that Vi,j = e|i-j|τ for i, j = 1,.,., T.
In practice, this formula is numerically intractable when the Gaussian vector’s dimension T reaches 150 or 3600 seconds. In that case, we propose to consider a GM process as a realization of its time-continuous counterpart known as an Ornstein-Uhlenbeck process, which is described below.
Definition 3: Ornstein-Uhlenbeck process
A standardized Ornstein-Uhlenbeck process is a continuous-time Gaussian process with mean 0, variance 1, and auto-correlation function equal to e−|Δt|, where .
The following theorems document how K-factors of GM processes can be bounded by the K-factors of standardized Ornstein-Uhlenbeck processes.
Theorem 2
If Yt is a standardized GM process with time correlation τ and Zt is a standardized Ornstein-Uhlenbeck process, then, as shown in Equation (10):
10
Demonstration
Let us consider the time-continuous process as equal to a time-rescaled Zt as shown in Equation (11):
11
For any is a GM process with a correlation time equal to τ because of the two following arguments:
is a standardized Gaussian random variable for any ,
because Zt is an Ornstein-Uhlenbeck process.
Because maximum is not necessarily reached on discrete times t, as in Equation (12):
12
Theorem 3
If Y1,t and Y2,t are independent standardized GM processes with time correlation τ, then Z1,t and Z2,t are independent standardized Ornstein-Uhlenbeck processes. As per Equation (13):
13
Demonstration
The demonstration relies on time rescaling as shown in Theorem 2.
The final step is to derive the K-factor for a standardized Ornstein-Uhlenbeck process. As noted by Azaïs & Wschebor (2002), the distribution of maxt∈[0,T/τ] |Zt| and is explicit as was described in De Long (1981). This theoretical result is illustrated in Figure 3. A total of 104 Monte Carlo simulations of maxt=0, …,T |Yt| and were performed (with a correlation time equal to 1000 seconds and T equal to 3600 seconds). As expected, the empirical quantiles of order (1 – IR) are very close to the GMK factors evaluated by Equation (4) in De Long (1981) for any IR from 0.01 to 0.9.
Unexpectedly, some numerical instabilities occur for IRs lower than 10−7 as shown in Figure 4. One can circumvent these instabilities by using the asymptotic formula provided by De Long (1981 p. 2205): for large values of T /τ and K, shown here in Equation (14):
14
where Γ is the Gamma function and Yt is a q-dimensional vector of independent GM processes with correlation time τ (here q equals 1, typically for the vertical component, or 2, for the horizontal component). As shown in Figure 4, the asymptotic formula presents a reliable evaluation of the GMK-factors when the exact formula cannot be used because of numerical problems (i.e., when GMK factors are >6 as shown in Figure 4). Likewise, the findings shown in Figure 4 reveal that the asymptotic formula is highly accurate once the K-factor is >4, which is within the typical range of interest for SBAS applications; this is the case even when T / τ = 1. Thus, this asymptotic formula is fully adequate for SBAS applications.
4 APPLICATIONS
This section describes the application of GMK-factors in typical SBAS contexts.
First of all, it is important to recognize that we have not deliberately fixed the correlation time τ. Theorem 1 provides some degrees of freedom on the choice of τ: this choice is safe with respect to IRs once the GM process with correlation time τ is overbounding the SFI process. Consequently, it is not necessary to perform a precise fitting of τ with the SFI time series. The choice of τ could be based on a trade-off between SBAS availability and simplicity at the system/user level under a given integrity constraint.
Our first application concerns the aviation approaches currently supported by SBAS. Shown in Table 1 are the integrity requirements for these various aviation approaches.
Figure 5 displays the GMK-factors obtained for these aviation approaches and compared them with the analogous MOPS K-factors. As expected, GMK-factors decrease when the correlation time τ increases. For PA-V and NPA-H approaches, GMK-factors are close to the MOPS K-factors when the correlation time τ is high (typically >1 hour). However, the MOPS K-factors seem a bit optimistic because the GMK-factors remain higher than the MOPS K-factors even when τ = 2 hours.
Because of the generalizability of the GMK method, it is also possible to derive K-factors for new SBAS applications (Figure 6) including
aviation Special Authorization CAT I approach (vertical IR equals to 10−5/150 s); and
maritime (horizontal IR equal to 10−5/ 3h).
Finally, the figures below provide GMK-factors for any correlation time τ up to two hours and IR between 10−4 and 10−9 over one hour (T = 3600 seconds) for:
a single component, typically for vertical or cross-track errors, as shown in Figure 7;
and the norm of two independent components, typically for horizontal errors, as shown in Figure 8.
Another possible application for GMK-factors could be for integrity monitoring tests at the system level, for example, SBAS position domain checks or the GBAS DSIGMA/CDD tests described by Pervan et al. (2017). These monitoring tests involve detection thresholds that are defined similarly to the K-factor. For example, in fault-free conditions, the statistical test is expected to remain under the detection threshold throughout the exposure time with a probability that is greater than that of a false alarm (PFA). Pervan et al. (2017) propose an analytic method based on the threshold-crossing problem to derive the PFA given a fixed threshold. The GMK method could provide an alternative evaluation of this probability. It would be interesting to compare these methods both numerically and theoretically with particular attention to the assumptions associated with the time correlation pattern.
5 CONCLUSION AND PERSPECTIVES
This paper describes a rigorous method to compute K-factors for any IR (e.g., 10−5, 10−7, and others) and time interval T (e.g., 150 seconds, 1 hour, or 3 hours) that is consistent with a common assumption of stationarity and Gaussianity of the safety indexes.
The proposed method was divided into two steps. The first step was the overbounding of the SFI PSD using a GM process with correlation time τ. We demonstrate that overbounding guarantees that a K-factor computed for a GM process bounds IR even if the SFI is not a GM process. In practice, this overbounding step could be performed once for all at system level and documented in the standards. The second step is the computation of K-factors for a GM process with correlation time τ. To the best of our knowledge, no explicit formulations are available for GM processes; thus, we used theoretical results that focus on the maximum of its time-continuous counterpart known as the Ornstein-Uhlenbeck process. We provide easy-to-implement formulations for IR related to errors in one (e.g., vertical or horizontal cross-track errors) and two dimensions (e.g., horizontal errors). These formulations have been validated with respect to alternative simulations and formulations. These factors are then applied to several new integrity needs, for example, aviation Special Authorization CAT I, maritime, and any IR per hour between 10−4 and 10−9.
The new method proposed for K-factor computation has several strengths. It relies on rigorous theoretical grounding with solid probability justification and is thus applicable to any integrity requirement. No restrictive assumptions are needed on the time correlation pattern of the errors. Finally, the method is easy-to-implement and applicable to any IR and time interval T.
The main limitation of the method is probably the assumption of normality of the reduced errors. Methods that may be used to relax this assumption are currently under investigation, with a specific interest in the methods emerging that address the integrity of the Kalman algorithm. Other topics for future study include the characterization of correlation times for potential SBAS applications, including those needed to support maritime, unmanned aircraft, rail, and road travel.
HOW TO CITE THIS ARTICLE
Antic, J., Maliet, O., Trilles, S. (2023). SBAS protection levels with Gauss-Markov k-factors for any integrity target. NAVIGATION, 70(3). https://doi.org/10.33012/navi.594
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.