Rethinking Star Selection in Celestial Navigation

Abstract

In celestial navigation the altitude (elevation) angles to multiple celestial bodies are measured; these measurements are then used to compute the position of the user on the surface of the Earth. Methods described in the literature include the classical “altitude-intercept” algorithm as well as direct and iterative least-squares solutions for over determined situations. While it seems rather obvious that the user should select bright stars scattered across the sky, there appears to be no established results on the level of performance that is achievable based upon the number of stars sighted nor how the “best” set of stars might be selected from those visible. This paper addresses both of these issues by examining the performance of celestial navigation noting its similarity to the performance of GNSS systems; specifically, modern results on GDOP for GNSS are adapted to this classical celestial navigation problem. INTRODUCTION In the world of GNSS, position accuracy is often described by the geometric dilution of precision, or GDOP [1]. This measure, a function of GNSS constellation geometry (specifically the azimuths and elevation angles to the satellites employed in the solution), is a condensed version of the covariance matrix of the errors in the position and time estimates. Combining the GDOP value and an estimate of the user range error allows one to establish the 95% confidence ellipsoid of position. Many papers in the navigation literature have considered GDOP: early investigations computed the GDOP as a function of time and location (on the surface of the Earth) to show the quality of the position performance achievable (e.g. [1, p.211]); the mathematics of the GDOP measure have been considered (e.g. [2]); and lower bounds on GDOP have been established as a function of number of satellites (e.g. [3]). A related problem is satellite selection for best GDOP. Specifically, in these days of multiple, full constellations, there could be times when there are many more satellites available than the receiver needs or is able to handle (e.g. [4], where GBAS channel bandwidth limitation was an issue); hence, some sort of satellite selection process is used, selecting those satellites that provide good (or best) performance. For example, early GPS receivers employed only 4 satellites to provide a solution while 6 or more satellites were typically available; multiple authors considered the best way to choose those 4 (e.g. [5]). Since then, many papers in the navigation literature have considered methods to find the best, or nearly best, subset of satellites both from single and multiple constellations (e.g. [6–14]). It is also possible to extend these approaches to horizontal DOP (HDOP) for navigating on the surface of the Earth. For 2-D positioning on the surface of the earth, one can also consider measuring altitudes of celestial bodies (elevation angles, relative to an artificial or visible horizon, measured with a sextant), to locate the user [15–18]. Assuming multiple measurements of star elevations, solving for position can be cast as an overdetermined set of equations; hence, it is natural to employ a least squares solution. Whatever the position units, latitude and longitude or East and North in a local coordinate frame, these equations are nonlinear so can be solved via an iterative, linearized approach (just like in GNSS) [19, 20]. Presumably a user is interested in the “most accurate” position estimate, which would suggest that the geometry of the celestial bodies chosen for measurement is important. In fact, guidance provided in Bowditch [15, p.276] is to “. . . Predict expected altitudes and azimuths for up to eight bodies when preparing to take celestial sights. Choose the stars and planets that give the best bearing spread. Try to select bodies with a predicted altitude between 30 degrees and 70 degrees.” (Underline added for emphasis.) While this guidance seems reasonable, it is a bit ambiguous as to what constitutes a “best bearing spread,” and the general interpretation has been to choose visible celestial bodies that are fairly evenly spread around the horizon in azimuth. As we show below, this is a sufficient condition for the best geometry, but it is not a necessary condition, so can be overly restrictive. The contents of this paper are: • We very briefly review the fundamental methods of celestial navigation, highlight a mathematically precise definition of HDOP for the celestial navigation problem similar to that used in GNSS applications, and note the similarity of the classical altitude-intercept method to current GNSS positioning algorithms. • We minimize this HDOP expression over the choice of the stars’ locations; the results are a set of necessary conditions on the stars’ azimuths and an achievable lower bound to HDOP as a function of the number of stars employed (parameterized by m). These conditions demonstrate the limitations of the conventional wisdom on star selection (and, hence, motivate the title of this paper). Specifically, these conditions do not, in fact, require that the m celestial observations be separated by 360/m degrees of azimuth in order to minimize HDOP. No longer requiring an even spread around the azimuth circle provides added flexibility to a navigator, especially if there is lack of a visible horizon in some direction (due to land, fog, or lighting conditions). • We present a simple star selection algorithm for HDOP minimization. Of the roughly 9000 stars visible to the naked eye, 58 stars (including Polaris) plus 4 planets (Venus, Mars, Jupiter, Saturn), plus the Sun and Moon, have been selected for tabulation in each year’s production of the Nautical Almanac [21]. These were chosen based upon visibility and ease of identification. Imagine if approximately one-half of these were visible (say 30); an exhaustive search for those 8 celestial objects that minimize HDOP is a “30 choose 8” problem, involving an evaluation of 5.8 million possible subsets of objects for minimum HDOP. Exploiting our knowledge of the HDOP formulation, we present a real time algorithm with significantly reduced complexity for selecting a near-optimum subset of stars for position determination. • Finally, we include a real-world example, taken from tabulated visible celestial bodies in the 2019 Nautical Almanac, and show the utility of our results. We highlight that the celestial bodies selected in the optimal solution are, in fact, not necessarily equally spread in azimuth. REVIEW OF CELESTIAL NAVIGATION Celestial Navigation to a mariner is considered to be the art of calculating a 2-D position on the surface of the earth with the aid of the sun, moon, planets, and major stars. For aviators it typically means calculating that 2-D position relative to a known (barometric) altitude above the earth. In either case it usually means the navigator is obtaining lines of position (LOPs) by measuring the altitudes of celestial bodies relative to an actual or artificial horizon, typically by using a sextant. The general strategy in celestial navigation is to (1) choose a set of reasonable celestial bodies for measurement, based on visibility, azimuth, and altitude, (2) measure observed altitudes of that set of celestial bodies at known times, (3) correct altitude measurements for effects such as refraction, height of eye, or parallax, and (4) determine each celestial body’s location on the celestial sphere at the observation time (e.g. declination and Greenwich Hour Angle). Each celestial body’s altitude measurement creates a “circle of equal altitude” on the surface of the earth. The intersection of these circular lines of position (LOPs) results in a unique celestial fix. For a more local “small-scale” position solution, and with observed altitudes below 75 degrees, one can draw these “circular arcs” as actual lines of position, and this is what is typically done in maritime navigation. The classical method used to solve for position at sea is called the “altitude-intercept” method, also known as the Marcq St. Hilaire method (after the French navigator who developed the technique). In this method we use known celestial body locations plus an assumed position (AP), presumably reasonably close to the actual position, to calculate the computed altitudes one would expect at known times at that AP. This can be accomplished using the Nautical Almanac and Publication 229 or https://aa.usno.navy.mil/data/docs/celnavtable.php to determine both the computed altitudes and azimuth angles. The difference between the computed and observed altitudes (measured in minutes in arc), combined with the azimuth, is used to plot each celestial line of position (LOP). More specifically, the line of position produced is oriented perpendicular to the azimuth toward each celestial body, and is initially drawn through the AP. The final LOP is shifted towards or away from the celestial body in the amount of the arc difference between the computed and observed altitudes (recognizing that one minute of arc is approximately one nautical mile). If the computed altitude is greater than the observed altitude, we shift the LOP away from the celestial body, along the azimuth line; if the computed altitude is less than observed altitude, we shift the LOP towards the celestial body, along the azimuth line. This is done for each celestial body, and the final fix is computed from the intersection of all LOPs. If the lines did not intersect in a single point, the “estimation of the ship’s position from the somewhat chaotic image of a number of position lines (is) left to the ‘keen judgement’ of the observer” [20]. The lack of a common intersection in the overdetermined case is familiar to those in the GNSS community; the usual response is to reformulate the problem in a least squares sense, minimizing the residual error for each LOP. Since the equations are nonlinear, this involves choosing an initial assumed position, linearizing th