Geometric Approximations of Heavy-Tail Effects for Chi-Square Integrity Monitors Jason H. Rife and John Scott Parker, Tufts University

BIOGRAPHY Jason Rife is an Associate Professor in the Department of Mechanical Engineering at Tufts University in Medford, Massachusetts. He directs the Automated Systems and Robotics Laboratory (ASAR), which applies theory and experiment to characterize integrity of autonomous vehicle systems. He received his B.S. in Mechanical and Aerospace Engineering from Cornell University and his M.S. and Ph.D. degrees in Mechanical Engineering from Stanford University. John Scott Parker is currently a Ph.D. candidate in the Department of Mechanical Engineering at Tufts University in Medford, MA. As a research assistant in the Automated Systems and Robotics Laboratory (ASAR), he is currently conducting research on sensor noise analyses both for GPS and THz-frequency navigation. He received his B.A. from Boston University in Physics and his M.S. from Tufts University in Mechanical Engineering. ABSTRACT This paper considers the verification of integrity monitors used commonly in safety-of-life navigation systems. The paper introduces three conservative bounds (also known as overbounds) that can be used to compute missed-detection performance for nominally chi-square integrity monitors subject to random noise with non-Gaussian heavy tails. The three proposed overbounds rely on simple geometric shapes – cones, cylinders, and spheres – to ensure that missed-detection probability can be computed conservatively, reliably, and efficiently. The overbounds provide conservative estimates of missed-detection performance even when the exact form of the noise distribution is not known, so long as the noise CDF is radially symmetric and bounded above and below. Monte Carlo simulations suggest that the spherical overbound performs best for faults inducing a bias of moderate or large magnitude. The cylindrical overbound performs best when the fault-induced bias is small in magnitude. INTRODUCTION In safety-of-life navigation applications, integrity monitors are critical for detecting navigation-system faults that might otherwise cause hazardously misleading information. In aviation applications, for instance, it is critical to detect navigationsystem faults that might cause airplanes to collide with each other or to crash on landing. To protect aviation applications reliant on the Global Navigation Satellite System (GNSS), engineers have developed a number of GNSS integrity-monitoring systems. These systems include Receiver Autonomous Integrity Monitoring (RAIM) [1]-[2], Satellite-Based Augmentation Systems (SBAS) [3], and Ground-Based Augmentation Systems (GBAS) [4],[5]. As part of their functionality, RAIM, GBAS and SBAS all scan for faulted GNSS satellites producing anomalous navigation signals. A subset of RAIM, GBAS and SBAS monitoring capabilities have been implemented using quadratic monitors that operate on vectors of input signal input. An example of such a vector-based monitor is Signal Deformation Monitoring (SDM), which scans a bank of correlators to identify anomalous GNSS waveforms that might result in positioning errors of tens or hundreds of meters [6]-[12]. Integrity monitors like SDM analyze the vector input’s magnitude (squared), a quantity which typically exhibits noise with an approximately chi-square probability density function (PDF). Modeling a monitor’s noise PDF is important to verifying that the monitor meets its design requirements, including its specification for the maximum allowed missed-detection probability Pmd [13]. Given that the exact form of the noise distribution is only known approximately, it is common to verify that the design meets its specification by using an overbound, which is a conservative approximation of the actual noise PDF [14],[15]. As one example, prior work on overbounding has

1

considered how to convolve multiple (not necessarily Gaussian) error PDFs of the same form [16],[17]. As another example, the authors have conducted recent research on overbounding an unknown covariance for Gaussian-distributed SDM noise, an uncertainty that can cause the PDF of the monitor statistic to diverge from its expected chi-square model [18]-[23]. A significant gap in the research literature concerns conservative bounding of PDF heavy-tails for vector monitors with nominally chi-square noise distributions. Here the term heavy-tails refer to PDFs exhibiting a higher-than-expected probability of occurrence for large errors [24]. Heavy-tails are observed commonly in analyzing GNSS signals [25]-[28], including SDM signals [29]. The influence of heavy tails on the characterization of confidence bounds for GNSS position estimates, also known as protection levels, has been carefully studied [24],[30],[31]. Little is known about the impact of heavy tails on integrity-monitor performance, however, particularly for vector monitors with nominally chi-square noise distributions. The key contribution of this paper is to introduce three new approaches for conservatively estimating the missed-detection probability Pmd for nominally chi-square monitors subject to heavy-tail noise distributions. Before introducing these three approaches, the paper first reviews the fundamentals of chi-square monitor operation and performance verification. Next, the paper reviews standard methods for warping one PDF to convert it into another PDF; these methods are applied to warping a radially symmetric non-Gaussian (i.e. heavy-tail) distribution into a radially symmetric, identity-covariance Gaussian distribution. As a next step, the paper proposes new methods for approximating the threshold surface that results from the warping process. Three geometric approximations are considered, including a conical approximation, a cylindrical approximation, and a spherical approximation. Methods are introduced to enable efficient computation of Pmd using these three approximate thresholds, and these methods are compared using a Monte Carlo simulation. Following the discussion of these results, a brief summary concludes the paper. BACKGROUND: CHI-SQUARE MONITORS In this paper, the term chi-square integrity monitor is used to describe a test that identifies when a random input vector x   N becomes anomalously large [19]. If the mean value of the x vector is nominally small, then the vector’s elements should rarely be large, unless a fault occurs. A convenient approach to identify a fault is therefore to compute a quadratic monitor statistic s,

s  xT x .

(1)

The scalar s is compared to a threshold at each moment in time to detect anomalies. The monitor alerts that a fault may be present whenever s exceeds the threshold T.

s  T  alert

(2)

Ideally, if the elements of the random vector x are conditioned to be independent and identically distributed Gaussians with unit variance, then the distribution of the monitor statistic s has a chi-square distribution. This ideal behavior is the reason chisquare monitors are so named. If, as considered in this paper, the random vector inputs are only approximately Gaussian distributed, then the monitor statistic, also, will only be approximately chi-square distributed. Non-Gaussian effects can make it difficult to assess monitor performance accurately, which in turn makes it a challenge to verify that the system satisfies its design specifications In designing a chi-square integrity monitor for safety-of-life applications, a key performance metric is the probability of missed detection. The goal of the monitor is to trigger an alert whenever a fault is present, so that an operator can take action to preserve overall system safety. Due to random noise on the monitor statistic, however, detection will occasionally fail. The missed-detection probability Pmd quantifies the likelihood of a detection failure given that a fault is present. To demonstrate that a monitor is effective, its Pmd must be assessed and compared to an upper bound Preq, a requirement derived from an overall system-risk assessment. The monitor is considered to meet the specification so long as the assessed risk is smaller than the requirement:

Pmd  Preq .

(3)

An accurate assessment of Pmd can be computationally intractable and is not necessarily required to show the monitor meets its requirement. Rather, it is sufficient to approximate Pmd conservatively using an upper bound Pmd . Given that an upper bound exists such that Pmd  Pmd , the design requirement is met so long as

2

Pmd  Preq .

(4)

