Limits on GNSS Performance at High Latitudes Limits on GNSS Performance at High Latitudes

of a collaborative eﬀort


Introduction
GNSS performance is often described by the Geometric Dilution of Precision, GDOP, a function of the receiver's and satellites' geometry.While users often think of satellites completely covering the sky, the satellites' orbital planes control where they might occur; even when averaged over time and longitude, at all latitudes there are some portions of the sky which never contain satellites.As an example, Figure 1 shows a sky view of the GPS satellite trajectories for a 24-hour period here in Reston VA (the site of ION ITM 2018, 38.9581 • N, 77.3563 • W ) showing a significant no-satellite zone (or "bald spot") in the northern sky.This characteristic is especially significant at higher latitudes, such as might be experienced by vehicles on polar transit routes in these times of climate change.For example, for users near the poles the observed satellites never achieve elevation above approximately 45 degrees.Since the vertical accuracy is heavily determined by the contribution of high elevation satellites, performance is thus limited.To demonstrate this numerically, Figure 2 shows the resulting GDOP for the best 7 GPS satellite position solution for a receiver at different latitudes (5 degree increments), averaged over this same day.It appears that once the latitude exceeds about 55 • the performance decreases considerable (the GDOP increases).
In this paper we characterize this loss of performance with increasing latitude through the development of a lower bound to GDOP characterized as a function of the number of satellites and the receiver's latitude.This paper proceeds as follows: • First, the GDOP performance criterion is briefly reviewed.
• Next, a mathematical characterization of the no-satellite zone, the bald spot, is presented parameterized by the receiver's latitude.
• The GDOP bound from [1] is reviewed, citing the impact of the no-satellite zone on it.
• The GDOP bound is modified to take into account the no-satellite zone.
• Finally, some real data is presented showing the utility of the analytical results.

GDOP
GNSS receivers convert the measured satellite pseudoranges into estimates of the position and clock offset of the receiver.A common implementation of the solution algorithm is an iterative, linearized least squares method.Assuming that pseudoranges from 4 or more non-coplanar satellites are measured, the direction cosines matrix G is formed and used to solve an over-determined set of equations to solve for position and receiver time offset.Using a local East (e), North (n), and Up (u) coordinate frame, for m satellites this matrix is in which (e k , n k , u k ) is the unit vector pointing toward the k th satellite from the receiver's location.
Since the pseudoranges themselves are noisy, the resulting estimates of position and time are random variables.To describe the accuracy of this solution, it is common to describe it statistically via the error covariance matrix, equal to G T G −1 scaled by the User Range Error, URE [2].Rather than considering the individual elements of this covariance matrix, users frequently reduce it to a scalar performance indicator.The most common of these is the Geometric Dilution of Precision, or GDOP, defined as In words this is the square root of the sum of the variances of the four solution terms (three position variables and the time offset) without the URE scaling, a measure of the size of the error ellipsoid.Other common reductions of the covariance matrix are the Position (PDOP), Horizontal (HDOP), Vertical (VDOP), and Time (TDOP) portions of GDOP.
It is well known that the GDOP is a function of the satellite geometry; with only a few visible satellites in poor sky locations the GDOP can become quite large.However, for a future with multiple, fully occupied GNSS constellations, it is expected that receivers would select those satellites to track so as to achieve the best possible performance; see, for example, [3,4].Hence, we think that an understanding of both how small the GDOP can be (a lower bound) and the characteristics of the constellations that meet that bound are of value in the satellite selection process and, generally, in understanding GNSS performance.Investigating the best possible GNSS satellite constellation with respect to GDOP is not a new problem: • It is known that the GDOP is a monotonically decreasing function of the set of satellites employed in the solution [5].Specifically, adding an additional satellite reduces the GDOP independent of the sky location of this extra satellite, even if it lies directly atop one of the original satellites.This monotonicity extends to multiple constellations except when the new satellite is the first of its constellation [6].
• Minimum GDOP for 4 satellites, with reference to optimizing the tetrahedron formed by their locations, has been considered by multiple authors, see e.g.[7,8].
• In a recent paper [1] these authors were able to develop an achievable lower bound to GDOP for terrestrial applications (i.e.satellites restricted to be above the horizon).Letting m represent the number of satellites, this bound is The constellations that achieve this bound consist of approximately 29% of the satellites at zenith and the remaining 71% "balanced" around the horizon (balance having a specific mathematical definition, see also [9]).
In that same work the bounds were extended to PDOP, to non-zero mask angle, and to multiple GNSS constellations.In [9] we showed, by example using real ephemeris data, that good GDOP/PDOP performance resulted from constellations similar to the "best" constellations resulting from the bounds.Recently these bounds were extended to allow for an additional range measurement (e.g. an altimeter) [10].
To date the results on GNSS performance and the development of satellite selection algorithms have NOT taken into account where on the Earth the user is.While not often discussed, the orbital planes of the satellites limit those parts of the sky in which satellites might be visible.For example, users near either of the poles see no satellites at elevation above about 45 degrees.Combining this with a mask angle of 5 or 10 degrees leaves a narrow ring of satellite elevations, limiting vertical accuracy.Even at lower latitudes some portions of the sky are off limits: • At the equator no satellites are visible at the horizon for ±30 degrees of azimuth around each pole; fortunately, this no-satellite swath only has maximum elevation of about 15 degrees.
• At latitude 40 degrees north (near the conference site in Reston) the no-satellite portion of the sky rises from the north pole reaching an elevation of 60 degrees and azimuth of −40 to 40 degrees.
Clearly these regions in the sky where no satellite could possibly be situated will impact possible system performance.

