Statistical characterization of composite protection ... - Semantic Scholar
GPS Solut DOI 10.1007/s10291-010-0188-2
ORIGINAL ARTICLE
Statistical characterization of composite protection levels for GPS Dean Bruckner • Frank van Graas Trent Skidmore
•
Received: 14 October 2009 / Accepted: 20 September 2010 Springer-Verlag 2010
Abstract A statistical characterization is presented of Global Positioning System (GPS) user range error as a normally distributed random variable with non-zero mean over the length of the aircraft precision approach operation, correlated from one GPS measurement epoch to another and from one satellite to another. This leads directly to modeling GPS error in the position domain as multivariate normal with non-zero mean. Based on this model, a vertical composite protection level VPLc and a horizontal composite protection level HPLc are each calculated as scalar values from a univariate normal distribution displaced from the origin by the worst-case position domain bias combination possible, given the maximum possible individual satellite bias magnitudes and the satellite geometry. A method is then presented by which exact values—that is, values accurate to a user-defined error tolerance and consistent with statistical assumptions—of VPLc and HPLc are obtained, and by which computationally efficient approximations may be evaluated. A statistical quadratic form under the multivariate normal distribution is used to derive a new class of protection levels based on the probability enclosed within a radius defined in two or more dimensions. A central Chi-square representation of this quadratic form is also presented and is incorporated into a six-step computational procedure for the two-dimensional composite radial protection level RPLc. This procedure is extended to the spherical protection level (SPLc) and the ellipsoidal protection level (EPLc).
D. Bruckner (&) F. van Graas T. Skidmore Ohio University, Athens, OH, USA e-mail: [email protected] URL: www.ohio.edu/avionics
Keywords VPL HPL Composite protection level Radial protection level Ellipsoidal protection level Conspiring bias
Introduction Setting a protection level (PL) is a widely accepted integrity methodology for navigation systems, particularly in safety-intensive aviation applications. A PL is commonly defined as the maximum navigation system error (NSE) from the true position that can occur in a specified direction such as the vertical, ‘‘guaranteed’’ to a sufficiently high level of probability. This section introduces competing approaches to setting PLs and provides context for the statistical characterization of the present effort. Protection levels in aviation safety Integrity specifications in the GPS-based Local Area Augmentation System (LAAS) and Wide Area Augmentation System (WAAS) call for calculating vertical, horizontal and lateral (cross-track) PLs, which are abbreviated as VPL, HPL and LPL, respectively (RTCA 2004, 2006, 2007). GPS Receiver Autonomous Integrity Monitoring (RAIM) algorithms likewise compute a VPL or HPL. In these types of applications, the PL is continually compared to an appropriate Alert Limit (AL) throughout the operation. If the PL exceeds the AL, an alarm is annunciated to the user. This typically results in the navigation system being declared unusable for the intended purpose until the PL again drops below the AL. In the RAIM context, system design emphasis is usually placed on detecting a satellite malfunction, typically marked by a growing pseudorange bias, before that bias becomes large enough to pose a
123
GPS Solut
serious problem. In the present effort, however, it is assumed that a certain level of pseudorange bias exists in normal GPS operation. This point has been illustrated in recent years in various studies of residual group delay bias in GPS antennas and antenna arrays (Murphy et al. 2007). The composite PL is an attempt to treat this bias explicitly to improve overall differential GPS performance. Integrity for differential GPS systems such as WAAS and LAAS has traditionally been assured by over-bounding both the bias-like errors, i.e. persistent and slowly varying but not necessarily growing non-zero mean errors, and noise-like errors seen in GPS measurements by using a single closed-form probability distribution. Typically, this has been the normal distribution. A different approach has also been advanced, one that explicitly assigns probability distributions separately to noise and to biases. The PL that results from this separate treatment is a composite of noise and bias bounds. Previously defined composite vertical protection level definitions, such as VPL–C (van Graas et al. 2004) and VPLbias (Rife and Pullen 2005b), are prominent examples. In the present effort, a statistical method for computing a new vertical composite protection (VPLc) and horizontal composite protection level (HPLc) is defined. A new class of composite protection levels is also developed from a quadratic statistical form. This class includes the twodimensional radial protection level (RPLc), as well as the three-dimensional spherical protection level (SPLc) and ellipsoidal protection level (EPLc). Previously defined composite PLs mentioned above are seen to be approximations to the PLs presented here. As others have done before, this effort examines PLs only under the H0, or fault-free, hypothesis (Rife and Pullen 2005a; Shively and Hsiao 2004). This is appropriate for the analysis of a generic PL formulation not closely tied to the configuration of any particular system such as LAAS or WAAS. It is also appropriate in order to permit a focus on improving system accuracy performance through error reduction. The H0 hypothesis governing the present analysis is simply normal GPS system operation, with the additional proviso that some level of pseudorange bias exists in normal GPS operation, as asserted above. The mere presence of such bias, therefore, does not necessarily lead to the rejection of this H0 hypothesis and the subsequent declaration of a system fault. The present effort is motivated largely by repeated examination of slowly varying biases in flight test data recorded in differential GPS-based precision approach and landing operations over a 15-year period (van Graas et al. 2004, Skidmore et al. 1996). The presence of persistent position domain biases that remain essentially constant over an individual aircraft precision approach but vary between approaches, as well as the shortcomings in
123
traditional PL techniques in adequately representing these biases, justify exploring a new technique. Such a technique must be able to separate differential GPS spatial and temporal errors cleanly into Gaussian and non-Gaussian components that may each be bounded. The authors present such techniques elsewhere, along with flight test results that demonstrate significant decreases in the size of composite PLs when compared to the traditional PLs used in LAAS (Bruckner et al. 2010a). Conceptually, decorrelation between reference station and aircraft observations due to localized tropospheric and ionospheric effects, if not removed by ionospheric-free processing, can be accounted for by inflating the bias term bound. Yet another approach could be taken using Extreme Value Theory (EVT). This method provides advantages in error model analysis for unknown distributions. However, since an EVT-produced probability distribution is of an unknown type, it does not seem useful for setting PLs in the high-integrity precision approach and landing context. For this reason, the EVT method is not considered or further explained here. A case study illustration We present a case study illustration of a general statistical model previously developed by others (Ruben 1962; Ruben 1963; Kotz et al. 1967; Sheil and O’Muircheartaigh 1977). The model is applied here to the linearized solution method of the GPS equations in a new way. It is known that the modeling of GPS pseudorange error as unbiased, uncorrelated, Gaussian random variables (RVs) has great practical value. Here, that model is extended by redefining the error on each pseudorange as a Gaussian RV correlated with the Gaussian RVs that describe errors on the other pseudoranges. This is intended to model both the observed correlation in biases and the observed absence of correlation in noise between satellites. These Gaussian RVs are then employed in a quadratic form to define composite protection levels that treat noise and bias separately for the differential GPS-based precision approach and landing operation. The advantage of this method, which is not well known in the navigation community, is readily apparent in the resulting definition of protection levels that are less conservative than those constructed for the unbiased model. This in turn benefits system availability and continuity. The novelty of the present effort is threefold: the use of a quadratic form to define composite PLs that treat noise and bias separately in the context of differential GPS for precision approach and landing; the expressions provided for propagating noise using non-traditional weighting matrices; and the definition of range domain error for a particular satellite as a piecewise fit over time
GPS Solut
of a succession of Nðli ; r2i Þ-distributed RVs, which are correlated with the RVs representing the other satellites. It is hoped that placing these equations for composite PLs into a single document will make this approach more accessible to the navigation community, particularly when evaluating the effectiveness of computationally efficient approximations to these more exact PL expressions. Due to lack of space, such approximations are presented elsewhere by the authors in a follow-on paper (Bruckner et al. 2010b).
Designation of random variables The random variables are briefly defined in this section. A few well-known expressions such as the GPS measurement model are briefly restated to highlight novel points in the new technique.
It is assumed that a differential GPS system is operating normally. A linearized GPS measurement model is also assumed. Although this derivation is not limited to LAAS applications, it begins with the same GPS measurement model as defined for that system (RTCA 2004). The east– north–up locally level coordinate frame assumed here is a slight modification to that definition, in which the LAAS horizontal coordinate axes are either aligned with the runway or are arbitrary. The GPS measurement model is stated briefly as follows: ð1Þ
where Dy is the (M 9 1) vector of ranging measurements, minus expected ranges based on known satellite positions, assumed user position and user clock offset; H, the (M 9 4) observation matrix, in which the first three elements of each row form a unit vector from the user to the satellite in the locally level coordinate frame, and the fourth element is 1; Dx, the (4 9 1) vector from the user’s assumed position to the user’s true position in the locally level frame, including time as a fourth dimension; and e, the (M 9 1) vector of ranging measurement errors, to be defined further below. The difference between assumed user position and true user position may be estimated as follows, through the use of multiple linear regression techniques that employ weighted least squares: D^ x ¼ S Dy where
Weighting matrix W 2 2 r1 0 6 0 r22 6 W1 ¼ 6 . .. . . 4 .. . . 0 0
ð3Þ is defined as follows: 3 0 0 7 7 .. 7 . 5
ð2Þ
ð4Þ
r2M
where r2i is the variance of the unbiased and uncorrelated Gaussian range domain error of the ith satellite measurements that results from ground receiver noise and multipath, airborne receiver noise and multipath, and residual delays from ionospheric and tropospheric gradients between the reference station and the user. Matrix S is the projection matrix from the range domain to the position domain. The first modification to the GPS measurement model is to characterize user range domain error vector e using M normally distributed, discrete time random processes: ui;k Nðli ; r2i Þ
The GPS measurement model
Dy ¼ H Dx þ e
S ¼ ðHT W HÞ1 HT W
ð5Þ
where ui;k is the ith random sequence that represents a continuous time random variable, uniformly sampled at period DT seconds; li , the mean of the ith random sequence; r2i , the variance of the ith random sequence; i = 1, 2,…, M is the satellite index; and k = …, -2, -1, 0, 1, 2 … is the sample index and time variable. These processes are assumed to be discrete time Gaussian stationary with constant mean and variance over the duration of an aircraft precision approach and landing operation. Accordingly, the derivations presented in this section are equally valid for each GPS measurement epoch. Therefore, only the current epoch, where sample subscript k = 0, is considered. The sample subscript k is suppressed. Range domain random variables Now let M-dimensional real random vector u consist of correlated RVs u1 ; u2 ; . . .; uM , and let these RVs characterize range domain error for the duration of the precision approach and landing operation. It follows from the definition of ui;k above that each of these RVs is normally distributed with parameters mean li and variance r2i . These parameters are themselves constant over the operation, but they may differ from one satellite to another during this time interval. The assumption of correlation models underlying correlation of pseudorange biases between satellites, such as can occur due to group delay when the observed angle between two satellites, as apparent to a user’s antenna, is small. The pseudorange noise of each satellite is still assumed to be uncorrelated with that of other satellites. In the formation of PLs, the assumption of
123
GPS Solut
correlated biases between satellites is taken a step further by intentionally correlating them perfectly to produce the maximum possible position domain bias to make up for not knowing the probability distribution of the biases. When compared to the LAAS error model presented above, this definition is seen to be a departure from conventional treatment of range domain user measurement error, in which u1 ; u2 ; . . .; uM are considered uncorrelated, zero mean, identically distributed normal variables. The conventional model has been widely accepted not because it fits GPS behavior perfectly, but rather because it has worked without imposing a great deal of complexity for minimal added performance gain (Misra and Enge 2004). However, new GPS user requirements are pressing the system to its limits and, as noise filtering is increasingly effective, the presence of biases is receiving greater recognition. Therefore, long-standing assumptions are being reexamined by researchers in this field. This paper constitutes one such reexamination. The model proposed here is motivated by the desire to make practical use of the phenomenon of biases in user GPS measurements that persist essentially unchanged over the five- to ten-minute duration of a precision approach, but may change noticeably by the beginning of the next approach. Indeed, GPS range domain measurements are known to exhibit both correlation over time and correlation between simultaneous measurements from different satellites (Misra and Enge 2004). Three alternate range domain error models were considered besides the one presented here. In the first case, treating range domain error as a linear process driven by white noise has the virtues of constant mean and variance, thus ensuring both wide-sense stationarity and normally distributed error in both the range domain and position domain. Furthermore, realizations of this process may legitimately take on successive values distant from the mean for an extended period of time, relative to the short duration of a precision approach. This behavior can therefore model a slowly changing bias, while remaining assured of convergence (Chatfield 2004). However, the ability to translate the current deviation from the constant process mean, i.e. the instantaneous bias, into the position domain, and there to make statistical sense of it, remains out of reach. A second model was described in (Skidmore et al. 1996), in which the convolution of uniform and Gaussian distributions represents the addition of bias-like and noise-like range domain errors, respectively. This description is straightforward, but does not yield a normal distribution in the position domain. A third model, defined in (Rife et al. 2006), accounts for range domain biases by convolving an identical pair of distributions, e.g. normal distributions when applied to GPS PLs, with each other, where one member of the pair is offset from the origin by
123
the positive magnitude of the range domain bias, and the other by the negative magnitude. The present effort can be viewed as a limiting case of the paired distribution approach. If the identical distributions represent a probability assignment of 50% to the positive maximum of the bias and 50% to the negative maximum, accomplished through the sifting property of convolution with the delta function, the approach in this paper is to assign 100% of the probability to one bound, whichever one has a greater magnitude, and zero to the other. If the two bounds have the same magnitude, then which one is selected does not matter. In any case, the PL’s symmetry about zero ensures a worst-case bound. In summary, the user range error u is defined here as normally distributed, but over the limited time horizon of a typical aircraft precision approach and landing operation. Mathematically, this model represents a piecewise fit of a succession of Nðli ; r2i Þ-distributed RVs over time. The parameters li and r2i for the ith satellite may change from one subdomain, or time segment, to another. However, within the kth subdomain, li and r2i are assumed to be constant. The quantities li and r2i may also be described as piecewise continuous functions, each of which consists of a series of steps. These steps occur on transition from one subdomain to the next. Given the frequent jumps in protection levels seen during a typical precision approach as currently implemented, discontinuities of bounded magnitude in li and r2i at subdomain boundaries do not appear to present a problem to aircraft operations. The value of bias li for a given satellite is generally unknown, but is contained within the symmetric bounds that are described above and that may decrease over time in some processing techniques. Such techniques may also estimate the bias value itself with increasing accuracy over time. Range domain errors are assumed to be correlated between satellites, also as described above. Correlation across time does not enter directly into the PL calculation performed at each GPS measurement epoch and is not considered further. Since correlated variables are already assumed in the PL derivation that follows, the assumption of correlation here does not invalidate it. One drawback with this assumption is that the position domain Dilution of Precision (DOP) values are now correlated and their usefulness as individual scalar parameters declines. Position domain random variables The position domain error, defined as the difference between estimated position offset D^x and true position offset Dx from an assumed position, may be characterized as a four-dimensional real random vector q with
GPS Solut
component RVs x, y, z and t. These are east, north and up in a locally level coordinate system and user receiver clock offset from GPS time. The random vector q has a quadvariate normal distribution, i.e. multivariate normal with four dimensions, inherited from multivariate normal vector u via linear transformation, as defined 2 3 2 3 u1 x 6 u2 7 6y7 6 7 7 q ¼ S u or 6 ð6Þ . 7 4z5 ¼ S6 4 .. 5 t uM where projection matrix S, defined in (3), is a deterministic linear transformation. Then RV x, which is the position domain error in the east direction, is written as x¼
M X
s 1 ; i ui
ð7Þ
i¼1
where s1 ; i is the element in the first row and ith column of matrix S. The RVs y, z and t are defined similarly. It is assumed that the elements of matrix S are known with certainty. Therefore, any actual geometry errors due to ephemeris inaccuracies, user position error and the like may be projected onto the line of sight between the assumed user position and the calculated satellite position, and then combined with the pseudorange bias error. Thus, the elements of matrix S may be treated deterministically as constants instead of as additional random parameters. In LAAS development, a great deal of effort was devoted to ensuring that large ephemeris errors would not result in differences in the projection of errors between ground and air. Whenever this algorithm is implemented, the reference station should be assigned to monitor this error source. Mean and variance of vector q The mean vector l and variance–covariance matrix R of vector q may be derived using full matrix notation and beginning directly with RVs x, y, z and t (Misra and Enge 2004; Therrien 1992). They may also be computed using the simple linear transformations of the GPS measurement model and the well-known theorem on the numerical characteristics of a linear combination of RVs (Wackerly et al. 2002, see Theorem 5.12). Covariance propagation from the range domain to the position domain is also treated extensively in a number of texts (Leick 2004) and is not repeated here. Modeling noise propagation: unweighted and weighted least squares If noise sigma values are not considered identical for all the satellites, then the propagation of non-identically distributed noise must be modeled when constructing the
variance–covariance matrix, which in turn is used to set the protection level. Such is the case in this statistical model. This modeling happens automatically in the equation for VPL for the LAAS H0 hypothesis, in which the least squares projection matrix S is formed using weighting terms 1 r2i and multiplied by the squares of the corresponding elements of the S matrix. The equation for VPLH0 assumes no bias and is thus appropriate to illustrate the propagation of noise; this equation is as follows (RTCA 2004), vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M uX VPL ð8Þ ¼ K t s2 r 2 H0
ffmd
3;i i
i¼1
where Kffmd is the fault-free missed-detection sigma multiplier. If, however, the unweighted least squares method is used for position solutions rather than the weighted method, then multiplying r2i by the squares of the elements of the unweighted S matrix to create desired terms of variance–covariance matrix R no longer correctly models the propagation of non-uniform noise. In this case, the revised variance–covariance matrix is instead obtained using notation G for the unweighted projection matrix (Suddapalli 2004), R0 ¼ G Q GT
ð9Þ
where the variance–covariance matrix of range domain error u is 2 2 3 r1 0 0 6 0 r22 0 7 6 7 Q¼6 . ð10Þ .. 7 .. . . 4 .. . 5 . . 0
0
r2M
and G ¼ ðHT HÞ1 HT
ð11Þ
As described above, noise is assumed to be uncorrelated from one satellite to another. In this case, the corrected VPL equation for the bias-free case would use the (3, 3) element of matrix R0 as pffiffiffiffiffiffiffiffiffiffiffi ð12Þ VPLH0 ¼ Kffmd R0 3 ; 3 If position solutions are obtained using a weighting scheme other than either unweighted least squares or least squares weighted by 1 r2i , then an additional correction is required. For the case of least squares weighted by 1=li , for example, the inverse of the weighting matrix is 2 3 l1 0 0 6 0 l2 0 7 6 7 W1 ¼ 6 . ð13Þ .. 7 .. . . 4 .. . 5 . . 0
0
lM
123
GPS Solut
and the S matrix, as before, is defined as T
1
S ¼ ðH W HÞ HT W
ð14Þ
distribution in the position domain with mean b and standard deviation r, defined as Zq
n o 1 pffiffiffiffiffiffi exp ðr bÞ2 =ð2r2 Þ dr r 2p
Then, the revised covariance matrix that correctly models the propagation of noise through the least squares method weighted by 1=li is
p¼
R00 ¼ S Q ST
In this PL context, position domain bias b is set to the bias bound, as in Fig. 1. The original VPLc value, computed as b ? Kr and shown in Fig. 1, is an artifact of the protection level for the zero-mean distribution from which it was derived. Because the PL must be symmetric about the origin, the integral of the probability density function with integration limits (b ? Kr) and b ? Kr contains a significantly larger probability content than for the original integration limits -Kr and ?Kr. Note that the interval [-Kr, Kr] for the zero-mean, i.e. unbiased case is the same size as the interval [b - Kr, b ? Kr] for the biased case, even though it has been translated from the origin. This translated interval, however, cannot be used as a protection level due to its asymmetry about the origin. If instead the value of q can be found numerically by varying the integration limits until the desired value of p is achieved, then an exact probability may be enclosed by a new VPLc, to the specified error tolerance and consistent with the assumptions. The same approach is applied to define HPLc. The quantity HPLc, like the original HPL, is defined under a univariate normal distribution as the projection of horizontal errors onto the major axis of the error ellipse in the horizontal plane. Because the semimajor axis length of an ellipse is its largest radius, the HPL and HPLc are worstcase measures of generic lateral (cross-track) error and are used when the orientation of a runway is unknown. The HPL is not properly a radial measure, in which a specified probability is enclosed within a given radius. In fact, an HPL interpreted as a radial bound will always contain a smaller probability than the same HPL value interpreted as a linear bound. Figure 2 illustrates the errors
ð15Þ
To avoid a weighting matrix that is not positive definite, minimum and maximum limits are placed on the values of li . Use of these limits prevents algorithm instability while not noticeably affecting results. Since noise and bias are treated separately in forming the composite PL, the noise propagation techniques described in this section remain valid beyond the bias-free case shown here. Bias uncertainty and formation of the protection level An estimator B^i ðtÞ may be developed for the bias li in pseudorange measurements from the ith satellite at time t. The notation B, more familiar than l in GPS contexts, is used for range domain bias hereafter. Although the estimator B^i ðtÞ is not specified in this paper, it is assumed that whatever sampling distribution is selected constitutes a 100% confidence interval for B^i ðtÞ, centered at zero with upper and lower bounds Bmax ; i . One approach assumes a uniform sampling distribution (Skidmore et al. 1996), while another would assume a sinusoidal distribution more heavily weighted toward the upper and lower bounds. Because the distribution is here assumed to be unknown except for the stipulation that it is contained within bounds Bmax ; i , the non-Gaussian portion of the composite PL is formed using a worst-case vector addition of all the Bmax ; i . The bias bound and the Gaussian noise bound are added to form the PL. The advantage of treating the noise and bias separately is primarily that the noise bound now contains only thermal noise and is thus much smaller than when other terms are included. Furthermore, the decrease in the noise bound typically exceeds the increase in the bias bound, resulting in a net decrease to the size of the PL.
q
ð16Þ
Statistical method for exact VPLc and HPLc When viewed within the statistical framework described above, the ‘‘original’’ composite VPLc expressions noted in the literature are actually conservative approximations for a true value that can be derived statistically. This true value q is the magnitude of the upper and lower integration limits, which are symmetric about zero, of the integral of the probability density function of a normal distribution having a non-zero mean and a given probability p. Consider such a
123
Fig. 1 Probability density function for a biased univariate normal distribution in the position domain, with original and new VPLc shown. Mean value b is set to the bias bound
GPS Solut
Fig. 2 Errors enclosed within ‘‘linear’’ HPL but not enclosed within ‘‘radial’’ HPL of the same magnitude
projected onto the ‘‘linear’’ HPL but not enclosed within the ‘‘radial’’ HPL, using data from a Monte-Carlo simulation created by the authors. These points between the straight line and circle, plotted with a larger symbol than the other points, are shown for the [B - 3r, B ? 3r] case for illustration purposes. As the sigma multiplier K increases, the probability of error in this region decreases. Conversely, a radial protection level that contains the same probability content as a linear protection level will always be a little larger. Radial protection levels are not defined for aircraft precision approach, simply because an airport is a world full of rectangles, right angles, straight lines and, geometrically speaking, planes. However, the radial protection level concept may be of significant value in other airborne operations, particularly in relative navigation applications such as uninhabited aerial systems (UASs) or spacecraft flying in close formation. Building on the characterization of position domain error as multivariate normal with nonzero mean, the following section defines a new class of radial protection levels for use in such applications.
Use of quadratic form for K-dimensional radial protection levels A statistical quadratic form can be used to implement a radial integration limit for probability content under a multivariate normal distribution. As will be seen, this is a
major benefit from the assumption that position domain error is normally distributed. Ruben derived an infinite linear combination representation for the probability density function of a quadratic form for a finite number of correlated normally distributed RVs (Ruben 1962, 1963). He also provided estimates of error and a proof of convergence. This derivation is applicable to the present case, since the position domain RVs x, y, z and t are correlated in general by transformation through S, even if the underlying range domain RVs u1 ; u2 ; . . .; uM were assumed to be uncorrelated (here they are not). What this means practically is that the probability content of a K-dimensional fixed ellipsoid R* under a multivariate normal distribution can be evaluated numerically to an arbitrary level of accuracy, regardless of the ellipsoid’s size, orientation and offset from the distribution’s center. A sphere, circle and line—already the basis of existing GPS error characterizations—are all special cases of the ellipsoid R*. More importantly, the ellipsoidal integration volume enables the use of new, more flexible error characterizations for evaluating protection levels in GPS. The terms of Ruben’s infinite series are Chi-square distribution functions, each with a scale parameter of arbitrary size. This derivation was later reformulated, incorporating Ruben’s work as a special case along with several similar expansions and adding estimates of error and proof of convergence for the reformulation (Kotz et al. 1967). Finally, an efficient algorithm was derived for the central Chi-square expansion in these previous works (Sheil and O’Muircheartaigh 1977). The derivation appearing below closely follows the presentation of that algorithm. For the sake of generality, the derivation is carried out for the general K-dimensional case.
Central Chi-square representation of quadratic form Suppose that a K-dimensional random vector q has a multivariate normal distribution with mean vector b and variance–covariance matrix R, which is assumed to be nondegenerate. Individual elements of q are correlated. The quadratic form ðq þ aÞT C ðq þ aÞ in effect establishes an ellipsoidal region of integration, centered at point a, and with shape and orientation determined by the symmetric positive definite, or semi-definite, square matrix C. For K = 3, defining the matrix C as the 3 9 3 identity matrix establishes a spherical region of integration. Similarly, a circular region of integration is established by defining matrix C as the 2 9 2 identity matrix. Reduced rank realizations of matrix C may also be defined if protection levels are desired for all K dimensions.
123
GPS Solut
Linear transformations are performed as follows: T
T
q b ¼ L D v; a þ b ¼ L D b
ð17Þ
so that the individual elements of vector v are uncorrelated and therefore independent, identically distributed RVs with the common standard normal distribution. Matrix L is the upper triangular matrix derived from the variance–covariance matrix R ¼ LT L. The eigenvalue decomposition is then applied to matrix L C LT to produce matrix A ¼ diagðai Þ, whose main diagonal elements are the eigenvalues ai ; and matrix D, whose columns are the eigenvectors. The probability may then be reformulated (Sheil and O’Muircheartaigh 1977): P fðq þ aÞT C ðq þ aÞ sg ¼ P fðv þ bÞT A ðv þ bÞ sg ð18Þ Here, the value of scalar s ¼ r 2 , where r is the radius of the region of integration. Let f ðm; wÞ be the probability density of a central Chi-square distributed RV having m degrees of freedom. Similarly, let Fðm; wÞ be the corresponding cumulative distribution function. The probability of the quadratic forms shown in the previous equation may be represented as follows (Ruben 1962, see Equation 1.13; Kotz et al. 1967, see Eq. 93): P fðq þ aÞT C ðq þ aÞ sg 8 1 < P c Fðm þ 2k; s=fÞ; k ¼ k¼0 : 0;
if s [ 0
ð19Þ
if s 0
Here, m is also the rank of matrix C, and f is a distribution-free constant. This series is uniformly absolutely convergent for all s [ 0 as long as 0 \ f \ 2 minðai Þ (Sheil and O’Muircheartaigh 1977). In the present implementation, the value of f ¼ 29=32 is also borrowed from this reference. Thus, the series is a mixture representation such that all coefficients ck are greater than zero, and their sum is unity over the interval from k = 0 to infinity. Because the v2 functions are monotonically decreasing, truncation error after L terms is bounded according to the following expression: ! 1 L X X ck Fðm þ 2k; w=fÞ\ 1 ck Fðm þ 2L; w=fÞ k¼ Lþ1
k¼ 0
ð20Þ The computation of the expansion of the series coefficients and v2 distribution functions is available to the reader (Sheil and O’Muircheartaigh 1977). A Matlab implementation of this algorithm by Alan Genz is available at http://www.math.wsu.edu/faculty/genz.
123
Using noise and conspiring bias to set radial protection levels Multidimensional composite PLs may now be constructed. These are described below in two, three and four dimensions. RPLc The worst-case combination of range domain bias and noise is used to form a two-dimensional composite radial protection level (RPLc) for a given probability Ppl. The range domain bias bound for the ith satellite, Bmax ; i , which always has a positive value, is used. The scalar value RPLc is computed using the following steps: 1.
2.
A 2-dimensional position domain variance–covariance matrix Rxy is constructed from the first two columns and rows of matrix R. An M-dimensional range domain bias vector Bmax is created with elements Bmax ; i and by assigning a positive or negative sign to each element. Initially, all of the signs of Bmax are positive. This signed vector is mapped into the position domain to obtain the bias vector b in the x–y plane, as follows: b ¼ S Bmax
ð21Þ
Although the largest magnitude of b with respect to any of the position domain axes may be easily obtained by assigning the same signs to Bmax as those found in the corresponding row of S (effectively taking the absolute value of each), the value of b that produces the largest sum of bias and noise between the axes may be found only through a search that systematically varies the signs of the elements of vector b while maximizing r. If the elements of Bmax are all the same, the orientations of the bias cloud, i.e. the scatter plot of all possible bias combinations, and noise ellipse are correlated, and the search may possibly be shortened by adding logical conditions. Otherwise, options for limiting the search while obtaining the minimum RPLc that bounds the desired probability are few. 3. The v2 series representation previously described is evaluated using quantities b, Rxy, and Ppl from the preceding steps, along with dimension K = 2, ellipsoid shape matrix C equal to a (2 9 2) identity matrix, maximum error in probability c specified by the user (must be smaller than 1 - Ppl), and initial protection level radius r = 1 m. This step yields the probability Pc for the given r value. 4. Step 3 is iterated while varying r. First, r is changed either in the increment of ?0.1 if Pc \ Ppl or in the
GPS Solut
5.
6.
increment of -0.1 if Pc [ Ppl, until the difference Pc - Ppl changes sign. Then, this iteration is repeated using successively smaller increments of r = (0.1) j, j = 3, 5, 7, 9 (Hsu 2006). After all iterations are complete, the resulting value of r obtained for the given Ppl is stored, along with the signs of the elements of Bmax that produced it. The sign of one element of Bmax is reversed, and steps 2 through 4 are repeated. The pattern of sign changes is such that all possible combinations of signs of Bmax elements are obtained. If the current value of r is larger than the previously stored maximum value of r, the maximum value is replaced with the current value, along with the sign combination that produced it. When step 5 is completed for all sign combinations of the elements of Bmax , protection level RPLc is set to the largest value observed for r during step 5.
It is clear that this procedure obtains the minimum RPLc value that bounds the combination of conspiring range domain biases that, when mapped into the horizontal plane, creates a worst-case offset for any position domain noise ellipsoid. Thus, the RPLc value is the lower error bound for all simpler but more conservative approximations, because it takes no simplifying shortcuts. The RPLc value is exact in the sense that it corresponds to probability value Ppl with arbitrarily small probability error tolerance ±c, both of which are determined by the user. At the same time, the RPLc technique produces an upper bound on radial horizontal error that is consistent with assumptions made regarding the distribution of noise vector u and the validity of the bias bound Bmax . Because of the nested loops implied in steps 3 through 5 above, the computational inefficiency of the procedure may make real-time implementation a challenge. The value of this technique is primarily as a baseline to evaluate the accuracy of more efficient approximations.
A composite ellipsoidal protection level may likewise be computed using a slightly different C. For example, to double the safety margin in the vertical dimension while leaving the horizontal dimensions unchanged, one could use the following matrix: 3 2 1 0 0 0 60 1 0 07 7 ð23Þ C¼6 40 0 2 05 0 0 0 0 This technique may also be used to vary the shape of the EPL in along-track or cross-track dimensions, or both, depending on the application. If precise time synchronization is essential, the fourth dimension of C might be used as well. Precise aerial mapping by UASs flying in close formation over populated areas is one example application that might benefit from such a reexamination of protection levels beyond VPL and HPL.
Calculating VPLc and HPLc The procedures for obtaining the one-dimensional composite expressions VPLc and HPLc are similar to those for the radial protection levels, with steps for computing VPLc included below as an example. Once the semimajor axis length and angle are found for the horizontal plane, these same steps yield HPLc. Since both bias and noise errors map directly either to the vertical axis or to the rotated horizontal axis, the determination of biases requires no search. As mentioned previously, the probability integral is performed in one dimension. The steps for VPLc are as follows: 1.
SPLc and EPLc The quadratic form described above may be adapted to compute a composite spherical protection level simply by setting the dimension to M = 3, using a (3 9 3) identity matrix for C, and using the first three rows and columns of matrices R and S. Alternately, one can employ the original dimensions M = 4, leave the other matrices as they are, and use a reduced rank version of C. This positive semi-definite matrix, obviously still symmetric, is shown as follows: 2 3 1 0 0 0 60 1 0 07 7 C¼6 ð22Þ 40 0 1 05 0 0 0 0
Total system sigma values for each satellite are placed in an (M 9 1) column vector R and then mapped onto the vertical axis by multiplying R by the third row of the matrix S, which corresponds to the vertical dimension. No variance–covariance matrix is needed. The expression for sigma is as follows: rz ¼ S3 ; 1:M R
2.
ð24Þ
The M-dimensional range domain bias vector Bmax is created from elements Bmax ; i as in the RPLc procedure, except that the signs of Bmax are equated to the corresponding signs in the third row of S, which is for the vertical dimension z. As before, this has the effect of taking the absolute value because positive signs are multiplied by positive signs and negative signs by negative signs. The column vector Bmax is mapped onto the vertical axis by multiplying by the third row of S, as in the previous step, as follows:
123
GPS Solut
bz ¼ S3 ; 1:M Bmax 3.
4.
ð25Þ
Next, the inverse function of the cumulative distribution function of the standard normal distribution is used with input arguments rz , bz and Ppl as in Eq. (16). This involves use of a preexisting numerical routine such as the Matlab functions norminv or erf. Finally, the integration limit q obtained in the previous step, as defined in Eq. (16), is assigned as the value of VPLc.
Summary A statistical characterization of composite horizontal protection levels VPLc and HPLc under the H0 hypothesis is provided, consistent with three key assumptions: 1.
2.
3.
A certain amount of range domain bias is present in user measurements in normal GPS operation, without necessarily rejecting the H0 hypothesis. Range domain error for the ith satellite is modeled as a Gaussian RV which has non-zero mean li and variance r2i and is correlated with the Gaussian RVs representing the other satellites, in a succession of probability distributions that are piecewise continuous over time. The bias represented by the mean is constant within a subdomain, but may vary between subdomains. This is also true of the variance. The bias is estimated by B^i ðtÞ and is bounded by parameter Bmax; i . Bias bounds, but not necessarily bias values, always conspire between satellites to produce the worst-case position domain value for PL computation. The position domain error is modeled as multivariate normal with non-zero mean.
Based on these assumptions, composite protection levels VPLc and HPLc are implemented as univariate normal distributions with non-zero means. A simple method is presented by which exact values—that is, values accurate to a user-defined error tolerance and consistent with statistical assumptions—of VPLc and HPLc are obtained, and by which computationally efficient approximations may be evaluated. A statistical quadratic form under the multivariate normal distribution is then used to derive a new class of protection levels based on the probability enclosed within a radius defined in two or more dimensions. A central Chi-square representation of this quadratic form is also presented and is incorporated into a six-step computational procedure for two-dimensional composite radial protection level RPLc. This procedure is extended to the
123
three-dimensional composite SPLc and EPLc and can be used for the fourth dimension of time as well. It is hoped that placing these equations and procedures into a single document will facilitate further exploration of composite PLs. Computationally efficient approximations to exact values VPLc and HPLc are presented along with a comparison of effectiveness that uses a Monte-Carlo simulation in Bruckner et al. (2010b). The approximations are also applied and compared to existing LAAS protection level methods using flight data recorded by an Ohio University test aircraft in Bruckner et al. (2010a). Those two publications elaborate on the assumptions presented in this paper and underscore both the novelty and effectiveness of this approach. Acknowledgments This research was supported, in part, by the Federal Aviation Administration under Aviation Cooperative Research Agreement 98-G-002.
References Bruckner D, van Graas F, Skidmore T (2010a) Algorithm and flight test results to exchange code noise and multipath for biases in dual frequency differential GPS for precision approach. Navigation 57(3):1–17 Bruckner D, van Graas F, Skidmore T (2010b) Approximations to composite GPS protection levels for aircraft precision approach and landing. GPS Solutions (in press) Chatfield C (2004) The analysis of time series. Chapman and Hall/ CRC, New York Hsu D (2006) Spatial error analysis: a unified application-oriented treatment. Wiley-IEEE Press, London Kotz S, Johnson N, Boyd D (1967) Series representations of quadratic forms in normal variables, I. central case. Ann Math Stat 38(3):823–837 Leick A (2004) GPS satellite surveying, 3rd edn. Wiley, New York Misra P, Enge P (2004) Global positioning system: signals, measurements and performance. Ganga-Jamuna Press, Lincoln Murphy T, Geren P, Pankaskie T (2007) GPS antenna group delay variation induced errors in GNSS based precision approach and landing systems. In: Proceedings of the institute of navigation GPS 2974–2989, September 2007 Rife J, Pullen S (2005a) Impact of measurement biases on availability for category III LAAS. In: Proceedings of the institute of navigation Annual Mtg 759–773, June 2005 Rife J, Pullen S (2005b) Impact of measurement biases on availability for category III LAAS. Navigation 50(4):215–228 Rife J, Pullen S, Enge P, Pervan B (2006) Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Trans Aerospace Electron Syst 42(4):1386–1395 RTCA (2004) Minimum Aviation Systems Performance Standards (MASPS) for the Local Area Augmentation System (LAAS), RTCA DO-245A RTCA (2006) Minimum Operational Performance Standards (MOPS) for GPS Wide Area Augmentation System (WAAS) Airborne Equipment, RTCA DO-229D RTCA (2007) Minimum Operational Performance Standards (MOPS) for GPS LAAS Airborne Equipment, RTCA DO-253B
GPS Solut Ruben H (1962) Probability content of regions under spherical normal distributions, IV: the distribution of homogenous and non-homogenous quadratic functions of normal variables. Ann Math Stat 33(2):542–570 Ruben H (1963) A new result on the distribution of quadratic forms. Ann Math Stat 34(4):1582–1584 Sheil J, O’Muircheartaigh I (1977) AS 106: the distribution of nonnegative quadratic forms in normal variables. Appl Stat 26(1):92–98 Shively C, Hsiao T (2004) Availability enhancements for cat IIIb LAAS. Navigation 51(1):45–47 Skidmore T, van Graas F, Liu F (1996) Ohio University LAAS monitoring: design and performance. Proceedings of the institute of navigation GPS 1615–1621, September 1996 Suddapalli R (2004) Aircraft position integrity for differential satellite-based navigation in the presence of both bias and noise errors. M.S.E.E. Thesis, Ohio University Therrien C (1992) Discrete random signals and statistical signal processing. Prentice-Hall, Englewood Cliffs van Graas F, Krishnan V, Suddapalli R, Skidmore T (2004) Conspiring biases in the local area augmentation system. In: Proceedings of the institute of navigation Annual Mtg 300–307, June 2004 Wackerly D, Mendenhall W, Scheaffer R (2002) Mathematical statistics with applications. Duxbury, Pacific Grove
Author Biographies Dr. Dean Bruckner is the Assistant Director (Technical) of the Avionics Engineering Center at Ohio University in Athens, Ohio. He received his Ph.D. degree in Electrical Engineering from Ohio University in 2010 with a dissertation that included flight test results for dual-frequency GPS precision approach and landing. During a career in the United States Coast Guard, he earned a B.S.E.E. from the Coast Guard Academy and an M.S.E.E. from the Naval Postgraduate School, both with honors and with specialization in electronic navigation. He served as Commanding Officer of the Coast Guard Loran-C transmitter station near Istanbul, Turkey, and was the nationwide electronics support manager for Coast Guard ships, boats, aids to navigation and shore facilities. He has recently provided navigation and engineering expertise to the F-35 Joint Strike Fighter navigation sensor test program.
Dr. Frank van Graas is the Fritz J. and Dolores H. Russ Professor of Electrical Engineering at Ohio University. He earned his Ph.D. from Ohio University in 1988. He holds a BSEE and MSEE degrees with specialization in avionics from Delft University of Technology, The Netherlands. His extensive involvement in satellite navigation included conducting the first real-time GPS attitude and heading flight experiment on a DC-3 in 1991, and the first kinematic dual-frequency GPS autoland flight tests using NASA Langley Research Center’s Boeing 737 in 1993. He received the Johannes Kepler Award for ‘‘sustained and significant contributions to satellite navigation’’ from the Satellite Division of the Institute of Navigation in 1996, and served as President of the Institute of Navigation from 1998 to 1999. At Ohio University’s Avionics Engineering Center, his roles include Principal Investigator for the FAAsponsored Global Positioning System (GPS) research grant to investigate the Local Area Augmentation System for Aircraft Precision Approach & Landing and Aircraft Surface Movement Guidance. He has lectured extensively throughout the United States, Canada, Europe and Russia for such organizations as the International Air Transport Association, the NATO Advisory Group on Aerospace Research and Development (AGARD), Navtech Seminars, Inc., and the International Society for Optical Engineering (SPIE). Dr. Trent Skidmore is a Senior Research Engineer and Adjunct Assistant Professor in the Avionics Engineering Center at Ohio University. He is involved in a variety of projects involving applications of GPS, including the use of differential GPS for aircraft precision approach applications. Recent projects include the Department of Defense Joint Precision Approach and Landing System (JPALS), the Automated Aerial Refueling (AAR) program, as well as the Federal Aviation Administration Local Area Augmentation System (LAAS). Trent received his B.Sc. in 1986 (with Honors) and M.Sc. in 1988 from Michigan Technological University and his Ph.D. from Ohio University in 1991, all in electrical engineering.
123
ORIGINAL ARTICLE
Statistical characterization of composite protection levels for GPS Dean Bruckner • Frank van Graas Trent Skidmore
•
Received: 14 October 2009 / Accepted: 20 September 2010 Springer-Verlag 2010
Abstract A statistical characterization is presented of Global Positioning System (GPS) user range error as a normally distributed random variable with non-zero mean over the length of the aircraft precision approach operation, correlated from one GPS measurement epoch to another and from one satellite to another. This leads directly to modeling GPS error in the position domain as multivariate normal with non-zero mean. Based on this model, a vertical composite protection level VPLc and a horizontal composite protection level HPLc are each calculated as scalar values from a univariate normal distribution displaced from the origin by the worst-case position domain bias combination possible, given the maximum possible individual satellite bias magnitudes and the satellite geometry. A method is then presented by which exact values—that is, values accurate to a user-defined error tolerance and consistent with statistical assumptions—of VPLc and HPLc are obtained, and by which computationally efficient approximations may be evaluated. A statistical quadratic form under the multivariate normal distribution is used to derive a new class of protection levels based on the probability enclosed within a radius defined in two or more dimensions. A central Chi-square representation of this quadratic form is also presented and is incorporated into a six-step computational procedure for the two-dimensional composite radial protection level RPLc. This procedure is extended to the spherical protection level (SPLc) and the ellipsoidal protection level (EPLc).
D. Bruckner (&) F. van Graas T. Skidmore Ohio University, Athens, OH, USA e-mail: [email protected] URL: www.ohio.edu/avionics
Keywords VPL HPL Composite protection level Radial protection level Ellipsoidal protection level Conspiring bias
Introduction Setting a protection level (PL) is a widely accepted integrity methodology for navigation systems, particularly in safety-intensive aviation applications. A PL is commonly defined as the maximum navigation system error (NSE) from the true position that can occur in a specified direction such as the vertical, ‘‘guaranteed’’ to a sufficiently high level of probability. This section introduces competing approaches to setting PLs and provides context for the statistical characterization of the present effort. Protection levels in aviation safety Integrity specifications in the GPS-based Local Area Augmentation System (LAAS) and Wide Area Augmentation System (WAAS) call for calculating vertical, horizontal and lateral (cross-track) PLs, which are abbreviated as VPL, HPL and LPL, respectively (RTCA 2004, 2006, 2007). GPS Receiver Autonomous Integrity Monitoring (RAIM) algorithms likewise compute a VPL or HPL. In these types of applications, the PL is continually compared to an appropriate Alert Limit (AL) throughout the operation. If the PL exceeds the AL, an alarm is annunciated to the user. This typically results in the navigation system being declared unusable for the intended purpose until the PL again drops below the AL. In the RAIM context, system design emphasis is usually placed on detecting a satellite malfunction, typically marked by a growing pseudorange bias, before that bias becomes large enough to pose a
123
GPS Solut
serious problem. In the present effort, however, it is assumed that a certain level of pseudorange bias exists in normal GPS operation. This point has been illustrated in recent years in various studies of residual group delay bias in GPS antennas and antenna arrays (Murphy et al. 2007). The composite PL is an attempt to treat this bias explicitly to improve overall differential GPS performance. Integrity for differential GPS systems such as WAAS and LAAS has traditionally been assured by over-bounding both the bias-like errors, i.e. persistent and slowly varying but not necessarily growing non-zero mean errors, and noise-like errors seen in GPS measurements by using a single closed-form probability distribution. Typically, this has been the normal distribution. A different approach has also been advanced, one that explicitly assigns probability distributions separately to noise and to biases. The PL that results from this separate treatment is a composite of noise and bias bounds. Previously defined composite vertical protection level definitions, such as VPL–C (van Graas et al. 2004) and VPLbias (Rife and Pullen 2005b), are prominent examples. In the present effort, a statistical method for computing a new vertical composite protection (VPLc) and horizontal composite protection level (HPLc) is defined. A new class of composite protection levels is also developed from a quadratic statistical form. This class includes the twodimensional radial protection level (RPLc), as well as the three-dimensional spherical protection level (SPLc) and ellipsoidal protection level (EPLc). Previously defined composite PLs mentioned above are seen to be approximations to the PLs presented here. As others have done before, this effort examines PLs only under the H0, or fault-free, hypothesis (Rife and Pullen 2005a; Shively and Hsiao 2004). This is appropriate for the analysis of a generic PL formulation not closely tied to the configuration of any particular system such as LAAS or WAAS. It is also appropriate in order to permit a focus on improving system accuracy performance through error reduction. The H0 hypothesis governing the present analysis is simply normal GPS system operation, with the additional proviso that some level of pseudorange bias exists in normal GPS operation, as asserted above. The mere presence of such bias, therefore, does not necessarily lead to the rejection of this H0 hypothesis and the subsequent declaration of a system fault. The present effort is motivated largely by repeated examination of slowly varying biases in flight test data recorded in differential GPS-based precision approach and landing operations over a 15-year period (van Graas et al. 2004, Skidmore et al. 1996). The presence of persistent position domain biases that remain essentially constant over an individual aircraft precision approach but vary between approaches, as well as the shortcomings in
123
traditional PL techniques in adequately representing these biases, justify exploring a new technique. Such a technique must be able to separate differential GPS spatial and temporal errors cleanly into Gaussian and non-Gaussian components that may each be bounded. The authors present such techniques elsewhere, along with flight test results that demonstrate significant decreases in the size of composite PLs when compared to the traditional PLs used in LAAS (Bruckner et al. 2010a). Conceptually, decorrelation between reference station and aircraft observations due to localized tropospheric and ionospheric effects, if not removed by ionospheric-free processing, can be accounted for by inflating the bias term bound. Yet another approach could be taken using Extreme Value Theory (EVT). This method provides advantages in error model analysis for unknown distributions. However, since an EVT-produced probability distribution is of an unknown type, it does not seem useful for setting PLs in the high-integrity precision approach and landing context. For this reason, the EVT method is not considered or further explained here. A case study illustration We present a case study illustration of a general statistical model previously developed by others (Ruben 1962; Ruben 1963; Kotz et al. 1967; Sheil and O’Muircheartaigh 1977). The model is applied here to the linearized solution method of the GPS equations in a new way. It is known that the modeling of GPS pseudorange error as unbiased, uncorrelated, Gaussian random variables (RVs) has great practical value. Here, that model is extended by redefining the error on each pseudorange as a Gaussian RV correlated with the Gaussian RVs that describe errors on the other pseudoranges. This is intended to model both the observed correlation in biases and the observed absence of correlation in noise between satellites. These Gaussian RVs are then employed in a quadratic form to define composite protection levels that treat noise and bias separately for the differential GPS-based precision approach and landing operation. The advantage of this method, which is not well known in the navigation community, is readily apparent in the resulting definition of protection levels that are less conservative than those constructed for the unbiased model. This in turn benefits system availability and continuity. The novelty of the present effort is threefold: the use of a quadratic form to define composite PLs that treat noise and bias separately in the context of differential GPS for precision approach and landing; the expressions provided for propagating noise using non-traditional weighting matrices; and the definition of range domain error for a particular satellite as a piecewise fit over time
GPS Solut
of a succession of Nðli ; r2i Þ-distributed RVs, which are correlated with the RVs representing the other satellites. It is hoped that placing these equations for composite PLs into a single document will make this approach more accessible to the navigation community, particularly when evaluating the effectiveness of computationally efficient approximations to these more exact PL expressions. Due to lack of space, such approximations are presented elsewhere by the authors in a follow-on paper (Bruckner et al. 2010b).
Designation of random variables The random variables are briefly defined in this section. A few well-known expressions such as the GPS measurement model are briefly restated to highlight novel points in the new technique.
It is assumed that a differential GPS system is operating normally. A linearized GPS measurement model is also assumed. Although this derivation is not limited to LAAS applications, it begins with the same GPS measurement model as defined for that system (RTCA 2004). The east– north–up locally level coordinate frame assumed here is a slight modification to that definition, in which the LAAS horizontal coordinate axes are either aligned with the runway or are arbitrary. The GPS measurement model is stated briefly as follows: ð1Þ
where Dy is the (M 9 1) vector of ranging measurements, minus expected ranges based on known satellite positions, assumed user position and user clock offset; H, the (M 9 4) observation matrix, in which the first three elements of each row form a unit vector from the user to the satellite in the locally level coordinate frame, and the fourth element is 1; Dx, the (4 9 1) vector from the user’s assumed position to the user’s true position in the locally level frame, including time as a fourth dimension; and e, the (M 9 1) vector of ranging measurement errors, to be defined further below. The difference between assumed user position and true user position may be estimated as follows, through the use of multiple linear regression techniques that employ weighted least squares: D^ x ¼ S Dy where
Weighting matrix W 2 2 r1 0 6 0 r22 6 W1 ¼ 6 . .. . . 4 .. . . 0 0
ð3Þ is defined as follows: 3 0 0 7 7 .. 7 . 5
ð2Þ
ð4Þ
r2M
where r2i is the variance of the unbiased and uncorrelated Gaussian range domain error of the ith satellite measurements that results from ground receiver noise and multipath, airborne receiver noise and multipath, and residual delays from ionospheric and tropospheric gradients between the reference station and the user. Matrix S is the projection matrix from the range domain to the position domain. The first modification to the GPS measurement model is to characterize user range domain error vector e using M normally distributed, discrete time random processes: ui;k Nðli ; r2i Þ
The GPS measurement model
Dy ¼ H Dx þ e
S ¼ ðHT W HÞ1 HT W
ð5Þ
where ui;k is the ith random sequence that represents a continuous time random variable, uniformly sampled at period DT seconds; li , the mean of the ith random sequence; r2i , the variance of the ith random sequence; i = 1, 2,…, M is the satellite index; and k = …, -2, -1, 0, 1, 2 … is the sample index and time variable. These processes are assumed to be discrete time Gaussian stationary with constant mean and variance over the duration of an aircraft precision approach and landing operation. Accordingly, the derivations presented in this section are equally valid for each GPS measurement epoch. Therefore, only the current epoch, where sample subscript k = 0, is considered. The sample subscript k is suppressed. Range domain random variables Now let M-dimensional real random vector u consist of correlated RVs u1 ; u2 ; . . .; uM , and let these RVs characterize range domain error for the duration of the precision approach and landing operation. It follows from the definition of ui;k above that each of these RVs is normally distributed with parameters mean li and variance r2i . These parameters are themselves constant over the operation, but they may differ from one satellite to another during this time interval. The assumption of correlation models underlying correlation of pseudorange biases between satellites, such as can occur due to group delay when the observed angle between two satellites, as apparent to a user’s antenna, is small. The pseudorange noise of each satellite is still assumed to be uncorrelated with that of other satellites. In the formation of PLs, the assumption of
123
GPS Solut
correlated biases between satellites is taken a step further by intentionally correlating them perfectly to produce the maximum possible position domain bias to make up for not knowing the probability distribution of the biases. When compared to the LAAS error model presented above, this definition is seen to be a departure from conventional treatment of range domain user measurement error, in which u1 ; u2 ; . . .; uM are considered uncorrelated, zero mean, identically distributed normal variables. The conventional model has been widely accepted not because it fits GPS behavior perfectly, but rather because it has worked without imposing a great deal of complexity for minimal added performance gain (Misra and Enge 2004). However, new GPS user requirements are pressing the system to its limits and, as noise filtering is increasingly effective, the presence of biases is receiving greater recognition. Therefore, long-standing assumptions are being reexamined by researchers in this field. This paper constitutes one such reexamination. The model proposed here is motivated by the desire to make practical use of the phenomenon of biases in user GPS measurements that persist essentially unchanged over the five- to ten-minute duration of a precision approach, but may change noticeably by the beginning of the next approach. Indeed, GPS range domain measurements are known to exhibit both correlation over time and correlation between simultaneous measurements from different satellites (Misra and Enge 2004). Three alternate range domain error models were considered besides the one presented here. In the first case, treating range domain error as a linear process driven by white noise has the virtues of constant mean and variance, thus ensuring both wide-sense stationarity and normally distributed error in both the range domain and position domain. Furthermore, realizations of this process may legitimately take on successive values distant from the mean for an extended period of time, relative to the short duration of a precision approach. This behavior can therefore model a slowly changing bias, while remaining assured of convergence (Chatfield 2004). However, the ability to translate the current deviation from the constant process mean, i.e. the instantaneous bias, into the position domain, and there to make statistical sense of it, remains out of reach. A second model was described in (Skidmore et al. 1996), in which the convolution of uniform and Gaussian distributions represents the addition of bias-like and noise-like range domain errors, respectively. This description is straightforward, but does not yield a normal distribution in the position domain. A third model, defined in (Rife et al. 2006), accounts for range domain biases by convolving an identical pair of distributions, e.g. normal distributions when applied to GPS PLs, with each other, where one member of the pair is offset from the origin by
123
the positive magnitude of the range domain bias, and the other by the negative magnitude. The present effort can be viewed as a limiting case of the paired distribution approach. If the identical distributions represent a probability assignment of 50% to the positive maximum of the bias and 50% to the negative maximum, accomplished through the sifting property of convolution with the delta function, the approach in this paper is to assign 100% of the probability to one bound, whichever one has a greater magnitude, and zero to the other. If the two bounds have the same magnitude, then which one is selected does not matter. In any case, the PL’s symmetry about zero ensures a worst-case bound. In summary, the user range error u is defined here as normally distributed, but over the limited time horizon of a typical aircraft precision approach and landing operation. Mathematically, this model represents a piecewise fit of a succession of Nðli ; r2i Þ-distributed RVs over time. The parameters li and r2i for the ith satellite may change from one subdomain, or time segment, to another. However, within the kth subdomain, li and r2i are assumed to be constant. The quantities li and r2i may also be described as piecewise continuous functions, each of which consists of a series of steps. These steps occur on transition from one subdomain to the next. Given the frequent jumps in protection levels seen during a typical precision approach as currently implemented, discontinuities of bounded magnitude in li and r2i at subdomain boundaries do not appear to present a problem to aircraft operations. The value of bias li for a given satellite is generally unknown, but is contained within the symmetric bounds that are described above and that may decrease over time in some processing techniques. Such techniques may also estimate the bias value itself with increasing accuracy over time. Range domain errors are assumed to be correlated between satellites, also as described above. Correlation across time does not enter directly into the PL calculation performed at each GPS measurement epoch and is not considered further. Since correlated variables are already assumed in the PL derivation that follows, the assumption of correlation here does not invalidate it. One drawback with this assumption is that the position domain Dilution of Precision (DOP) values are now correlated and their usefulness as individual scalar parameters declines. Position domain random variables The position domain error, defined as the difference between estimated position offset D^x and true position offset Dx from an assumed position, may be characterized as a four-dimensional real random vector q with
GPS Solut
component RVs x, y, z and t. These are east, north and up in a locally level coordinate system and user receiver clock offset from GPS time. The random vector q has a quadvariate normal distribution, i.e. multivariate normal with four dimensions, inherited from multivariate normal vector u via linear transformation, as defined 2 3 2 3 u1 x 6 u2 7 6y7 6 7 7 q ¼ S u or 6 ð6Þ . 7 4z5 ¼ S6 4 .. 5 t uM where projection matrix S, defined in (3), is a deterministic linear transformation. Then RV x, which is the position domain error in the east direction, is written as x¼
M X
s 1 ; i ui
ð7Þ
i¼1
where s1 ; i is the element in the first row and ith column of matrix S. The RVs y, z and t are defined similarly. It is assumed that the elements of matrix S are known with certainty. Therefore, any actual geometry errors due to ephemeris inaccuracies, user position error and the like may be projected onto the line of sight between the assumed user position and the calculated satellite position, and then combined with the pseudorange bias error. Thus, the elements of matrix S may be treated deterministically as constants instead of as additional random parameters. In LAAS development, a great deal of effort was devoted to ensuring that large ephemeris errors would not result in differences in the projection of errors between ground and air. Whenever this algorithm is implemented, the reference station should be assigned to monitor this error source. Mean and variance of vector q The mean vector l and variance–covariance matrix R of vector q may be derived using full matrix notation and beginning directly with RVs x, y, z and t (Misra and Enge 2004; Therrien 1992). They may also be computed using the simple linear transformations of the GPS measurement model and the well-known theorem on the numerical characteristics of a linear combination of RVs (Wackerly et al. 2002, see Theorem 5.12). Covariance propagation from the range domain to the position domain is also treated extensively in a number of texts (Leick 2004) and is not repeated here. Modeling noise propagation: unweighted and weighted least squares If noise sigma values are not considered identical for all the satellites, then the propagation of non-identically distributed noise must be modeled when constructing the
variance–covariance matrix, which in turn is used to set the protection level. Such is the case in this statistical model. This modeling happens automatically in the equation for VPL for the LAAS H0 hypothesis, in which the least squares projection matrix S is formed using weighting terms 1 r2i and multiplied by the squares of the corresponding elements of the S matrix. The equation for VPLH0 assumes no bias and is thus appropriate to illustrate the propagation of noise; this equation is as follows (RTCA 2004), vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi u M uX VPL ð8Þ ¼ K t s2 r 2 H0
ffmd
3;i i
i¼1
where Kffmd is the fault-free missed-detection sigma multiplier. If, however, the unweighted least squares method is used for position solutions rather than the weighted method, then multiplying r2i by the squares of the elements of the unweighted S matrix to create desired terms of variance–covariance matrix R no longer correctly models the propagation of non-uniform noise. In this case, the revised variance–covariance matrix is instead obtained using notation G for the unweighted projection matrix (Suddapalli 2004), R0 ¼ G Q GT
ð9Þ
where the variance–covariance matrix of range domain error u is 2 2 3 r1 0 0 6 0 r22 0 7 6 7 Q¼6 . ð10Þ .. 7 .. . . 4 .. . 5 . . 0
0
r2M
and G ¼ ðHT HÞ1 HT
ð11Þ
As described above, noise is assumed to be uncorrelated from one satellite to another. In this case, the corrected VPL equation for the bias-free case would use the (3, 3) element of matrix R0 as pffiffiffiffiffiffiffiffiffiffiffi ð12Þ VPLH0 ¼ Kffmd R0 3 ; 3 If position solutions are obtained using a weighting scheme other than either unweighted least squares or least squares weighted by 1 r2i , then an additional correction is required. For the case of least squares weighted by 1=li , for example, the inverse of the weighting matrix is 2 3 l1 0 0 6 0 l2 0 7 6 7 W1 ¼ 6 . ð13Þ .. 7 .. . . 4 .. . 5 . . 0
0
lM
123
GPS Solut
and the S matrix, as before, is defined as T
1
S ¼ ðH W HÞ HT W
ð14Þ
distribution in the position domain with mean b and standard deviation r, defined as Zq
n o 1 pffiffiffiffiffiffi exp ðr bÞ2 =ð2r2 Þ dr r 2p
Then, the revised covariance matrix that correctly models the propagation of noise through the least squares method weighted by 1=li is
p¼
R00 ¼ S Q ST
In this PL context, position domain bias b is set to the bias bound, as in Fig. 1. The original VPLc value, computed as b ? Kr and shown in Fig. 1, is an artifact of the protection level for the zero-mean distribution from which it was derived. Because the PL must be symmetric about the origin, the integral of the probability density function with integration limits (b ? Kr) and b ? Kr contains a significantly larger probability content than for the original integration limits -Kr and ?Kr. Note that the interval [-Kr, Kr] for the zero-mean, i.e. unbiased case is the same size as the interval [b - Kr, b ? Kr] for the biased case, even though it has been translated from the origin. This translated interval, however, cannot be used as a protection level due to its asymmetry about the origin. If instead the value of q can be found numerically by varying the integration limits until the desired value of p is achieved, then an exact probability may be enclosed by a new VPLc, to the specified error tolerance and consistent with the assumptions. The same approach is applied to define HPLc. The quantity HPLc, like the original HPL, is defined under a univariate normal distribution as the projection of horizontal errors onto the major axis of the error ellipse in the horizontal plane. Because the semimajor axis length of an ellipse is its largest radius, the HPL and HPLc are worstcase measures of generic lateral (cross-track) error and are used when the orientation of a runway is unknown. The HPL is not properly a radial measure, in which a specified probability is enclosed within a given radius. In fact, an HPL interpreted as a radial bound will always contain a smaller probability than the same HPL value interpreted as a linear bound. Figure 2 illustrates the errors
ð15Þ
To avoid a weighting matrix that is not positive definite, minimum and maximum limits are placed on the values of li . Use of these limits prevents algorithm instability while not noticeably affecting results. Since noise and bias are treated separately in forming the composite PL, the noise propagation techniques described in this section remain valid beyond the bias-free case shown here. Bias uncertainty and formation of the protection level An estimator B^i ðtÞ may be developed for the bias li in pseudorange measurements from the ith satellite at time t. The notation B, more familiar than l in GPS contexts, is used for range domain bias hereafter. Although the estimator B^i ðtÞ is not specified in this paper, it is assumed that whatever sampling distribution is selected constitutes a 100% confidence interval for B^i ðtÞ, centered at zero with upper and lower bounds Bmax ; i . One approach assumes a uniform sampling distribution (Skidmore et al. 1996), while another would assume a sinusoidal distribution more heavily weighted toward the upper and lower bounds. Because the distribution is here assumed to be unknown except for the stipulation that it is contained within bounds Bmax ; i , the non-Gaussian portion of the composite PL is formed using a worst-case vector addition of all the Bmax ; i . The bias bound and the Gaussian noise bound are added to form the PL. The advantage of treating the noise and bias separately is primarily that the noise bound now contains only thermal noise and is thus much smaller than when other terms are included. Furthermore, the decrease in the noise bound typically exceeds the increase in the bias bound, resulting in a net decrease to the size of the PL.
q
ð16Þ
Statistical method for exact VPLc and HPLc When viewed within the statistical framework described above, the ‘‘original’’ composite VPLc expressions noted in the literature are actually conservative approximations for a true value that can be derived statistically. This true value q is the magnitude of the upper and lower integration limits, which are symmetric about zero, of the integral of the probability density function of a normal distribution having a non-zero mean and a given probability p. Consider such a
123
Fig. 1 Probability density function for a biased univariate normal distribution in the position domain, with original and new VPLc shown. Mean value b is set to the bias bound
GPS Solut
Fig. 2 Errors enclosed within ‘‘linear’’ HPL but not enclosed within ‘‘radial’’ HPL of the same magnitude
projected onto the ‘‘linear’’ HPL but not enclosed within the ‘‘radial’’ HPL, using data from a Monte-Carlo simulation created by the authors. These points between the straight line and circle, plotted with a larger symbol than the other points, are shown for the [B - 3r, B ? 3r] case for illustration purposes. As the sigma multiplier K increases, the probability of error in this region decreases. Conversely, a radial protection level that contains the same probability content as a linear protection level will always be a little larger. Radial protection levels are not defined for aircraft precision approach, simply because an airport is a world full of rectangles, right angles, straight lines and, geometrically speaking, planes. However, the radial protection level concept may be of significant value in other airborne operations, particularly in relative navigation applications such as uninhabited aerial systems (UASs) or spacecraft flying in close formation. Building on the characterization of position domain error as multivariate normal with nonzero mean, the following section defines a new class of radial protection levels for use in such applications.
Use of quadratic form for K-dimensional radial protection levels A statistical quadratic form can be used to implement a radial integration limit for probability content under a multivariate normal distribution. As will be seen, this is a
major benefit from the assumption that position domain error is normally distributed. Ruben derived an infinite linear combination representation for the probability density function of a quadratic form for a finite number of correlated normally distributed RVs (Ruben 1962, 1963). He also provided estimates of error and a proof of convergence. This derivation is applicable to the present case, since the position domain RVs x, y, z and t are correlated in general by transformation through S, even if the underlying range domain RVs u1 ; u2 ; . . .; uM were assumed to be uncorrelated (here they are not). What this means practically is that the probability content of a K-dimensional fixed ellipsoid R* under a multivariate normal distribution can be evaluated numerically to an arbitrary level of accuracy, regardless of the ellipsoid’s size, orientation and offset from the distribution’s center. A sphere, circle and line—already the basis of existing GPS error characterizations—are all special cases of the ellipsoid R*. More importantly, the ellipsoidal integration volume enables the use of new, more flexible error characterizations for evaluating protection levels in GPS. The terms of Ruben’s infinite series are Chi-square distribution functions, each with a scale parameter of arbitrary size. This derivation was later reformulated, incorporating Ruben’s work as a special case along with several similar expansions and adding estimates of error and proof of convergence for the reformulation (Kotz et al. 1967). Finally, an efficient algorithm was derived for the central Chi-square expansion in these previous works (Sheil and O’Muircheartaigh 1977). The derivation appearing below closely follows the presentation of that algorithm. For the sake of generality, the derivation is carried out for the general K-dimensional case.
Central Chi-square representation of quadratic form Suppose that a K-dimensional random vector q has a multivariate normal distribution with mean vector b and variance–covariance matrix R, which is assumed to be nondegenerate. Individual elements of q are correlated. The quadratic form ðq þ aÞT C ðq þ aÞ in effect establishes an ellipsoidal region of integration, centered at point a, and with shape and orientation determined by the symmetric positive definite, or semi-definite, square matrix C. For K = 3, defining the matrix C as the 3 9 3 identity matrix establishes a spherical region of integration. Similarly, a circular region of integration is established by defining matrix C as the 2 9 2 identity matrix. Reduced rank realizations of matrix C may also be defined if protection levels are desired for all K dimensions.
123
GPS Solut
Linear transformations are performed as follows: T
T
q b ¼ L D v; a þ b ¼ L D b
ð17Þ
so that the individual elements of vector v are uncorrelated and therefore independent, identically distributed RVs with the common standard normal distribution. Matrix L is the upper triangular matrix derived from the variance–covariance matrix R ¼ LT L. The eigenvalue decomposition is then applied to matrix L C LT to produce matrix A ¼ diagðai Þ, whose main diagonal elements are the eigenvalues ai ; and matrix D, whose columns are the eigenvectors. The probability may then be reformulated (Sheil and O’Muircheartaigh 1977): P fðq þ aÞT C ðq þ aÞ sg ¼ P fðv þ bÞT A ðv þ bÞ sg ð18Þ Here, the value of scalar s ¼ r 2 , where r is the radius of the region of integration. Let f ðm; wÞ be the probability density of a central Chi-square distributed RV having m degrees of freedom. Similarly, let Fðm; wÞ be the corresponding cumulative distribution function. The probability of the quadratic forms shown in the previous equation may be represented as follows (Ruben 1962, see Equation 1.13; Kotz et al. 1967, see Eq. 93): P fðq þ aÞT C ðq þ aÞ sg 8 1 < P c Fðm þ 2k; s=fÞ; k ¼ k¼0 : 0;
if s [ 0
ð19Þ
if s 0
Here, m is also the rank of matrix C, and f is a distribution-free constant. This series is uniformly absolutely convergent for all s [ 0 as long as 0 \ f \ 2 minðai Þ (Sheil and O’Muircheartaigh 1977). In the present implementation, the value of f ¼ 29=32 is also borrowed from this reference. Thus, the series is a mixture representation such that all coefficients ck are greater than zero, and their sum is unity over the interval from k = 0 to infinity. Because the v2 functions are monotonically decreasing, truncation error after L terms is bounded according to the following expression: ! 1 L X X ck Fðm þ 2k; w=fÞ\ 1 ck Fðm þ 2L; w=fÞ k¼ Lþ1
k¼ 0
ð20Þ The computation of the expansion of the series coefficients and v2 distribution functions is available to the reader (Sheil and O’Muircheartaigh 1977). A Matlab implementation of this algorithm by Alan Genz is available at http://www.math.wsu.edu/faculty/genz.
123
Using noise and conspiring bias to set radial protection levels Multidimensional composite PLs may now be constructed. These are described below in two, three and four dimensions. RPLc The worst-case combination of range domain bias and noise is used to form a two-dimensional composite radial protection level (RPLc) for a given probability Ppl. The range domain bias bound for the ith satellite, Bmax ; i , which always has a positive value, is used. The scalar value RPLc is computed using the following steps: 1.
2.
A 2-dimensional position domain variance–covariance matrix Rxy is constructed from the first two columns and rows of matrix R. An M-dimensional range domain bias vector Bmax is created with elements Bmax ; i and by assigning a positive or negative sign to each element. Initially, all of the signs of Bmax are positive. This signed vector is mapped into the position domain to obtain the bias vector b in the x–y plane, as follows: b ¼ S Bmax
ð21Þ
Although the largest magnitude of b with respect to any of the position domain axes may be easily obtained by assigning the same signs to Bmax as those found in the corresponding row of S (effectively taking the absolute value of each), the value of b that produces the largest sum of bias and noise between the axes may be found only through a search that systematically varies the signs of the elements of vector b while maximizing r. If the elements of Bmax are all the same, the orientations of the bias cloud, i.e. the scatter plot of all possible bias combinations, and noise ellipse are correlated, and the search may possibly be shortened by adding logical conditions. Otherwise, options for limiting the search while obtaining the minimum RPLc that bounds the desired probability are few. 3. The v2 series representation previously described is evaluated using quantities b, Rxy, and Ppl from the preceding steps, along with dimension K = 2, ellipsoid shape matrix C equal to a (2 9 2) identity matrix, maximum error in probability c specified by the user (must be smaller than 1 - Ppl), and initial protection level radius r = 1 m. This step yields the probability Pc for the given r value. 4. Step 3 is iterated while varying r. First, r is changed either in the increment of ?0.1 if Pc \ Ppl or in the
GPS Solut
5.
6.
increment of -0.1 if Pc [ Ppl, until the difference Pc - Ppl changes sign. Then, this iteration is repeated using successively smaller increments of r = (0.1) j, j = 3, 5, 7, 9 (Hsu 2006). After all iterations are complete, the resulting value of r obtained for the given Ppl is stored, along with the signs of the elements of Bmax that produced it. The sign of one element of Bmax is reversed, and steps 2 through 4 are repeated. The pattern of sign changes is such that all possible combinations of signs of Bmax elements are obtained. If the current value of r is larger than the previously stored maximum value of r, the maximum value is replaced with the current value, along with the sign combination that produced it. When step 5 is completed for all sign combinations of the elements of Bmax , protection level RPLc is set to the largest value observed for r during step 5.
It is clear that this procedure obtains the minimum RPLc value that bounds the combination of conspiring range domain biases that, when mapped into the horizontal plane, creates a worst-case offset for any position domain noise ellipsoid. Thus, the RPLc value is the lower error bound for all simpler but more conservative approximations, because it takes no simplifying shortcuts. The RPLc value is exact in the sense that it corresponds to probability value Ppl with arbitrarily small probability error tolerance ±c, both of which are determined by the user. At the same time, the RPLc technique produces an upper bound on radial horizontal error that is consistent with assumptions made regarding the distribution of noise vector u and the validity of the bias bound Bmax . Because of the nested loops implied in steps 3 through 5 above, the computational inefficiency of the procedure may make real-time implementation a challenge. The value of this technique is primarily as a baseline to evaluate the accuracy of more efficient approximations.
A composite ellipsoidal protection level may likewise be computed using a slightly different C. For example, to double the safety margin in the vertical dimension while leaving the horizontal dimensions unchanged, one could use the following matrix: 3 2 1 0 0 0 60 1 0 07 7 ð23Þ C¼6 40 0 2 05 0 0 0 0 This technique may also be used to vary the shape of the EPL in along-track or cross-track dimensions, or both, depending on the application. If precise time synchronization is essential, the fourth dimension of C might be used as well. Precise aerial mapping by UASs flying in close formation over populated areas is one example application that might benefit from such a reexamination of protection levels beyond VPL and HPL.
Calculating VPLc and HPLc The procedures for obtaining the one-dimensional composite expressions VPLc and HPLc are similar to those for the radial protection levels, with steps for computing VPLc included below as an example. Once the semimajor axis length and angle are found for the horizontal plane, these same steps yield HPLc. Since both bias and noise errors map directly either to the vertical axis or to the rotated horizontal axis, the determination of biases requires no search. As mentioned previously, the probability integral is performed in one dimension. The steps for VPLc are as follows: 1.
SPLc and EPLc The quadratic form described above may be adapted to compute a composite spherical protection level simply by setting the dimension to M = 3, using a (3 9 3) identity matrix for C, and using the first three rows and columns of matrices R and S. Alternately, one can employ the original dimensions M = 4, leave the other matrices as they are, and use a reduced rank version of C. This positive semi-definite matrix, obviously still symmetric, is shown as follows: 2 3 1 0 0 0 60 1 0 07 7 C¼6 ð22Þ 40 0 1 05 0 0 0 0
Total system sigma values for each satellite are placed in an (M 9 1) column vector R and then mapped onto the vertical axis by multiplying R by the third row of the matrix S, which corresponds to the vertical dimension. No variance–covariance matrix is needed. The expression for sigma is as follows: rz ¼ S3 ; 1:M R
2.
ð24Þ
The M-dimensional range domain bias vector Bmax is created from elements Bmax ; i as in the RPLc procedure, except that the signs of Bmax are equated to the corresponding signs in the third row of S, which is for the vertical dimension z. As before, this has the effect of taking the absolute value because positive signs are multiplied by positive signs and negative signs by negative signs. The column vector Bmax is mapped onto the vertical axis by multiplying by the third row of S, as in the previous step, as follows:
123
GPS Solut
bz ¼ S3 ; 1:M Bmax 3.
4.
ð25Þ
Next, the inverse function of the cumulative distribution function of the standard normal distribution is used with input arguments rz , bz and Ppl as in Eq. (16). This involves use of a preexisting numerical routine such as the Matlab functions norminv or erf. Finally, the integration limit q obtained in the previous step, as defined in Eq. (16), is assigned as the value of VPLc.
Summary A statistical characterization of composite horizontal protection levels VPLc and HPLc under the H0 hypothesis is provided, consistent with three key assumptions: 1.
2.
3.
A certain amount of range domain bias is present in user measurements in normal GPS operation, without necessarily rejecting the H0 hypothesis. Range domain error for the ith satellite is modeled as a Gaussian RV which has non-zero mean li and variance r2i and is correlated with the Gaussian RVs representing the other satellites, in a succession of probability distributions that are piecewise continuous over time. The bias represented by the mean is constant within a subdomain, but may vary between subdomains. This is also true of the variance. The bias is estimated by B^i ðtÞ and is bounded by parameter Bmax; i . Bias bounds, but not necessarily bias values, always conspire between satellites to produce the worst-case position domain value for PL computation. The position domain error is modeled as multivariate normal with non-zero mean.
Based on these assumptions, composite protection levels VPLc and HPLc are implemented as univariate normal distributions with non-zero means. A simple method is presented by which exact values—that is, values accurate to a user-defined error tolerance and consistent with statistical assumptions—of VPLc and HPLc are obtained, and by which computationally efficient approximations may be evaluated. A statistical quadratic form under the multivariate normal distribution is then used to derive a new class of protection levels based on the probability enclosed within a radius defined in two or more dimensions. A central Chi-square representation of this quadratic form is also presented and is incorporated into a six-step computational procedure for two-dimensional composite radial protection level RPLc. This procedure is extended to the
123
three-dimensional composite SPLc and EPLc and can be used for the fourth dimension of time as well. It is hoped that placing these equations and procedures into a single document will facilitate further exploration of composite PLs. Computationally efficient approximations to exact values VPLc and HPLc are presented along with a comparison of effectiveness that uses a Monte-Carlo simulation in Bruckner et al. (2010b). The approximations are also applied and compared to existing LAAS protection level methods using flight data recorded by an Ohio University test aircraft in Bruckner et al. (2010a). Those two publications elaborate on the assumptions presented in this paper and underscore both the novelty and effectiveness of this approach. Acknowledgments This research was supported, in part, by the Federal Aviation Administration under Aviation Cooperative Research Agreement 98-G-002.
References Bruckner D, van Graas F, Skidmore T (2010a) Algorithm and flight test results to exchange code noise and multipath for biases in dual frequency differential GPS for precision approach. Navigation 57(3):1–17 Bruckner D, van Graas F, Skidmore T (2010b) Approximations to composite GPS protection levels for aircraft precision approach and landing. GPS Solutions (in press) Chatfield C (2004) The analysis of time series. Chapman and Hall/ CRC, New York Hsu D (2006) Spatial error analysis: a unified application-oriented treatment. Wiley-IEEE Press, London Kotz S, Johnson N, Boyd D (1967) Series representations of quadratic forms in normal variables, I. central case. Ann Math Stat 38(3):823–837 Leick A (2004) GPS satellite surveying, 3rd edn. Wiley, New York Misra P, Enge P (2004) Global positioning system: signals, measurements and performance. Ganga-Jamuna Press, Lincoln Murphy T, Geren P, Pankaskie T (2007) GPS antenna group delay variation induced errors in GNSS based precision approach and landing systems. In: Proceedings of the institute of navigation GPS 2974–2989, September 2007 Rife J, Pullen S (2005a) Impact of measurement biases on availability for category III LAAS. In: Proceedings of the institute of navigation Annual Mtg 759–773, June 2005 Rife J, Pullen S (2005b) Impact of measurement biases on availability for category III LAAS. Navigation 50(4):215–228 Rife J, Pullen S, Enge P, Pervan B (2006) Paired overbounding for nonideal LAAS and WAAS error distributions. IEEE Trans Aerospace Electron Syst 42(4):1386–1395 RTCA (2004) Minimum Aviation Systems Performance Standards (MASPS) for the Local Area Augmentation System (LAAS), RTCA DO-245A RTCA (2006) Minimum Operational Performance Standards (MOPS) for GPS Wide Area Augmentation System (WAAS) Airborne Equipment, RTCA DO-229D RTCA (2007) Minimum Operational Performance Standards (MOPS) for GPS LAAS Airborne Equipment, RTCA DO-253B
GPS Solut Ruben H (1962) Probability content of regions under spherical normal distributions, IV: the distribution of homogenous and non-homogenous quadratic functions of normal variables. Ann Math Stat 33(2):542–570 Ruben H (1963) A new result on the distribution of quadratic forms. Ann Math Stat 34(4):1582–1584 Sheil J, O’Muircheartaigh I (1977) AS 106: the distribution of nonnegative quadratic forms in normal variables. Appl Stat 26(1):92–98 Shively C, Hsiao T (2004) Availability enhancements for cat IIIb LAAS. Navigation 51(1):45–47 Skidmore T, van Graas F, Liu F (1996) Ohio University LAAS monitoring: design and performance. Proceedings of the institute of navigation GPS 1615–1621, September 1996 Suddapalli R (2004) Aircraft position integrity for differential satellite-based navigation in the presence of both bias and noise errors. M.S.E.E. Thesis, Ohio University Therrien C (1992) Discrete random signals and statistical signal processing. Prentice-Hall, Englewood Cliffs van Graas F, Krishnan V, Suddapalli R, Skidmore T (2004) Conspiring biases in the local area augmentation system. In: Proceedings of the institute of navigation Annual Mtg 300–307, June 2004 Wackerly D, Mendenhall W, Scheaffer R (2002) Mathematical statistics with applications. Duxbury, Pacific Grove
Author Biographies Dr. Dean Bruckner is the Assistant Director (Technical) of the Avionics Engineering Center at Ohio University in Athens, Ohio. He received his Ph.D. degree in Electrical Engineering from Ohio University in 2010 with a dissertation that included flight test results for dual-frequency GPS precision approach and landing. During a career in the United States Coast Guard, he earned a B.S.E.E. from the Coast Guard Academy and an M.S.E.E. from the Naval Postgraduate School, both with honors and with specialization in electronic navigation. He served as Commanding Officer of the Coast Guard Loran-C transmitter station near Istanbul, Turkey, and was the nationwide electronics support manager for Coast Guard ships, boats, aids to navigation and shore facilities. He has recently provided navigation and engineering expertise to the F-35 Joint Strike Fighter navigation sensor test program.
Dr. Frank van Graas is the Fritz J. and Dolores H. Russ Professor of Electrical Engineering at Ohio University. He earned his Ph.D. from Ohio University in 1988. He holds a BSEE and MSEE degrees with specialization in avionics from Delft University of Technology, The Netherlands. His extensive involvement in satellite navigation included conducting the first real-time GPS attitude and heading flight experiment on a DC-3 in 1991, and the first kinematic dual-frequency GPS autoland flight tests using NASA Langley Research Center’s Boeing 737 in 1993. He received the Johannes Kepler Award for ‘‘sustained and significant contributions to satellite navigation’’ from the Satellite Division of the Institute of Navigation in 1996, and served as President of the Institute of Navigation from 1998 to 1999. At Ohio University’s Avionics Engineering Center, his roles include Principal Investigator for the FAAsponsored Global Positioning System (GPS) research grant to investigate the Local Area Augmentation System for Aircraft Precision Approach & Landing and Aircraft Surface Movement Guidance. He has lectured extensively throughout the United States, Canada, Europe and Russia for such organizations as the International Air Transport Association, the NATO Advisory Group on Aerospace Research and Development (AGARD), Navtech Seminars, Inc., and the International Society for Optical Engineering (SPIE). Dr. Trent Skidmore is a Senior Research Engineer and Adjunct Assistant Professor in the Avionics Engineering Center at Ohio University. He is involved in a variety of projects involving applications of GPS, including the use of differential GPS for aircraft precision approach applications. Recent projects include the Department of Defense Joint Precision Approach and Landing System (JPALS), the Automated Aerial Refueling (AAR) program, as well as the Federal Aviation Administration Local Area Augmentation System (LAAS). Trent received his B.Sc. in 1986 (with Honors) and M.Sc. in 1988 from Michigan Technological University and his Ph.D. from Ohio University in 1991, all in electrical engineering.
123