Representing the actual Pmd with a conservative approximation is an example of an overbound [19]. A carefully selected overbound Pmd can greatly simplify analysis, particularly when the performance requirement must be evaluated repeatedly, thousands of times, to evaluate different monitor configurations or environmental configurations. Of course the overbounding approximation cannot be too overly conservative. An approximation is most useful if it is tight, meaning that the difference between Pmd and Pmd is small. If the difference is too large, then it becomes very difficult to satisfy the requirement (4). As an extreme example, Pmd is known to be less than one (since all probabilities lie between 0 and 1 inclusive); thus a Pmd overbound of 1 is a conservative approximation, but not a useful one (as this loose approximation implies a worst case where the monitor fails to detect each and every fault). Thus there is subtlety in selecting an overbounding strategy for any particular application. The overbound must be i) guaranteed conservative, ii) simple to implement in order to streamline repeated calculations, and iii) as tight as possible, so that the approximation does not cause the system to fail its design requirements. The primary theme of this article is to compare different methods of conservatively approximating Pmd for the quadratic monitor statistic (1) given that the random vector input x is non-Gaussian distributed. According to (2), a missed-detection event occurs when a fault is present yet the monitor statistic s fails to exceed the threshold T. Thus, the probability of missed detection can be obtained by integrating ps(s), the PDF for the monitor statistic, over all values of s below the threshold, given that a fault is present.

Pmd (T ) 

ps  s | fault  ds



(5)

m T

Using (1), Pmd can be equivalently expressed in terms of the fault-case distribution for the random input vector x, labeled px(x | fault).

Pmd (T ) 



px  x | fault  dVx

(6)

xT x T

The second form is obviously more complicated than the first, in that (5) is an integral in one dimension, whereas (6) is an integral over an N-dimensional volume. For this reason, it is desirable to employ the first form when possible to visualize and compute Pmd. If the random vector x is well characterized and conditioned so its elements are all independent and Gaussiandistributed with zero mean and unit variance, then (5) can be evaluated as a chi-square cumulative distribution function (CDF). However, in the more general case when the distribution of x is uncertain or non-Gaussian, then no simplified form of the distribution (5) exists, so Pmd must instead be modeled using (6). Because this paper considers random vectors x with nonGaussian PDFs, integral (6) must be evaluated. For this reason, methods are sought to streamline the computation of (6). Computing the Pmd integral in higher dimensions is challenging, particularly when Pmd is small (as is desired for a good monitor). Accurate direct numerical integration of (6) is very slow, in large part because roundoff errors quickly grow large when summing probability values spanning many orders of magnitude. Likewise, accurate probabilistic integration [32] is also very slow. In pursuit of computational efficiency, it is convenient to seek a solution that is approximate yet still guaranteed to be an upper bound for the integral (6). As long as the approximation is provably conservative, it can be used to verify that the monitor satisfies requirement (4). One mechanism of assuring a conservative probability model is to define the approximate missed-detection probability using an expanded volume of integration, without changing the integrand. Consider the case when the true Pmd is described by an integral similar to (6), evaluated over more general volume  x . Pmd (T ) 



px  x  dVx

(7)

x x

As with all missed detection analyses, a fault is assumed to be present; accordingly, the conditional probability notation in (6) has been suppressed for compactness. The above integral (7) can be approximated as

3

px  x  dVx ,



Pmd (T ) 

(8)

x x

where the original volume of integration  x is replaced with a larger volume  x that contains the first, such that  x   x . It is easy to show that Pmd  Pmd when the new volume of integration contains the old. To show this, expand the integral of (8) into two parts.



px  x  dVx 

x x



px  x  dVx 

x x



px  x  dVx

(9)

x y /  y

The final term in the equation represents the region of  x not included in the original volume  x . This second term is positive or zero, since probability density is by definition greater than or equal zero. Thus,  x   x implies that



px  x  dVx 

x x



px  x  dVx ,

(10)

x x

so Pmd  Pmd . Selecting the form for the approximation  x is a difficult design decision. The key requirement for safety-critical applications, is that the approximation must be conservative, to ensure Pmd  Pmd , but ideally the choice of  x would greatly simplify Pmd computation and enable robust bounding even when the noise PDF is not known precisely. Importantly, the overbound Pmd should also be a tight bound, as excess inflation can result in the monitor failing requirement (4). An extremely useful technique to enhance the computational efficiency of evaluating integral (8), when the PDF px(x) is nonGaussian, involves converting the integral to a Gaussian form. As such, techniques for warping a non-Gaussian distribution into a Gaussian distribution are described in the following section. THE WARPING FUNCTION

The idea of transforming (or warping) one distribution into another distribution is a well-understood concept in probability theory. Warping can be used, for example, to convert a set of samples created by a uniform random number generator into a set of samples characterized by an entirely different PDF (such as a chi-square PDF). A useful principle in converting one distribution to another is to match their CDFs [33]. The matching process integrates two PDFs to obtain CDFs, and finds the thresholds that result in matching probabilities for both distributions. For example, consider two random variables, x and y, which are respectively non-Gaussian and Gaussian distributed. Labeling the PDFs of these random variables png(x) and pg(y), their respective CDFs Png ( X ) and Pg (Y ) can be defined for thresholds X and Y as follows.



Png ( X ) 

png  x  dx

x X

Pg (Y ) 



pg  y  dy

(11)

y Y

By setting these two CDFs equal, we can convert the coordinates for a non-Gaussian CDF into coordinates for a Gaussian CDF.

Y  Pg1  Png ( X ) 

(12)

In this paper, functions of the form of (12) are called warping functions. The notation f will be used to represent the warping function, which simplifies (12) to be

Y  f (X ) .

(13)

4

Using the warping function, the non-Gaussian CDF from (11) can be computed from a Gaussian PDF as pg  y  dy .



Png ( X ) 

(14)

y f ( X )

Though the warping process stretches the entire PDF, the warping only appears in the limits of integration for the CDF. The same concept can be extended to higher dimensional vector spaces using integrals defined as a function of distance from the origin. This operation is particularly straightforward when considering radially symmetric distributions, as will be assumed in the remainder of this paper. Consider the case of mapping a radially symmetric, non-Gaussian distribution to a radially symmetric, Gaussian distribution. Let the non-Gaussian distribution be defined for a random vector x   N and the Gaussian distribution for a second vector y   N .

 p  x  dV

Png 

ng

x

Pg 

 p  y  dV

(15)

g

y

In these integrals the volume element is dV and the domain of integration is defined as  x for the non-Gaussian distribution and as  y for the Gaussian distribution. In order to map the non-Gaussian distribution into a Gaussian distribution, it is useful to warp only the radial coordinate, leaving all the angular components identical for both domains of integration  x and  y . As high-dimensional integrals are somewhat abstract, it is instructive to work through a specific two-dimensional example in more detail. For our specific example, consider a pair of two-dimensional distributions: one exponential and the other Gaussian, both with zero mean. The PDFs are pe  x  and p g  y  , respectively.

1  1L x e 2 L2  1 2  yT y  1 pg  y   e 2 2 2 pe  x  

(16)

Here the parameter L is the scaling parameter for the exponential distribution, and the parameter σ is the standard deviation of the Gaussian distribution. Both PDFs are radially symmetric, so the vector arguments x and y can be replaced with scalar radius values r and ρ, respectively. The CDFs can then be evaluated in polar coordinates to a threshold radius of R for the non-Gaussian distribution and of Rw for the Gaussian distribution.

 R   1 1 r    e L rdr d 2      0  r  0 2 L  2

Pe ( R) 

 Rw   1  1 2   2   Pg ( Rw )     e 2  d  d 2     0    0 2  2

(17)