Characterizing the Bald Spot
The first step in understanding the impact of the bald spot in the GNSS constellation is to characterize its size and location as a function of the receiver's latitude.Our approach is to fix the location of the receiver on the Earth at latitude α, to define the position of a typical satellite on its orbit, and to rotate the earth, tracking the satellite's position in the sky relative to the receiver (in coordinates of azimuth and elevation) so as to find the boundary of the bald spot.For mathematical simplicity, below we set the receiver's location at a convenient spot on earth and rotate the orbit around the earth instead of rotating the location underneath the orbit.
1.The Satellite's Location: Within a coordinate system defined by the orbital plane (assumed to go through the center of the earth), and assuming zero eccentricity, the satellite's positions are in which r s is the orbital radius and ν ranges over a full 360 degrees.Transforming to a coordinate frame with origin at the earth's center, inclination angle θ, and right ascension φ requires two rotation matrices: The idea is that we can trace out a satellites location by keeping θ fixed while ν and φ vary over 360 degrees (ν to trace out the satellite's movement in its orbit and φ to account for the earth's rotation).

The Receiver's Location:
Our goal is to find a satellite's elevation and azimuth for a point on the earth at a fixed latitude, α.Since we will allow the right ascension to vary, we can without loss of generality fix the longitude of this earth location so that the y component is zero; hence, its location in x, y, z is  By letting ν and φ both range over a full 360 degrees, the expressions above provide the possible elevation and azimuth angles for a satellite on a circular orbit of radius r s and inclination angle θ (assumed constant for all satellites) as seen by a receiver at latitude α on the surface of the Earth (for points h units above the spheroid we could add h to r e ). Figure 3 shows the results for a GPS satellite (r s = 26, 559.7 km, θ = 55 • ) and a receiver at the surface of the earth with various latitudes from the equator to the north pole.For these graphics the resolution on ν and φ was 5 • .As expected, receivers at the equator do not see satellites at the horizon directly to the north or south.With increasing latitude, the bald spot rises from the north pole, eventually becoming directly overhead.
Examining Figure 3 we make several points: • Satellite azimuths and elevations appear to have a continuous range except for one or two generally oval areas -this can be easily imagined from the figure and becomes more obvious with denser computation over ν and φ.
• Not obvious from the figures, but notable in computation, is that the edge of the northern bald spot results from ν = 90 • (the right image in Figure ?? repeats the plot, but with the ν = 90 • results plotted as blue crosses).Similarly, ν = 270 • generates the boundary of the southern bald spot, if it exists.
• This images are for a receiver at sea level; as the receiver altitude increases the bald spot grows.
Based on these comments, especially the second, we focus below on the mathematics for ν = 90 • .
which partially describes the extent (height) of the bald spot in the north/south direction.Figure 5 plots this result versus latitude α.The upper blue line is the elevation of the leading edge of the bald spot (using the Manipulating the algebra and trigonometry, this is equivalent to α = θ; the leading edge hits zenith when the latitude is equal to the satellite's inclination angle.
Similarly, the bald spot breaks off from the horizon when the trailing elevation equals 0 • , or

A Brief Review of [1]
Following [1] the GDOP under GNSS alone is determined by the elements of the symmetric matrix G T G with Finding the minimum of GDOP in [1] was a four step procedure: 1. Set the off-diagonal terms to zero.Strictly we cannot do this for all of these terms; however, setting the ones shown in red (all but the k u k since each u k ≥ 0) is desired.Hence, we want Effectively, these are conditions that the satellite constellation exhibit symmetry in multiple directions.In [1] these conditions, along with next one below, are used to define a "balanced" satellite constellation.See [9] for a more detailed discussion of balance.
2. Set the first two diagonal terms (in blue) to be equal to one another We note that it is possible to satisfy all 6 of these conditions for any m greater than 3 if the sky locations are unrestricted (beyond being above the horizon) [9].
3. The impact of the remaining off-diagonal term on GDOP (in black) is minimized by forcing all of the satellites to be at either zenith or on the horizon (elevations 90 • and 0 • , respectively); let β represent the fraction of the m satellites at zenith, then the GDOP bound reduces to 4. Finally, standard calculus provides the optimum value for β is ≈ 0.29 which yields the GNSS-only bound presented above in Eq. ( 1).
Including the bald spot we expect a similar result, specifically, of the form in which the numerator constant depends upon the latitude.

The Bald Spot and the GDOP Bound
To extend the results in [1] for restrictions on satellite locations we consider several cases, depending upon latitude: • α < 55 • : First, consider the case of smaller latitudes, α < θ.In this case the "bald spot" does not include zenith, so the results from [1] hold; the small portion(s) of the horizon restricted due to the no satellite zone does not inhibit the horizon satellites from forming a balanced constellation.Specifically, Eq. ( 1) holds so that C(α) = 3.45 and the performance is the blue part of the plot of C(α) = √ m × GDOP for this range of α shown in Figure 7.
• At the North Pole, α = 90 • : Arguments presented in [11] prove that the best GPS satellite constellation for a receiver at the pole is to have satellites limited to either the horizon or elevation θ p = 45.The blue line is the original bound for α < 55 • ; the red dotted line is a lower bound to the lower bound which clearly indicates that the performance gets worse (monotonically) as the no-satellite zone encompasses more of the high elevation region; the red dot is an achievable bound at the pole; the red line is achievable performance for the conjectured solution.
Assume that βm of the satellites are at this higher elevation, (1 − β)m are at the horizon, and that both sets are "balanced" as per [1].Then from [11] we have that β should be a root of the 4 th order polynomial Rooting the polynomial yields β = 0.41 (i.e.41% of the satellites should be at the highest elevation) and GDOP (90 a 12% increase in the GDOP coefficient.This appears as the red dot at α = 90 • in Figure 7.
• Latitude between 55 • and 90 • -Part 1: A non-achievable lower bound can be constructed by shrinking the no-satellite region and applying the results from [11].Specifically, consider Figure 8 which shows the situation at latitude 65 • ; the pink region is the no-satellite zone.Shrinking this to the circular cap about zenith (shown in red), we can compute the size of this cap using Eq. ( 2) and apply the results of [11].The result is the red dotted line in Figure 7.As an optimum satellite constellation for this smaller region would required balanced satellites in the pink region, the bound is generally unachievable.
• Latitude between 55 • and 90 To tighten the bound, consider Figures 9 and 10 (in these and other sky view figures below, the circle and diamond symbols are meant to show locations for satellites, not individual satellites; for example, the leftmost sky view in Figure 9 could be showing the locations of a 4 satellite constellation, but might also be showing locations for 8 satellites, two at each spot, the center sky view would be optimum for 6 satellites, 2 at zenith and 1 at each of the other 4 horizon locations): -The first of these shows possible optimum configurations at latitude 55 • ; 29% of the satellites would be at zenith while the remaining 61% are distributed evenly around the horizon (potentially at either 3 or 4 locations).
-The second of these shows possible optimum configurations at latitude 90 • ; now 41% of the satellites are distributed evenly at elevation 45.3 • (2 or 3 locations) while the remaining 59% are distributed around the horizon (evenly for the first two figures; in the last configuration with only two locations at elevation 45.3 • , north and south, there would need to be a higher percentage of satellites to east and west to balance them).Note that there are many other configurations that would achieve the same minimum GDOP at either of these latitudes (for example, arbitrary rotations of the patterns or more symmetric locations beyond three or four); these choices just help to motivate the constellations considered below.Specifically, as shown in Figure 11 we started with a generic configuration for latitudes between 55 • and 90 • : -Some fraction of the satellites are distributed amongst 4 locations on the boundary of the bald spot, due north, due south, and a pair symmetrically east/west; the remaining satellites appear at 6 locations on the horizon, again divided due north, due south, and to two symmetric east/west pairs of locations.
-The swaths of pink and blue in this diagram are meant to suggest that the azimuths of the east/west pairs (pairs of diamond symbols both on the horizon and the boundary) will be allowed to range over [0 • , ±180 • ] during an optimization step.
-Three azimuths (for the locations on the boundary the azimuth specifies the elevation, for the locations on the horizon the elevations are identically zero).
Constrained numerical minimization of the GDOP over these 10 parameters yields the performance shown as the solid red line in Figure 7; we conjecture that no better configuration of satellites exists.Further, an examination of the optimum choices for the 10 parameters indicates that the constellation is over parameterized.Specifically the fraction of satellites at horizon north goes to zero and the three azimuths do not change from their initialized values over the course of the optimization.Other (reasonable) initial values for the three azimuths also yield the same level of performance with different results for the optimized fractions.
This generic constellation can be altered in various ways without impacting the optimum performance.The simplest is shown as an inset in Figure 12: two locations on the bald spot boundary (due north and due south) and three locations on the horizon (due south and a symmetric east/west pair).Figure 12 itself shows the fractions of satellites at these locations, the color and symbol choices match the inset diagram; Figure 13 shows the azimuths: -The locations on the bald spot boundary starts with approximately 29% of the satellites due south (blue circles) with 0 due north (red circles) at latitude 55 • and eventually converge to a combined total of 41% on the boundary at latitude 90 • ..
-Interestingly, we don't need any satellites due north on the bald spot boundary (red circles) until we exceed approximately 70 • of latitude.
-At latitude 55 • the 71% of the satellites on the horizon are split evenly across the three locations and these locations are evenly distributed on the horizon (120 • spacing); as the latitude increases the optimum values for the fractions and the azimuths change somewhat.
In a similar fashion we optimized the performance of several other reductions of the generic configuration in Figure 11.While these yielded different sets of optimum parameters, we repeat that they all resulted in equivalent GDOP performance !! In other words, there are many possible equivalent minima for GDOP (recall that for latitudes below 55 • all that matters is that we put 29% of the satellites at zenith and the remaining 71% "balanced" around the horizon -these could be scattered into 3 locations, 4 locations, or more (see [9]).
Examples of this, for latitudes 60

Examples
To compare the developed bounds and conjectured optimum constellations to reality we generated GPS constellation data for a 24-hour period at multiple locations on the Earth (latitudes 60 • , 70 • , and 80 • N, 5 • spacing in longitude) and computed the GDOP for a 7 satellite solution, keeping the best configurations.Figure 17 shows an example from each of these latitudes that mimic the optimum configurations of the above discussion.Specifically, we notice how three (or four) satellites hug the no-satellite boundary and that the others are near the horizon in a somewhat balanced configuration.

Conclusions/Future
We have developed a lower bound to GDOP taking into account the latitude of the receiver and the fact that at higher latitudes the satellites are restricted to be lower in the sky.For parts of the range of latitudes (specifically, 0 to 55 • and at 90 • ) the bound is tight (i.e.achievable).For the remainder of the range we have a non-achievable bound and conjecture a second result to be the actual lower bound (both shown in Figure 7).
The lack of high satellites at higher elevation directly impacts vertical accuracy.As an example, Figure 18 repeats the GDOP coefficients of Figure 7, but separates out the horizontal, vertical, and time portions of the DOP (the GDOP is the square root of the sum of squares of these components).Quite clear from this figure is that the growth in GDOP for higher latitudes is primarily due to a loss of VDOP performance (the VDOP increases for latitudes above 55 • while the HDOP and TDOP remain relatively constant).

Figure 1 :Figure 2 :
Figure 1: Traces of the GPS satellites' orbits over the conference site for a 24-hour period in October 2017.
r e is the earth's radius.(Note that we could more correctly employ an oblate spheroid model for the Earth, yieldingx r =   a cos α 0 b sin α  with a and b being the equatorial and polar radii, respectively, but observe that there is no observable numerical difference in the results below.)Thisreceiver location forms a tangent plane on the surface of the earth.The unit vector this plane and is used below in computing elevation.The unit vector a = tangent plane and points north and is used below in computing azimuth.Finally, the unit vector b = orthogonal coordinate frame with n and a.

3 .
The Sky View's Terms: The difference vector (pointing from the the receiver to the satellite) isd = x i − x r =  r s (cos φ cos ν − sin φ cos θ sin ν) − r e cos α r s (sin φ cos ν + cos φ cos θ sin ν) r s sin θ sin ν − r e sin α   Let d represent the length of this vector d = |x i − x r | and d n , d a , and d b be its projections on the coordinate frame defined by the receiver's location d n = d • n , d a = d • a and d b = d • b The elevation and azimuth angles to the satellite are trigonometric functions of these projections elev = sin −1 d n d and azi = tan −1 d b d a 4. The Bald Spot:

Figure 3 :
Figure 3: Range of visible satellites versus latitude.

Figure 5 :
Figure 5: The extent of the northern bald spot's boundary at azimuths 0 • and 180 • .

rWe can solve for sin φ sin φ = r s sin α sin θ − r e r s cos α cos θ We also have cos φ = 1 − sin 2 φ 1   cos θ r 2 s cos 2 α cos 2 θ − r 2 s sin 2 α sin 2 θ + 2r s r e sin α sin θ − r 2 eFigure 6 Figure 6 :
Figure 6  plots this result versus latitude α.For lower latitudes the the intersection with the horizon starts at 30.3 • for a receiver at the equator, increases and then decreases smoothly with increasing latitude, finally disappearing (shown here as azimuth equal to zero) for latitudes above 53.4• .

3 •Figure 7 :
Figure 7: The coefficient of the GDOP, C(α) =√ m × GDOP , as a function of latitude α.The blue line is the original bound for α < 55 • ; the red dotted line is a lower bound to the lower bound which clearly indicates that the performance gets worse (monotonically) as the no-satellite zone encompasses more of the high elevation region; the red dot is an achievable bound at the pole; the red line is achievable performance for the conjectured solution.

Figure 8 :
Figure 8: The region for a lower bound to the lower bound.

Figure 12 :Figure 13 :
Figure 12: Allocation of satellites (as percentages) for the 5 location constellation versus latitude.
, and 16 (each figure shows three distinct constellations with the locations and fractions of satellites marked, the leftmost being the simplest configuration).

Figure 18 :
Figure 18: Separation of the GDOP into horizontal, vertical, and time portions vs latitude.