Because the radial limits of integration, R and Rw, are defined to be functions of angle, the angular limits of integration can be defined identically for both integrals (over the range from 0 to 2 ). The radial integrals (in parentheses above) can now be related by a warping function. Using the functions PR ( R) and PW ( Rw ) to indicate the non-Gaussian and Gaussian radial integrals, respectively, (17) can equivalently be written as follows.

5

2

Pe ( R) 

PR  R   d



 0 2

Pg ( Rw ) 



(18)

PW  Rw   d

 0

For the two-dimensional exponential and Gaussian PDFs, the radial integrals can be calculated analytically.   R  R  1 PR ( R)  1  1   e L    L  2 R   12 w2 PW ( Rw )  1  e    2

(19)

 1   2 

(20)

Setting these two radial integrals equal and solving for Rw  f ( R) gives the warping function

f ( R)   2

R  R  ln 1   . L  L

(21)

This warping function (21) is illustrated in Fig. 1.

5

4

3

2

1

0

0

1

2

3

4

5

R/L Fig. 1. Warping function. Each point on the solid blue curve represents a particular probability. For the 2D exponential CDF, the normalized radius that achieves a particular probability is R/L. For a 2D Gaussian CDF, the normalized radius that achieves the same probabilty is f(R)/σ.

6

Equating the two integrals from (17), the exponential CDF can be computed by integrating a Gaussian PDF, so long as the bound of integration is warped through f ( R) .

 f ( R   ) 1  1 2   2     e 2  d  d 2      0  r  0 2  2

Pe ( R) 

(22)

If the warped domain of integration described by f(R) can be approximated intelligently as a larger volume, then the Gaussian integral can be evaluated efficiently using standard software tools, for example using existing functions available in Matlab. That a non-Gaussian CDF can be constructed by evaluating a Gaussian CDF at a warped radius is the key point of the paper. To help emphasize this point, Fig. 2 shows the 2D exponential CDF of (19) constructed by evaluating the Gaussian CDF (20) with an argument equal to the warped radius from (21). In other words, PR ( R)  PW ( f ( R)) .

1 0.9 0.8 0.7 0.6 0.5 0.4 0.3 2D Gaussian: P R (R)

0.2

2D Exponential: P W (R)

0.1 0

Gaussian with Warped Argument: P R (f(R))

0

0.5

1

1.5

2

2.5

3

3.5

4

4.5

5

R Fig. 2. Warping converts a Gaussian CDF into an Exponential CDF.

Note that, in this example, the non-Gaussian PDF was selected purposefully to be a two-dimensional exponential PDF to ensure the integrals of (17) were analytically calculable. In general, the warping function has no analytic expression and must be evaluated numerically. Extrapolating from this two-dimensional example, it is possible to define a general approach for higher dimensional integrals and arbitrary non-Gaussian distributions. By analogy to the two-dimensional example, the warping function f(R) can be defined to relate two radial integrals.

f ( R)  PW1  PR ( R) 

(23)

In N dimensions, these radial integrals are the following.

7

PR ( R) 



pr  r  r N 1dr

rR

PW ( Rw ) 



(24)

p     N 1d 

  Rw

The key to defining the two radial integrals is to restrict both the non-Gaussian and Gaussian PDFs to be radially symmetric functions pr  r  and p    , where the radial coordinates r  x and   y refer to the non-Gaussian and Gaussian PDFs, respectively. The original probability integrals (15) can be rewritten in terms of the radial CDFs of (24) and a set of N-1 angles, labeled 1 ,, N 1 . After integrating in the radial direction, the remaining angular integration element is

dS  sin N  2 1  sin N 3 2  sin N  2  d1d2  dN  2 dN 1

,

(25)

and the integrals (15) take the following form.

Png   PR ( R 1 , , N 1 )dS 

Pg   PW ( Rw 1 , , N 1 )dS

(26)



Here the angular domain of integration is defined as

  1 ,, N  2  j  [0,  ]  N 1   0, 2  .

(27)

Substituting warping function (23) and equating the two probabilities in (26) gives the non-Gaussian CDF in terms of the Gaussian radial integral PW.





Png   PW f  R 1 ,, N 1   dS . 

(28)

This Gaussian integral can be approximated by conservatively inflating R 1 ,, N 1  so that the approximate volume of integration contains the original volume, as described for (8) above. In defining an approximation of the volume of integration, it is useful not only to consider computation efficiency but also robustness, in that the approximation can be defined to bound not just a particular non-Gaussian noise PDF, but a whole range of possible non-Gaussian noise PDFs. One way to ensure that an analysis bounds a wide range of possible non-Gaussian distributions is to bracket the possible values of the warped radius f(R) for each radius R. For instance, an upper bound f ( R) and a lower bound f ( R) might be defined, such that

f ( R)  f ( R)  f ( R) .

(29)

As long as an analysis can show that all warped radii on this range are conservatively bounded, then the shape of the actual noise distribution does not matter! This is equivalent to the conventional notion of paired CDF overbounding [15], where the exact form of a noise distribution may be unknown, so long as an overbound and an underbound CDF can be defined. A subtle detail here is the lower-bound radius f ( R) corresponds to the conservative overbound CDF. This detail can be observed by comparing Fig. 1 and Fig. 2. In Fig. 1, the warped radius f(R) for the heavy-tail exponential is less than R; however the CDF bound in Fig. 2 indicates that the exponential CDF is clearly more conservative (with greater probability in the tail) as compared to the Gaussian. To summarize briefly, we have now defined two key steps for conservative approximation of Pmd for integrity monitors with non-Gaussian vector inputs of N dimensions. One of these steps is warping, which converts the non-Gaussian integral to a

8

Gaussian integral. The other step is approximation of the integration volume, a step which can both streamline computation and promote bounding of an entire range of input distributions, as characterized by (29).

CHARACTERIZING THE WARPED THRESHOLD AS A SURFACE OF REVOLUTION In order to define approximate volumes of integration that enable robust and efficient Pmd computation, it is important that the approximate volume be a tight bound. To support tight bounding, this section characterizes the shape of the threshold for a chisquare monitor as being a surface of revolution. Though the shape of the warping function f(R) is arbitrary, this section shows that the warped threshold is still always a surface of revolution about a particular axis, one which lies in the direction of the fault-induced bias. As a starting point in our discussion of surfaces of revolution, consider that the original threshold surface is a sphere. To be precise, the original threshold surface encloses the volume xT x  T , as described by (6). This spherical volume is symmetric about any axis through the origin, including the bias axis, a term that will be used to describe the axis parallel to the mean vector μ x  E  x  . The importance of the bias axis traces back to the warping process, described by (23) and (24), equations which assume the non-Gaussian PDF is radially symmetric. When a nonzero mean is present, due for instance to the occurrence of a fault, then the distribution px(x) loses its radial symmetry about x  0 , even if the distribution remains radially symmetric about the mean μ x . The loss of radial symmetry means that the warping transformation, as described by (15), cannot be directly applied to the Pmd integral (6). Given that px(x) is radially symmetric about its mean, the radially symmetric form of the integral (15) can be recovered by introducing the change of variables v  x  μ x . This change of variables recovers a radially symmetric PDF

pv  r   px  v  μ x  ,

(30)

r v ,

(31)

where

and transforms (6) into an integral over the spherical volume  v , for which the boundary is

 v   v  x  μ x

x T

.

(32)

This boundary is constructed from all points x on a sphere of radius T about the origin of the x-coordinate system. These points are shifted by an offset μ x , which moves the sphere’s center away from the origin of the v-coordinate system. Integrating within this boundary gives Pmd (T ) 



pv  r  dVv .

(33)

v v

This integral can now be warped into a Gaussian form by CDF matching, as described by (15). The result is

Pmd (T ) 



pg    dVy .

(34)

y y

The integration volume  y for the Gaussian integral (34) has a surface  y obtained by scaling the radial distance of each point on the original threshold surface  v . In other words, the warped threshold surface is



 y  y  g  v  v



v   v .

(35)

9

The scale factor g  v



maps the radial coordinate of the boundary of integration from its original value at R  v to a new

value Rw  f  R  .

 f ( R ) / R, R  0 g ( R)   . R0  0,

(36)

Different scalings g can be introduced to represent the range of possible distributions described by (29). Combining (29) and (36) gives the following limits on g, written in terms of upper and lower bounds g ( R) and g ( R) .

g ( R)  f ( R)  g ( R )

(37)

An interesting characteristic of the warping function is that it merely translates points v   v along a ray from the origin (directly inwards along that ray for a heavy-tail distribution). The warping process does not change the direction of the ray. Because the warping is a simple scaling function, the warped threshold surface  y retains some of the structure of the original threshold  v . The structure of the spherical surface  v is perhaps best described as axially symmetric about a line through the origin O (at

v = 0) and the threshold center at v  μ x . This use of the term axial symmetry is distinct from our earlier use of the term radially symmetric to mean symmetric about any axis, as in the case of  x , which is radially symmetric about the point x = 0. The N-dimensional interpretation of axial symmetry is that, for any hyperplane orthogonal to the axis of symmetry, all the points at the intersection of that hyperplane and the threshold surface must be equidistant from the axis of symmetry. By extension, the threshold surface can be described by two coordinates – one along the bias axis and one normal to it. Moreover, the original or warped threshold surface can be constructed by taking a two-dimensional curve-of-revolution, as shown in Fig. 3, and sweeping it around the bias axis. The warped curve-of-revolution, in Fig. 3(b), is simply a scaled version of the original curve-of-revolution, in Fig. 3(a), where the scaling g, as described by (35), moves each point on the original curve inward toward the origin O.

Fig. 3. Curve-of-revolution for (a) original spherical threshold and for (b) warped surface.

10

More than one option exists in choosing the parameter pair used to represent a surface of revolution, like  v or  y . One choice might be to use two spherical coordinates: the radial coordinate r and the first angle 1   0,   . The advantage of this formulation is that it is relatively easy to understand how the curve-of-revolution can be swept through the full domain of all the remaining spherical-coordinate angles 2 ,, N 1 , where these domains are defined by (27). Though it is hard to visualize this process of sweeping out a curve-of-revolution beyond three-dimensions; the mathematics extends in a straightforward manner to higher dimensions. Another useful representation of the surfaces  v or  y is to quantify each in terms of a pair of orthogonal coordinates, one describing distance along the bias axis and one describing distance normal to the bias axis. Both terms are conveniently parameterized using the angle    0,   , which is an arc length around the original spherical threshold, as illustrated in Fig. 3(a). Parameterizing by θ has an advantage over parameterizing by 1 , in that each point around the curve-of-revolution corresponds uniquely to a particular value of θ. Let’s begin by defining the curve-of-revolution in the x-coordinate system, where the threshold surface is uniform distance from x = 0, and the bias-axis coordinate x1 and normal coordinate xn are

x1  T cos  xn  T sin 

.

(38)

Moving to the v-coordinate system, shifted by the bias μ x such that v  x  μ x , the spherical threshold can be expressed in terms of axial and normal coordinates, v1 and vn respectively, such that

v1  T cos   b vn  T sin 

.

(39)

Here the bias magnitude is defined as b  μ x . It should be noted that, after warping, the parameter θ no longer has a geometric interpretation as an angle. However, it is nonetheless a useful parameterization describing progress between the endpoints of the curve-of-revolution. A warped threshold, as illustrated in Fig. 3(b), can be constructed from the v-coordinates by a scaling operation, as described by (35). The resulting y-coordinates for the warped threshold include y1 in the bias direction and yn normal to it. y1  g  R 



T cos   b



yn  g  R  T sin 

(40)

Here R refers to length of the v vector, such that

R2  T  b2  2b T cos 

.(41)

For cases with a fault present, the direction of the fault-induced bias is arbitrary, given that the actual implementation of the threshold and alert in (2) is radially symmetric and makes no assumptions about the bias direction. Although the bias direction is considered arbitrary, the bias magnitude b will be constrained. In this paper, the magnitude b is assumed large enough to place the expected value of the monitor statistic outside the threshold surface. To be precise, it is assumed that

b T .

(42)

The implication of (42) is that v1  0 and y1  0 .

11

With the warped threshold surface  y characterized as a simple two-dimensional curve-of-revolution (through the coordinates

y1 and yn ), it is now possible to approximate a surface  y as a curves-of-revolution, too. As long as the approximate curveof-revolution encloses all values of y1 and yn on the warped threshold, then  y   y and an approximate Pmd can be defined from (34) as

Pmd (T ) 



pg    dVy

(43)

y y

with the guarantee that Pmd  Pmd . With this goal in mind, the following three sections will describe three approximations of the warped surface based on different simple geometries like cones, cylinders, and spheres – which all can be constructed from curves-of-revolution (in the form of wedges, rectangles, or semicircles).

CONICAL OVERBOUND The first upper-bound for Pmd is constructed by introducing an approximate volume of integration  y using an N-dimensional cone constructed as a surface-of-revolution about the bias axis. To provide a tight fit to the warped volume of integration  y , the cone is truncated at a minimum and a maximum radius, as shown in Fig. 4. The figure illustrates the curve-of-revolution for the truncated cone overbound as a thick green contour. In comparison a representative curve-of-revolution for the warped threshold is illustrated as a thin blue contour. The overbound surface is constructed as a tight fit by defining it to be tangent to the warped surface at three points: at both endpoints (where θ is 0 or π) and at a critical point in between the endpoints (a point that will be quantified below). For additional perspective, a three-dimensional view of the warped and approximate surfaces of revolution is shown in Fig. 5. The illustration shows the warped threshold surface  y , followed by a cutaway of the bottom half of the conical overbound

 y , followed by a depiction of the two surfaces nestled together, with the warped threshold surface fitting tightly inside the truncated cone. The most salient feature of the cone is its spreading angle relative to the bias axis. This maximum spreading angle, labeled here as * 1 is defined so that the cone is tangent to the warped surface at exactly one point over the range of the parameter    0,   , as shown in Fig. 4. This tangent point can be computed by noting that

1  argmax  tan 1 

(44)

  yn    y1 

(45)

*

  0, 

and expanding the tangent term as follows.

1  argmax 

*

 (0, )

The corresponding maximization problem can be solved using a calculus-based solution. Setting the derivative of the tangent to zero to look for extrema gives

   yn    y1

 0. 

(46)

Substituting (40) into the above equation and solving results in * , the extremum of θ.

 T  b

  cos 1 

*

  

(47)

12

Fig. 4. Curve-of-revolution for the conical overbound (thick green contour) and for the warped threshold surface (thin blue contour)

Fig. 5. Three dimensional representations of (a) warped threshold surface, (b) cutaway of overbounding volume defined as a truncated cone, (c) warped threshold surface nestled in overbounding volume.

Evaluating (40) and (45) at * , the maximum value of the spreading angle is



  .  b / T 1 

1  tan 1 

*

1

2

(48)

This interior extremum is unique, because the cosine function in (47) is monotone over the range    0,   . Although extremum can also occur at the endpoints of the domain, the endpoints turn out to be minima of 1 (and are handled separately).

13

As an interesting side-note, the same spreading angle * 1 describes the line that is tangent to both the original threshold (in v coordinates) and the warped threshold (in y coordinates). These matched values of * 1 are illustrated as dashed lines in both in Fig. 3(a) and Fig. 3(b). That the tangent condition for the original and warped threshold occur at the same angle * 1 is somewhat intuitive, because the scaling function g(R) only changes the distance of each point relative to the origin at O, without changing angle. The same result can be obtained analytically, by noting that (45) can be expanded using (40), to give





 yn  g  R  T cos   b ,  y1 g  R  T sin 

(49)

that the scaling factor cancels out of this ratio, and that (39) can be substituted into this expression to give  yn vn .  y1 v1

(50)

Once the N-dimensional cone is defined as a surface of revolution based on the spreading angle * 1 , a truncation of the cone can be defined to provide a tight bound at the inner and outer radii of the warped surface. In truncating the cone, the maximum and minimum radii are set to be the values f  R  for the two endpoints, where   0,   . To show that the two endpoints of the curve-of-revolution are the extreme values of f  R  , consider the following derivative, expanded using the chain rule.

f  R  R  f  R      R 

(51)

This derivative can be shown to be always greater than zero. As a first step, expand the first term using (23) to obtain

f  R  R



PW1  PR  PR ( R) PR

R

.

(52)

The two terms above are both non-negative. As for the first term, note that the function PW1  PR  is the inverse of the radial Gaussian integral, which monotonically increases with probability PR . Hence the derivative of PW1  PR  with respect to PR is strictly positive. The second term is the derivative of the non-Gaussian radial integral PR ( R) , given by (24). Taking the derivative of (24) yields a product of a PDF (non-negative) and the radius R (positive).

PR ( R)  pNG  R  R N 1 . R

(53)

Since both terms of (52) are non-negative, it is true that

f  R  R

0.

(54)

Returning to (51), we must still consider the derivative of R . The derivative in question can be obtained from (41).

 b T sin   R   R 

(55)

The parameters b and T are positive, the sine function is non-negative on the domain   [0,  ] , and the radius R is defined to be positive. Thus, it is true that

14

  R   0 , 

(56)

  f  R    0 . 

(57)

and hence from (51) and (54),

Because the warped radius f  R  never decreases on the domain   [0,  ] , its extrema can be found at the endpoints, with a maximum at    and a minimum at   0 . Evaluating R at its extrema using (41) allows the upper and lower bounds for

f  R  to be obtained in terms of b and T. The result is that



 



f  R    f b  T , f b  T  .  

(58)

If the specific form of the actual distribution is unknown but upper and lower bounded by (29), then (58) can be generalized to span this range of possible distributions. To this end, the inner and outer radius limits can be defined as ρi and ρo respectively, such that

f  R    i , o 

(59)

 .  f b  T 

(60)

with

i  f b  T o

Using the limits of integration defined by (48) and (60), the conical bound can be evaluated in the form of (43). Integrating the unit-variance, zero-mean Gaussian joint PDF inside the cone with spreading angle * 1 and truncated to an inner radius ρi and an outer radius ρo gives 1* o

Pmd (T )     0 

 p   

N 1

g

d  sin N  2 1d1dS .

(61)

i

 refers to the corresponding domain of Here dS refers to all of the angular elements except those associated with 1 ; the set 

integration for those angles. Removing the 1 from (25) gives dS  sin N 3 2  sin N  2  d2  dN  2 dN 1

(62)

   ,,   2 N  2  j  [0,  ]  N 1   0, 2  .

(63)

and from (27) gives

Since the PDF is radially symmetric, it is only a function of ρ, and so the 1 terms can be factored out of (61) to give Pmd (T )   1*  

o

 p   g

N 1

d  dS ,

(64)

dS i

where

15

1*

 1*    sin N  2 1d1 .

(65)

0

Equation (64) can be related to the chi-square distribution P 2 ( R 2 ; N ) with N degrees of freedom and threshold R2. R

P 2 ( R 2 ; N )       pg     N 1d  dS

(66)

dS 0

Note, it is assumed that the Gaussian distribution pg has an identity covariance matrix and zero mean. Evaluating (66) for two radii, an inner radius ρi and an outer radius ρo , gives o

  p   g

N 1

d  dS 

dS i

P 2 ( 02 ; N )  P 2 ( i2 ; N )   

.

(67)

Combining (64) with (67) gives

Pmd (T ) 

 1*    

P

2



( 02 ; N )  P 2 ( i2 ; N ) .

(68)

Evaluating this equation with the definitions of the function  from (65), the parameter * 1 from (48), and the parametersρi andρo from (60) produces the bounding value Pmd . For the most part, this functional form is very quick to compute since the only integrals needing evaluation are one-dimensional. The N-dimensional character of the bound is captured through P 2 , which is a conventional chi-square CDF of N degrees of freedom. The conventional chi-square CDF can be evaluated using standard tools in common analysis packages (for instance, using the chi2cdf command in Matlab).

CYLINDIRICAL OVERBOUND This section considers a second method for conservatively approximating the warped surface  y . Specifically, this section develops a bound using cylindrical coordinates. In this case, the approximate volume of integration for (43) is obtained using a box-shaped curve-of-revolution, as illustrated in in Fig. 6. The figure depicts the warped threshold surface (thin blue line) tightly bounded by the cylindrical overbound (thick green contour). The key parameters for the cylindrical bound are the upper and lower limits along the bias axis, which define the domain y1   y1 , y1  , and the outer limit normal to the bias axes, which defines the domain yn    yn , yn  . The lower-bound y1 can be obtained analytically. This lower-bound occurs where the warped radius f  R  is largest, which according to (57) occurs at θ = π, where





y1   f b  T .

(69)

This point occurs on the left side of the warped surface illustrated in Fig. 6. By contrast, neither the upper-limit y1 (on the bias axis) nor yn (on the normal axis) can be described in a closed-form. Both must be obtained by numerical methods. Practically speaking, the easiest way to find both values is to use the variable definitions of (40) to perform an exhaustive line search over the range    0,   . The upper bound y1 is obtained from the smallestpossible scaling function g  R  , because y1 is negative; by contrast, the upper bound yn is obtained from the largest-possible scaling function g  R  , because yn is positive on the range    0,   . Thus the line searches have the following form:

16





y1   min  g  R  b  T cos      0,   yn  max  g  R  T sin     0, 

.

(70)

The line search of (70) cannot easily be simplified using calculus-based optimization; the reason is that the derivative of the warping function f  R  is strongly dependent on the character of the non-Gaussian distribution, making the derivative of y1 difficult to analyze without requiring a root-finding method that must be evaluated over the entire region    0,   .

Fig. 6. Cross-section for the cylindrical overbound (thick green contour) and for the warped threshold surface (thin blue contour)

Fig. 7. Three dimensional representations of (a) warped threshold surface, (b) cutaway of overbounding volume defined as a cylinder, (c) warped threshold surface nestled in overbounding volume.

17

Since the lower bound y1 occurs at θ = π, it might seem intuitive for the upper bound y1 to occur at the opposite end of the curve-of-revolution (where θ = 0). Indeed, this is sometimes the case, as illustrated in Fig. 6. However, in some cases, the warped surface is non-convex such that y1 briefly increases as θ increases near θ = 0, and so y1 occurs off the bias axis. An example of such a case is described in the next section, titled Spherical Overbound. Once the line search for y1 and yn is complete, the approximate missed-detection probability (43) can be evaluated. Assessing limits of integration using (69) and (70), the approximate Pmd can be expressed as Pmd (T ) 

 y1    pg    dy1  dy2  dyn .      y 2  y12  yn2  y1

(71)

The Gaussian distribution pg has an identity covariance matrix, which in turn implies that the distribution is independent along each of the orthogonal axes identified in (71). Factoring out the independent integral along the bias axis, equation (71) takes the following form.

 y1  Pmd (T )    1 pg  y1  dy1   y   1       N 1

 N 1 pg  y2 , , y N  dy2  dyn   

(72)

In the above expression, the leading superscript on each Gaussian distribution pg is introduced to indicate its dimension; the domain of integration for the second term is a sphere of dimension N-1:





2

 N 1  y | yn2  y  y12 .

(73)

On inspection, it is evident that the first term of integral (72) is Gaussian distributed and that the second is chi-square distributed. If the notation Pnorm is used to describe the standard-normal CDF, as in

Pnorm  y  

y





1 2

1 x 2

e 2 dx ,

(74)

then integral (72) can be written equivalently as



   P  y ; N  1 .

Pmd (T )  Pnorm  y1   Pnorm y1

2

2 n

(75)

As with the conical bound, this cylindrical bound is for the most part straightforward to evaluate. The primary complexity is in conducting the exhaustive one-dimensional searches for the bounds of integration, as described by (70). Evaluating the standardnormal CDF using these bounds is trivial, as most mathematical analysis packages provide tools for quickly evaluating Gaussian and chi-square integrals (for instance, using the normcdf and chi2cdf commands in Matlab).

SPHERICAL OVERBOUND This section explores a third geometric overbound, based on the notion of replacing the warped integration volume with a hypersphere in N dimensions. This spherical volume has a radius τ and a center at  y1 , yn    m,0  . The key parameters τ and m are defined to be positive scalars. An illustration of the spherical overbound is shown in Fig. 8, where curves-of-revolution for the warped threshold (thin blue curve) and the spherical bound (thick green curve) are plotted together. For perspective, a three-dimensional surface of revolution is illustrated in Fig. 9. Note that the spherical bound is tight on the right side (near the origin O), where probability is highest; however, the spherical bound is very loose on the left side, where probability is low.

18

Of the two parameters describing the spherical bound, the easier to describe is the center parameter





m   f b T .



The argument of f b  T

 is positive since b 

(76)

T , an assumption originally expressed as (42).

Fig. 8. Cross-section for spherical overbound (thick green contour) and warped threshold surface (thin blue contour). The region shown here is much broader than that shown in prior figures illustrating the conical and cylindrical bounds; to provide scale, the boxed region in the lower right corner indicates the smaller region shown in Fig. 4 and Fig. 6.

Fig. 9. Three dimensional representations of (a) warped threshold surface, (b) cutaway of overbounding volume defined as a sphere, (c) warped threshold surface nestled in overbounding volume.

19

Equation (76) defines the spherical center in such a way as to ensure that the sphere passes through the point on the warped surface along the bias axis whereθ = 0. Unless this point is the highest value of y1 on the warped surface, the spherical surface defined by (76) cannot contain the warped surface. As such it should be checked that no points on the surface are higher than y1  0 , or equivalently it must be confirmed that     0 for all    0,   , where

    y1  0   y1   .

(77)

An example of a case that violates this condition is illustrated in Fig. 10. The illustration depicts a nonconvex warped surface in which    goes briefly negative on either side of θ = 0. This figure was generated from the warping function for a twodimensional exponential distribution, given by (21), for the case of L = 1, b = 7.7 and T = 49. This is a case when the bias falls relatively near the threshold. To exclude cases like the one in Fig. 10,    can be evaluated using the following expression, which expands y1 in terms of (40) and which scales by the most extreme distributions in the range of (29) and (37).



 

    g ex  R  b  T cos   f b  T



(78)

The most extreme scaling function gex occurs either at the upper or lower limit, depending on the value of θ as compared to the location of the tangent point at * .

 g  R    [0, * ) g ex  R    *  g  R    [  ,  ]

(79)

The value of * is given by (47).

Fig. 10. A Warped Surface with Negative ∆(θ) Values.

20

So long as     0 for all    0,   , a spherical bound can be constructed with radius τ obtained by the following onedimensional exhaustive search.

 g ex2  R  T sin 2 

  max    0, 



2   



     2 

(80)

The expression in the bracket is the radius needed to bound the warped surface at any coordinate θ. By maximizing this radius, we ensure a bound for all θ. The derivation of this expression for τ is introduced in a companion paper [34], as the derivation is somewhat lengthy. A key detail in implementing (80) is that the expression inside the bracket has a step transition at * , where the evaluation of gex, as described by (79), and of    , as described by (78) switches from lower bounding to upper bounding. Using the center and radius parameters m and τ, the conservative missed-detection probability (43) can be approximated over a spherical domain of integration. The form of this integral is

Pmd (T ) 



y  meˆ1

pg    dVy . 2

(81)

 2

The boundary of integration is a sphere of radius τ centered on the point meˆ1 , where eˆ1 is a unit vector pointing along the bias axis in the μ x direction. By construction, this integral has the form of a noncentral chi-square distribution. Using the



notation Pncx  2 ; N , m 2



to represent the noncentral chi-square distribution integrated to a radius τ, and parameterized by N

degrees of freedom (corresponding to an N dimensional Gaussian PDF) and by a noncentrality parameter of m2, the integral (81) can be written as

Pmd (T )  Pncx  2 ; N , m2  .

(82)

Once again, this bound is for the most part straightforward to evaluate. The primary complexity is in conducting the exhaustive one-dimensional searches to check the sign of  and to maximize τ. Evaluating the noncentral chi-square CDF is trivial, as most mathematical analysis packages provide tools for evaluating noncentral chi-square CDFs (for instance, using the ncx2cdf command in Matlab).

OVERBOUND COMPARISON This section compares the performance of the three geometric overbounds: the conical and cylindrical bounds (derived in this paper) as well as the spherical bound (derived in companion paper [34]). The primary question is one of performance. How well do the various overbounds approximate the actual Pmd of (6)? With this question as motivation, a Monte Carlo simulation [32] was implemented to evaluate the actual Pmd for the warping function converting a two-dimensional exponential distribution into a two-dimensional Gaussian distribution, as described by (21). The actual Pmd integral (6) was estimated statistically by sampling 4107 random vector samples with a two-dimensional exponential PDF as described by (16). In the expression for the exponential distribution, the scaling parameter L was arbitrarily set to a value of 1. The samples were then compared to a spherical threshold with a radius T  7 . The simulation considered fifteen values of bias magnitude b, evenly distributed on the range between 1.1 T and 2.5 T . Bias





values were all chosen to lie outside the threshold b / T  1 , since this is a requirement for all three geometric overbounds proposed in this paper. For each bias value, the Pmd integral (6) was computed simply by summing all of the random samples occurring inside the threshold and dividing by the total number of samples. Biases were applied by shifting the threshold (which results in the same Pmd as would shifting the random samples in the opposite direction). In the process of computing the actual Pmd integral, the Monte Carlo simulation was also applied to validate the geometric overbounds. To do this, random samples were counted inside each of the proposed overbounding volumes – the cone, the

21

cylinder and the sphere. An illustration of the Monte Carlo simulation is shown in Fig. 11. The figure shows the warped surface of integration (left side) and three volumes that conservatively approximate the warped surface. These include the conical bound (bottom), which looks like a truncated wedge in 2D; the cylindrical bound (right), which looks like a rectangle in 2D; and the spherical bound (top), which looks like a circle in 2D. In the illustration, the overbounding surfaces have been rotated away from the warping surface for visualization purposes only. Though conceptually the orientation of the bias is arbitrary, since the probability distribution and overbounding volumes are all radially symmetric, all numerical integrals were performed on the left side, using the same data points to ensure the consistency of the integral result.

Fig. 11. Monte Carlo simulation of 2D Gaussian distribution showing missed-detections inside warped threshold surface (right), conical bound (bottom), cylindrical bound (right), and spherical bound (top).

The Monte Carlo integration obtained Pmd as a function of normalized bias b / T , as shown in Fig. 12. The figure identifies the actual missed detection cases (inside warped threshold) as red crosses. The figure indicates the Monte Carlo results for the three overbounds using other symbols (triangle for conical bound, square for cylinder bound, circle for spherical bound). The figure also shows the analytical results for Pmd for the three overbound cases. The analytical results were obtained by evaluating (68) for the conical bound, (75) for the cylindrical bound, and (82) for the spherical bound. These analytical results are plotted as dashed lines in Fig. 12. The result of the verification study is that the Monte Carlo estimates of the overbound probabilities closely match theoretical predictions. The similarity between the results for the two methods of computing the overbound probability – one statistical and the other probabilistic – are evident from the fact that the theoretical curves (dashed lines in Fig. 12) closely represent the MonteCarlo markers of the matching color. The results of Fig. 12 can be used to consider the primary question regarding the relative quality of the three bounds, at least for the specific case of a two-dimensional exponential distribution. In order to compare the bounds, the Monte Carlo integral was

22

normalized by each Pmd overbound to give a Pmd ratio as a function of bias magnitude b, as plotted in Fig. 13. This ratio of Pmd values should approach one for a tight bound. For a conservative bound, the ratio should never exceed one. In other words, a good overbound should go as close to the correct answer as possible without going over, much like a good guess on the gameshow The Price is Right.

Fig. 12. Pmd as a function of Normalized Bias. Results computed from Monte Carlo simulation (markers) are compared to theoretical predictions (dashed lines).

Of the three proposed overbounds, performance is poorest for the conical bound (mauve triangles). The conical bound might be described as “adequate” for small biases, since the Pmd ratio exceeds 80% on the left side of Fig. 13. The performance of the conical bound deteriorates quickly with increasing bias, moreover, with the Pmd ratio dropping below 55% on the right-hand side of Fig. 13. This decay in performance for large biases is directly related to the shape of the bound on the side facing the origin. The inner edge of the bound has two “ears” that project into the region of high probability near the origin. These “ears” increasingly diverge from the shape of the warped surface as bias increases, causing the performance of the overbound to worsen. The cylindrical bound performs better than the conical bound at all bias values. The tightness of the cylindrical bound is impressive, particularly for small biases, where the Pmd ratio approaches 95% on the left side of Fig. 13. The Pmd ratio for the cylindrical bound decays slightly as bias magnitude increases (below 80% on the right side of the plot). Although the deterioration of the cylindrical bound is not as pronounced as that of the conical bound, a similar mechanism causes the issue. Like the conical bound the cylindrical bound also has corners that project into the high probability region near the origin of the PDF – corners that become increasingly pronounced at higher bias magnitudes. The spherical bound is somewhat different than the other two bounds. The spherical bound is remarkably tight for moderate and high bias magnitudes; however, the spherical bound was not assessed for cases where the minimum    was below zero, a result indicating nonconvexity of the warped surface (see Fig. 10). For the two-dimensional exponential PDF, with a threshold

23

T = 49, the nonconvexity disappears for values of b / T above approximately 1.5. The spherical Pmd bound is only plotted above this transition point, as can be seen both in Fig. 12 and Fig. 13. Where plotted, the spherical bound achieves a Pmd raio everywhere above 90% in Fig. 13. In interpreting Fig. 13, it is important to remember that the actual Pmd value (used as the numerator in the Pmd ratio) was computed using a Monte Carlo integrator, and that the accuracy of this statistical integrator becomes poor as bias becomes large (because few points are left in the volume of the warped surface, as might be inferred from Fig. 11). The higher uncertainty of the Pmd ratio is shown visually by the widening of the 95% confidence bounds, which are shown as gray regions in Fig. 13. An interesting characteristic of the spherical bound is that the performance of the bound actually appears to improve as the bias magnitude b increases. This trend is evident in the middle of Fig. 13, especially in the range b / T  1.5, 2 . In this region, the Pmd ratio appears to approach slowly toward its ideal value of one. For larger biases, the confidence bounds grow too wide, too quickly to make a clear conclusion about whether the Pmd ratio continues to converge toward one. More work is needed in the future to confirm whether the spherical bound does indeed converge to the true answer as b grows. There is reason to suspect that convergence might occur, since the spherical bound is designed to be a very tight fit for the actual warped surface where probability is highest (i.e., at points closest to the origin of the radially-symmetric PDF).

Fig. 13. Comparison of Pmd Ratio for Three Overbounds. The ratio is assessed as the actual threshold, assessed by Monte Carlo simulation, normalized by overbound, assessed by theory. Nintey-five percent confidence bounds for the Monte Carlo simulation (obtained using the binomial distribution) are ploted in gray.

24

CONCLUSION The motivation of this paper is to aid in verifying the design of integrity monitors that ensure the safety of mission-critical navigation systems. In particular, this paper introduces three new methods for conservatively estimating the missed-detection probability of quadratic integrity monitors with nominally chi-square distributions. For monitors with vector inputs whose noise distributions are not precisely Gaussian, monitor performance cannot be described precisely by a chi-square distribution. To ensure a conservative assessment of missed-detection performance, this paper proposes geometric approximations, called overbounds, which transform the actual non-Gaussian distribution into a Gaussian form. The three overbounds introduced use simple geometric configurations to ensure bounds may be computed efficiently and reliably. The three simplifying geometries include conical, cylindrical and spherical volumes defined in N dimensions. The resulting expressions for missed-detection probability can be easily evaluated using functions available in commonly used mathematical analysis packages. The overbounds provide a clear result even when the form of the actual distribution is not known precisely (i.e. when the actual distribution can only be constrained to lie between upper and lower bounds). The key performance characteristic of the three proposed overbounds is tightness-of-fit. In other words, the overbound should be as close to the actual value as possible, while still always erring on the side of conservatism (such that missed-detection risk is estimated equal or higher than its actual value). Monte Carlo simulations are used to evaluate tightness-of-fit for the case of modeling an input vector with noise described by a two-dimensional, radially symmetric exponential PDF. The results indicate that the spherical bound outperforms the other bounds (where the spherical bound can be computed). Where the spherical bound cannot be computed, the cylindrical bound performs best. The conical bound is consistently the worst performer of the three proposed overbounds. Nonetheless, the conical bound offers one distinct advantage over the other two overbounds, in that it can be expressed in an entirely analytic form.

ACKNOWLEDGEMENTS The author gratefully acknowledges the Federal Aviation Administration GBAS Program (DTFACT-15-C-00033) for supporting this research. The opinions discussed here are those of the author and do not necessarily represent those of the FAA or other affiliated agencies.

REFERENCES [1]

R.G. Brown, "Receiver autonomous integrity monitoring," Global Positioning System Theory and Applications, Progress in Astronautics and Aeronautics, vol. 163, B. Parkinson and J. Spilker, Eds. U.S.: AIAA, 1996, pp. 143-168.

[2]

J. Rife and S. Pullen, “Aviation applications,” GNSS Applications and Methods, S. Gleason and D. Gebre-Egziabher, Eds. Artech House, Norwood, MA: Artech House, 2009, pp. 245-268.

[3]

P. Enge, T. Walter, S. Pullen, C. Kee, Y.-C. Chao, Y.-J. Tsai, “Wide area augmentation of the Global Positioning System,” Proc. of the IEEE, vol. 84, no. 8, pp. 1063-1088, August 1996.

[4]

P. Enge., “Local area augmentation of GPS for precision approach of aircraft,” Proc. of the IEEE, vol. 87, No. 1, 1999, pp. 111-132.

[5]

T. Murphy and T. Imrich, “Implementation and operational use of ground-based augmentation systems (GBASs)—a component of the future air traffic management system,” Proc. of the IEEE, vol. 96, no. 12, pp. 1936-1957, December 2008.

[6]

R.E. Phelts, T. Walter, and P. Enge, “Toward real-time SQM for WAAS: improved detection techniques,” Proc. ION-GNSS, 2003, pp. 2739-2749.

[7]

A. Mitelman, R. Phelts, D. Akos, S. Pullen, P. Enge, “Signal deformations on nominally healthy GPS satellites,” Proc ION NTM, 2004.

[8]

F. Liu, M. Brenner, C.Y. Tang, “Signal deformation monitoring scheme implemented in a prototype local area augmentation system ground installation,” Proc. ION-GNSS, 2006, pp. 367-380.

[9]

T. Houston, F. Liu, and M. Brenner, “Real-time detection of cross-correlation for a precision approach ground based augmentation system,” Proc. ION GNSS, 2011, pp. 3012-3025.

[10] G. Wong, R. Phelts, T. Walter, and P. Enge, “Bounding errors caused by nominal GNSS signal deformations,” Proc. ION GNSS, 2011, pp. 2657-2664. [11] S. Gunawardena and F. van Graas, “Analysis of GPS pseudorange natural biases using a software receiver,” Proc. ION GNSS, 2012, pp. 2141-2149.

25

[12] G. Wong, Y.-H. Chen, R. Phelts, T. Walter, and P. Enge, “Mitigation of nominal signal deformations on dual-frequency WAAS position errors,” Proc. ION GNSS+, 2014, pp. 3129-3147. [13] S. Pullen and J. Rife, “Differential GNSS: Accuracy and Integrity,” GNSS Applications and Methods, S. Gleason and D. GebreEgziabher, Eds. Artech House, Norwood, MA: Artech House, 2009, pp. 87-120. [14] B. DeCleene, “Defining pseudorange integrity – overbounding,” Proc. ION-GPS, 2000, pp. 1916-1924. [15] J. Rife, S. Pullen, B. Pervan, B., and P. Enge, “Paired overbounding for nonideal LAAS and WAAS error distributions,” IEEE Trans. Aerospace and Electronic Systems, Vol. 42, No. 4, 2006, pp. 1386-1395. [16] J. Rife and D. Gebre-Egziabher, “Symmetric overbounding of correlated errors,” NAVIGATION, Vol. 54, No. 2, 2007, pp. 109-124. [17] W.G. Pulford, “A proof of the spherically symmetric overbounding theorem for linear systems,” NAVIGATION, Vol. 55, No. 4, 2008, pp. 283-292. [18] J. Rife, “Overbounding missed-detection probability for a chi-square monitor,” Proc. ION-GNSS, 2012, pp. 1348-1368. [19] J. Rife, “The effect of uncertain covariance on a chi-square integrity monitor,” NAVIGATION, Vol. 60, No. 4, 2013, pp. 291-303. [20] J. Rife and D. Schuldt, “Minimum volume ellipsoid scaled to contain a tangent sphere, with application to RAIM,” Proc. IEEE PLANS, 2014, Monterey, CA. [21] J. Rife, “Comparing performance bounds for chi-square monitors with parameter uncertainty,” IEEE Trans. Aerospace and Electronic Systems, Vol. 51, No. 3, 2015, pp. 2379-2389. [22] J. Rife. “Overbounding false-alarm probability for a chi-square monitor with natural biases,” NAVIGATION, accepted 2016. [23] J. Rife, “Robust chi-square monitor performance with noise covariance of unknown aspect-ratio,” Proc. ION-GNSS+, 2016. [24] R. Braff and C. Shively, “A method of over bounding ground based augmentation system (GBAS) heavy tail error distributions.” Journal of Navigation, Vol. 58, No. 1, 2005, pp. 83-103. [25] J. Warburton, LGF Sigma Pseudorange Ground Establishment, Overbounding and Monitoring: Support Data. Presented by the FAA William J. Hughes Technical Center to RTCA SC-159, December 2002. [26] T. Dautermann, C. Mayer, F. Antreich, A. Konovaltsev, B. Belabbas, and U. Kalberer, “Non-Gaussian error modeling for GBAS integrity assessment,” IEEE Trans. Aerospace and Electronic Systems, Vol. 48, No. 1, 2012, pp. 693-706. [27] P. Ober, D. Imparato, S. Verhagen, C. Tiberius, Christian, H. Veerman, A. Van Kleef, F. Wokke, A. Bos, and A. Mieremet, “Empirical integrity verification of GNSS and SBAS based on the extreme value theory,” NAVIGATION, Vol. 61, No. 1, 2014, pp. 23-38. [28] J. Larson and D. Gebre-Egziabher, “Analysis and utilization of extreme value theory for conservative overbounding," Proc. IEEE/ION PLANS, 2016. [29] Honeywell, Inc., Hazardously Misleading Information (HMI) for GAST-D Signal Deformation Monitor, in progress, August 2016. [30] J. Rife and B. Pervan, “Overbounding revisited: discrete-error distribution modeling for safety-critical GPS navigation,” IEEE Trans. Aerospace and Electronic Systems, Vol. 48, No. 2, 2012, pp. 1537-1551. [31] J. Lee, S. Pullen, and Per Enge. “Sigma overbounding using a position domain method for the Local Area Augmentation of GPS,” IEEE Trans. Aerospace and Electronic Systems Vol. 45, No. 4, 2009, pp. 1262-1274. [32] W.H. Press, Numerical Recipes 3rd edition: The Art of Scientific Computing. Cambridge university press, 2007. [33] L.C. Ludeman, Random Processes: Filtering, Estimation, and Detection. John Wiley & Sons, Inc., 2003. [34] J. Rife, “Spherical approximation to predict performance of a quadratic integrity monitors with non-Gaussian random inputs,” in preparation, 2016.

26