## Abstract

Planetary rover navigation frequently relies on dead reckoning and external infrastructure such as orbiting satellites. Celestial navigation techniques combine measurements of the Sun, stars, and gravity to provide autonomous absolute localization. This study examines the performance of digital star sextants (DSS)—a suite of sensors combining a star tracker and an inclinometer—on estimating position on the planetary surface. In particular, we discuss the estimation, calibration, and error analysis for an elevation-only measurement formulation that does not rely on ground-truth heading information. Field tests and Monte Carlo simulations provide validation of the proposed techniques. The real-world performance of the experimental system gives a mean single-orientation error of 296 m. The relative agreement between the predicted and observed error reveals a clear roadmap to help evaluate the impact of prospective sensor improvements on DSS performance.

## 1 INTRODUCTION

Most planetary rover missions rely, at least in part, on dead reckoning to navigate. By their nature, these methods accumulate errors over time resulting in position estimate drift. On Mars, this has been mitigated using orbital landmark tracking to occasionally provide absolute state estimation by identifying and matching landmarks visible to both the rover and an orbiting satellite (Li et al., 2004, 2005, 2007). Star-based celestial navigation provides an alternative method for localization that does not require orbiting infrastructure or extensive image databases.

The Mars Pathfinder missions used images of the Sun to measure heading (Eisenman et al., 2002) but instantaneous observations of the Sun are limited in their accuracy and the information they provide. In contrast, star trackers typically offer improved measurement accuracy and complete inertial attitude knowledge (Wertz et al., 2011). Inclinometers relate this attitude measurement to the local-horizontal frame permitting the direct calculation of sensor position (Enright et al., 2012). We refer to this particular system as a digital star sextant (DSS). This paper presents a suite of estimation equations based on a deterministic heading-independent forward measurement model for use with a DSS.

The recent study of DSS systems in literature has revolved around their calibration. Recent works in this field include that of Wei et al. (2019), Zhan et al. (2020), and Jovanovic and Enright (2017). These authors examined DSS calibration taking particular concern with accurate longitude and latitude estimates, but not heading. A DSS is capable of providing all three, but because the precision of star tracker heading measurements is often much better than available ground-truth, it is challenging to assess the heading errors in the calibration process. Ning and Fang (2009) conducted simulated trials under this principle which, along with Zhan et al. (2020), recognized that heading was unnecessary in the general position solution. Ning and Fang (2009) developed a measurement model that converted global position into star tracker and inclinometer measurements (which we refer to as a backwards measurement model) that was independent of heading for use in their Kalman filter. Zhan et al. (2020) developed a heading-independent measurement model based on star elevation that required a-priori information about global position. By contrast, the measurement model we present in this paper is deterministic.

In this paper, we present a DSS measurement model that also utilizes star elevation measurements to achieve heading independence. Together with this model, we include the derivation of a least-squares localization solution, a covariance model, and a calibration framework. Finally, we validate our models using field and simulated data. Although our results do not directly improve upon the accuracy reported by other sextant studies, our approach offers a more rigorous analysis of the calibration process and system error characteristics. Our covariance models directly relate the DSS performance to the sensor specifications and present a clear roadmap for improving DSS accuracy.

## 2 MATHEMATICAL FRAMEWORK

This section presents the governing mathematical framework for our localization solution. We include a general discussion of the celestial localization problem and measurement models for the DSS component sensors. The sensor models then allow us to derive our elevation-only localization estimate and a corresponding parametric model of system covariance. Due to the sensitivity of the DSS output to even small angular errors, our system models depend on a variety of calibration parameters for proper operation. We conclude this section with analysis and discussion of parameter observability and present a strategy for effective calibration.

### 2.1 General Celestial Localization

The primary components of our DSS system are a star tracker (used to measure and identify stars) and an inclinometer (used to measure gravity direction). Calculating planetary orientation corresponding to each measurement is also essential—for instance, the Earth’s rotation amounts to approximately 465 m/s on the equator—so a precise time reference is also needed. Although a deployed DSS system would require a standalone precision time reference, an onboard GPS receiver is an acceptable timing source for our experiments. Timekeeping on any extraterrestrial mission would not be possible using GPS, however, meter-level position accuracy requires relatively coarse millisecond-level clock accuracy. This synchronization requirement is not considered particularly challenging over standard radio links (Mills, 2016).

The localization process relies on relating quantities between a number of reference frames. The frames relevant to DSS operations are shown in Figure 1. Between important pairs of frames are notes explaining the corresponding rotational relationship. Beginning with star tracker frame , the rotation **C*** _{IS}* is derived from the star tracker attitude measurement, where is the inertial frame. The rotation from to the planetary fixed frame, , represents the rotation of the planetary body. For this paper we calculate

**C**

*from the observation time stamp (recorded as a UTC time), and the Earth Orientation Parameters provided by the International Earth Rotation and Reference Systems Service (IERS) using the IAU 2000/2006 reduction (Petit & Luzum, 2010). From to the topocentric frame, , requires knowledge of global position, which is the goal of the DSS. This solution can be obtained by solving for the remaining rotations.*

_{FI}Going from the star tracker’s frame, , to the inclinometer frame requires the mounting rotations, **C*** _{SA}*. The rotation between and relies on knowledge of the DSS roll and pitch angles as well as the sensor mounting geometry. The former are obtained from a gravity vector calculated from inclinometer measurements. The frame, shown in Figure 2, is defined as North-East-Down (NED) in which the

*x*-axis always points north, the

*y*-axis points east, and the -axis points down. The frame does not move relative to the frame if the DSS is stationary on the planet’s surface.

With these transformations defined, the rotation that defines global position can be expressed as:

1

This same rotation can also be expressed directly in terms of longitude *λ*, geodetic latitude *ϕ*, and heading *ψ*. Adopting the North-East-Down (NED) definition for the topocentric frame and using *s*_{(⋅)} and *c*_{(⋅)} as shorthand notations for sine and cosine, we have:

2

Equating this result to Equation (1) allows for the recovery of position and heading estimates from the elements *c _{i,j}* of

**C**

*:*

_{FN}3

### 2.2 Star Tracker Measurement Model

At a top level, the star tracker’s processing sequence consists of turning imaged stars on the detector plane into star vectors, , and then matching those star vectors to a catalog of inertial star vectors, **s*** _{I}*. From corresponding pairs of star vectors, the sensor computes its own inertial orientation, . The small error rotation between the measured orientation and true orientation is expressed as a three-element angular error vector,

**. Thus the approximate expression for the error rotation is:**

*δθ*4

In addition to the inertial attitude, **C*** _{SI}*, the star tracker provides the sensor-frame star vector measurements, , and the corresponding catalog identifiers for each observed star. When used within the atmosphere, these star vectors are refracted and, hence, are regarded as the apparent vectors.

To obtain airless star vectors on Earth, , atmospheric refraction corrections are calculated according to Bennett (2009) and applied to the vectors using the method outlined in Enright et al. (2012). The airless star vectors still contain measurement errors not associated with refraction effects.

### 2.3 Inclinometer Measurement Model

The inclinometer measurement model produces a rotation, **C*** _{AN}*, between the inclinometer and navigation frames, the latter of which is always level by definition. Biaxial inclinometers like the one used in this study are commonly constructed from two orthogonally mounted single-axis inclinometer subsensors (see Figure 3). These subsensors are nominally mounted with a 90° -rotation between them. We refer to the individual subsensor frames as and for the

*x*and

*y*axes, respectively. We define these frames as nominally aligned to one another, not rotated by 90°, as well as the inclinometer’s body frame, . The subsensor frames and are identified individually and may be misaligned relative to .

The navigation frame is defined with the -axis pointed down, however, many manufacturers define inclinometer sensor frames with the -axis in the nominal *up* direction. Introducing as a virtual frame representing a level inclinometer, we can account for the -axis change and any inclinometer mounting misalignment with a constant **C*** _{LN}* matrix:

5

Our choice of model is based on Jovanovic and Enright (2020) in which the gravity vector in is:

6

Each subsensor defines a sensing plane with normals constructed according to:

7

where *θ _{x}* and

*θ*are the inclinometer measurements. Note that

_{y}**n**

*and*

_{x}**n**

_{y}are defined in and , respectively. Since these frames, along with , are nominally aligned, no transformation is necessary in this derivation.

Expanding the cross product gives the inclinometer measurement model as:

8

We then use the model in Equation (8) to produce the rotation between and . If we define our rotation **C*** _{AL}* as:

9

then:

10

By equating Equation (8) to Equation (10), we arrive at:

11

### 2.4 Elevation-Only Measurement Formulation

We now present a novel heading-independent forward measurement model for the DSS. This localization solution relies on a least-squares minimization of the errors between observed and predicted elevation angles. The elevation angles are computed from the star and gravity vectors when expressed in a common reference frame. These angles are independent of the working frame of reference.

Expressed in a given frame, a star’s elevation, *ε _{i}*, is a function of the angle between the gravity and star vectors:

12

Thus the airless elevation can be computed directly in using the DSS measurements of star and gravity vectors. In , the corresponding inertial star vectors are corrected for stellar aberration^{1} and then rotated into . This gives the true (geometric) vector **s*** _{F,i}*. When working in this frame, the gravity vector,

**g**

*, encodes global position:*

_{F}13

To simplify the navigation solution, we replace the explicit elevation calculations with the modified error metric *e _{i}* based on the

**s**

^{T}

**g**terms alone:

14

This also allows us to solve for **g*** _{F}* rather than

*λ*and

*ϕ*explicitly.

Although each error component compares a single star vector with the gravity vector measurement, calculating position requires a batch of star measurements. Normally, these batches would reflect the multiple stars identified in a single star tracker image, but they can also include the vectors collated from multiple exposures collected at the same location.

An array of *N* star vectors can be compactly represented as:

15

which, when substituted into Equation (14), gives:

16

Solving for position now requires us to calculate the optimal **g*** _{F}* vector. As this is a unit vector, two stars would be sufficient to solve the component parts of this set of equations

^{2}. However, we will focus on a general case of three or more stars which we can solve as a least-squares problem. We can define a least-squares cost function

**g**

*and introduce a Lagrange multiplier, , to enforce a unity constraint on*

_{F}**g**

*:*

_{F}17

18

where:

19

Taking the derivative with respect to **g*** _{F}* and setting the left-hand side to zero gives:

20

21

where:

22

and is the optimal gravity vector.

This is a ridge regression problem that must be solved iteratively by varying to find a unity magnitude **g*** _{F}* vector. In many cases, the optimal Lagrange multiplier is zero or very close to it and converges after only a few iterations. On occasion, our implemented solver converges to a local minima in the solution space. Therefore, in practice, we simply set the Lagrange multiplier to zero to emphasize reliability over accuracy.

It is important to ensure that is invertible. If is set to zero, only needs to be positive definite for this to be the case.

23

Conversely, were positive definiteness not satisfied, then we must have:

24

for some nontrivial **a**. Recall that the rows of **S*** _{F}* are the star vectors and a minimum number of three star vectors is implied. This means that for Equation (24) to be satisfied,

**a**would need to be orthogonal to all the vectors simultaneously, something which is only possible if all the star vectors are coplanar. Therefore, is invertible if the constituent star vectors are not coplanar.

Combining the results of this section, we can say that the vector , which encodes global position, can be calculated from measurements according to:

25

### 2.5 Elevation-Only Covariance

To complement the position solution presented in Section 2.4, this section examines the error characteristics of the position estimates. Our goal with this derivation is to obtain the covariance **P** of **g*** in terms of the measurement error characteristics. For this derivation, we assume that both the star tracker and inclinometer measurements are corrupted by independent error components confined to small regions tangent to the measured vectors. Thus:

26

27

We also apply this concept to the estimated gravity vectors and *γ* values:

28

29

where ** δη, δγ, δg_{A}**, and

**are the error components.**

*δ*S_{S}The Equation (19) and Equation (21) relations can then be expanded, dropping the *F* and *S* subscripts for notational simplicity:

30

31

where * δγ* is a function of

**and**

*δ*g_{A}**, but**

*δ*S**is not.**

*γ*The covariance is then:

32

Recognizing that ** E[δη]** is zero and

**S**

^{+}

*γ*and

**g**are equivalent allows us to simplify the covariance expression to be:

33

The full expression for is:

34

and, dropping the second-order error term, results in further simplification:

35

Therefore:

36

The covariance between star measurements **P*** _{s}* is:

37

Each element of this matrix can be expressed in terms of contributions from individual stars:

38

The terms are scalars, so they can be reordered giving:

39

where:

40

It is reasonable to assume that star vector errors are uncorrelated, so:

41

and, hence, **P*** _{s}* is a diagonal matrix.

If we define the gravity vector error covariance, **R*** _{g}*, as:

42

and recognize from our definitions:

43

we can then substitute these results into Equation (33) and simplify the expression to be:

44

To eliminate the dependence on global position **g**, we recall that:

45

from which it follows that:

46

An important note is that the elements of *R _{s,i}* and

*R*represent the covariance of their respective vectors’ elements. They do not represent the rotational noise along each sensor’s principal axes that is typically provided by sensor manufacturers. These two quantities can be related using the QUEST measurement model (Crassidis et al., 2007; Shuster & Oh, 1981):

_{g}47

48

where σ* _{s}* is the star tracker centroid noise,

*f*is the star tracker focal length,

*∆*is the star tracker pixel pitch, and σ

*is the per-axis measurement noise of the inclinometer.*

_{g}We now wish to relate the covariance **P** of that we derived in the previous section to a localization error along the planet’s surface, σ* _{geo}*. In order to do this, we need to project

**P**onto the plane tangent to the planet’s surface at the estimated global position. This derivation is made easier since is normal for such a plane.

We accomplish this by constructing a projection matrix defined by two arbitrary linearly independent vectors orthogonal to , **a** and **b**. The projection matrix, *𝒫*, is then:

49

where:

50

The covariance **P** projected onto the plane normal to is:

51

If we look at the eigenvalues of **P**_{proj}, we find that one of them is zero. Total localization error, σ* _{geo}*, can now be estimated by taking the root sum of the non-zero eigenvalues

*l*

_{1}and

*l*

_{2}:

52

The eigenvalues can also be considered independently as the principal variances of the localization noise.

## 3 CALIBRATION METHODOLOGY

In developing our localization solution, we had to also be aware of the need for system calibration. Many of the component transforms used in the position solution must be measured through a calibration process. This section examines the calibration parameterization along with a detailed discussion on compatible testing procedures and cost-function metrics.

### 3.1 Cost Functions

The calibration cost function quantifies the difference between an estimated position and the known position. DSS calibration is performed at a single known location and at several orientations. The resulting residual calculation can be formulated in terms of several different quantities that include star vector elevation, global position, star tracker quaternions, and inclinometer measurements. Having tested these error formulations in practice, cost function optimizations based on star-vector elevation tend to exhibit fewer local minima and lower localization error than other approaches. Thus, during calibration we use the known position and Equation (16), varying the calibration parameters in order to minimize ||**e**^{2}||.

### 3.2 Calibration Parameters

This section introduces the calibration parameters of the DSS. These include the mounting rotation error, measurement timing error introduced by clock offsets and latency, and latitude bias introduced by errors in the planet’s gravitational model.

#### 3.2.1 Mounting Error

The mounting rotation errors between the star tracker and inclinometer are captured in the **C*** _{SA}* rotation used in the measurement model in Equation (21). In practical applications, this rotation represents both the nominal rotation and misalignment. The mounting rotation parameters correct misalignment between the inclinometer frame, , and individual inclinometer subsensors frames, and .

The resulting mounting angle model augments Equation (7) to be:

53

Both **C*** _{A,x}* and

**C**

*are rotation matrices with three degrees of freedom defined by the small angle rotation vectors and , respectively. The implication of this formulation is that any measurement bias in the*

_{A,y}*x*and

*y*subsensors is equivalent to the mounting error along the respective subsensor’s measurement axis.

#### 3.2.2 Timing

Timing information is used to determine planetary rotation such that any timing error could manifest as bulk longitude error. This error could be introduced by improperly synchronized clocks or, more likely, data latency between sensors and the central data acquisition computer. Such latency could result from processing time, transmit times, or poor clock resolution.

To calculate **C*** _{FI}* on Earth, we derive International Atomic Time from GPS time and use it to compute the IAU 2000/2006 reduction (Petit & Luzum, 2010). This reduction provides

**C**

*by accounting for variations in the Earth’s rotation and precession about its axis.*

_{FI}To simplify calculations, following an initial IAU 2000/2006 reduction, rotation was treated as a pure -rotation between and , neglecting the effects of polar motion, nutation, and precession over the acquisition interval.

54

where is the timing error and *t _{d}* is the duration of a sidereal day.

This assumption was tested by comparing the average difference in when calculated using the method in Equation (54) versus using the IAU 2000/2006 reduction. Our results are from 100 data points collected over the course of an hour. The result had an average difference of 0.0005° or an equivalent error at Earth’s equator of 55 m. Additionally, there is no measurable difference in the calibration cost function residual from this approximation.

Although it is possible to include the timing error parameter directly in the IAU 2000/2006 reduction, such a solution adds significant computation to the calibration optimization. For our DSS system, the simplification in Equation (54) provides a substantial speedup with little impact on the final residuals. A two-stage optimization would likely be a good compromise when calibrating a real DSS system.

#### 3.2.3 Latitude Bias

The DSS is dependent on Earth’s geodesy and how it is modeled. While GPS uses the WGS 84 model, the new EGM2008 standard offers a typical 1.1 to 1.3 arc second improvement in local gravity corrections over the entire globe relative to the WGS 84 standard Pavlis et al. (2012). To account for this, a latitude bias, , was added to allow for approximately one arc second of constant error:

55

Correcting longitudinal differences between EGM2008 and WGS 84 is handled by the timing bias.

#### 3.2.4 Temperature Correction

The manufacturer of the electrolytic inclinometer used in this paper (Jewell Instruments) suggests a correction model for each of the two subsensors of the biaxial inclinometer (Jewell Instruments, 2018) of:

56

where , , and are calibration parameters, is the measured temperature, and *T _{cal}* is a reference temperature at which the unit was calibrated. For a given test,

*T*is the temperature at which the inclinometer is calibrated. The details of inclinometer calibration are discussed in Jovanovic and Enright (2020).

_{cal}#### 3.2.5 Parameter Summary

Table 1 lists the calibration parameters as well as the equation number in which they are applied for convenience. All of these parameters have nominal values of zero except *S*, which has a nominal value of one.

## 4 CALIBRATION

DSS calibration was done using field data collected from a Sinclair Interplanetary ST-16RT star tracker and a Jewell Instruments Tuff Tilt digital D801 wide-angle inclinometer. The calibration procedure discussed in this section involves rotating the DSS through a number of orientations. Measurements from the inclinometer and star tracker are recorded at each orientation while static. The goal of calibration is to find a set of parameters that make location estimates of the DSS invariant to changes in orientation. Two mechanizations were considered, a strapdown mechanization in which a DSS is on the body of a rover allowed to incline by approximately 10°. The other was an articulate mechanization in which the DSS is on an articulate boom on the rover and held level. Each of these were evaluated against a common set of error metrics.

### 4.1 Test Hardware

For DSS field testing, we used a Jewell Instrument (formerly Applied Geomatics) Tuff Tilt digital electrolytic biaxial inclinometer and a Sinclair Interplanetary ST-16 star tracker. Key performance metrics for each are tabulated in Table 2 and Table 4.

### 4.2 Optimizer

For this paper, we chose to use the Matlab lsqnonlin. We compared this solver to the Matlab fmincon and fminsearch solvers, observing significant differences in the final solution. The lsqnonlin solver arrived at a solution faster than the other two, offering similar positioning estimates. All the solvers performed consistently and did not find better solutions which would be indicative of local minima influencing the solution.

### 4.3 Error Metrics

Inadequate calibration of a DSS at multiple orientations tends to produce position estimates that are clustered by orientation. To track this aspect of performance as well as overall error with respect to a known truth, we developed multiple error metrics. Each of these metrics were derived from error along the planet’s surface between any points (*λ, φ*) and (*λ**, *ϕ**), expressed as:

57

Given an orientation indexed by *k*, the error metrics considered are illustrated in Figure 4. For individual system measurements with a common orientation indexed by *m*, the error within each orientation is given by:

58

where and are the mean positions at orientation *k*.

The error between orientations is:

59

where *λ* and *ϕ* are true positions provided by GPS.

Finally, the total error is represented as:

60

### 4.4 Strapdown Mechanization

The strapdown calibration data collection involved changing the DSS heading and inclination in multiple combinations. This process, illustrated by Figure 5, used a motorized telescope mount to acquire data from multiple combinations of tip, tilt, and heading. The exact procedure was:

The DSS was placed approximately level and a set of system measurements were recorded.

The system was titled by approximately 10° and Step 1 was repeated.

The system was leveled and Step 2 was repeated for three other rotation axes approximately 90° apart.

The system was rotated about the zenith by approximately 90° and Steps 1, 2, and 3 were repeated.

Once four distinct headings were completed, data acquisition was complete.

Following the procedure described above, σ* _{w}* was 527 m while σ

*was 1.78 km (see Table 5). This level of error is expected based on the inclinometer validation described in Jovanovic and Enright (2020) in which pockets of the inclinometer’s validation space had residuals in excess of 0.01° (equivalently 1 km on Earth’s surface).*

_{b}### 4.5 Articulate Mechanization

Another popular approach to calibrating a DSS is to constrain the range of motion of the inclinometer (Wei et al., 2019; Ning & Fang, 2009; Zhan et al., 2020). We refer to this as articulate mechanization. This method mounts the DSS to a motorized platform to ensure that the inclinometer remains level despite any orientation changes experienced by the rover body. In this study, the motorized telescope mount used to tilt the platform is, instead, used to level the DSS using feedback from the inclinometer. This keeps the DSS z-axis within 0.02° of the local vertical. Over this range, the inclinometer has a standard deviation on repeated measurements of 0.001°.

This mechanization effectively removes systematic inclinometer errors from the DSS measurements, significantly improving performance. Extensive testing with several types of inclinometers reveals that accuracy nonuniformities of 0.01° are common across wide-angle inclinometers and are extremely challenging to remove without time-consuming and dense characterization.

The articulate mechanization relies on a modified calibration procedure during which the DSS is held level at six unique headings. This reduces the degrees of freedom captured in the calibration data and renders some calibration parameters redundant. Table 5 shows the performance for an articulate calibration data set using different parameter combinations. From these results, the detrimental effect of utilizing the inclination dependent parameters in an articulate mechanization are apparent. Their omission reduced systemic errors, σ* _{b}*, by a factor of two. Conversely, omitting any parameters in the strapdown mechanization increases σ

*.*

_{b}### 4.6 Parameter Observability

Calibration parameters generally reflect physical models of device operations. Although these models are distinct, their effects on the calibration cost function are sometimes similar. We can look at the parameter observability and separability to assess the extent to which distinct parameters have similar effects in the calibration models. The latter quality contributes to a well-behaved calibration with well-defined parameter effects. To do this, we analyzed the singular value decomposition (SVD) of the calibration Jacobian:

61

where **e** is the cost function and **x** is the parameter set. The SVD of **H** can be writtten:

62

where **U** and **V** are unitary matrices and Σ is a diagonal matrix. By these definitions, the values in Σ (known as singular values) represent the sensitivity of the cost function to combinations of parameters while **V** encodes those combinations.

Figure 6 is a visual representation of the elements of **V*** ^{T}* with the columns corresponding to parameters and the rows corresponding to singular values.

There is an unexpected relationship between and . Such a relationship would be expected in the presence of a correlation between temperature and, in this case, the *y*-axis inclinometer reading. Testing the correlation between these two factors shows a correlation coefficient of 0.63. This correlation is rather significant and is the likely explanation for the observed relationship between and . There is no clear cause for such a correlation other than coincidence. The test data was collected over a 35-minute period, during which the controlled temperature experienced two sudden drops. The first drop was 0.45°C and made up the majority of the 0.81°C variability in temperature during the test. This drop coincided with the inclinometer being moved from a *y*-axis inclination of 10° to one of −10°. The second temperature drop was 0.18°C and coincided with the inclinometer *y*-axis measurement going from 10° to 0°. The temperature drops were anomalous and likely caused by sudden changes to convective cooling such as sustained gusts of wind.

Additionally, the *r*_{x,2} and *r*_{y,1} parameters were coupled but separable. The cost function was expected to be the least sensitive to these parameters and this is borne out of the experiments. The omission of these parameters degrades global localization accuracy. For these reasons, we included these parameters with the understanding that potential qualitative effects of these parameters coupling would be minimal due to the cost function’s insensitivity to them.

With the exception of these two coupling parameters, we see from Figure 6 that the parameters were observable and well separated.

## 5 MONTE CARLO VALIDATION COVARIANCE MODEL

The covariance model was validated using a Monte Carlo simulation which varied the noise statistics and orientation of both the inclinometer and star tracker. Using the star tracker’s onboard star catalog, synthetic inclinometer and star tracker measurements were generated for global position, orientation, and time.

The parameters of the simulated star tracker are tabulated in Table 6. Complete star tracker images were not synthesized; rather, catalog stars were passed through a pinhole camera model and had centroid noise applied. Individual star centroid noise was not adjusted for each star’s visual magnitude. A conservative value indicative of dim stars was used.

In the first Monte Carlo simulation, the values of σ* _{g}* and σ

*varied together such that the angular noise in the measured star and gravity vectors were equivalent. This represents a scenario in which neither sensor acts as a substantial performance bottleneck. Figure 7 shows the simulation predictions with input and output noise values expressed in equivalent distance along the Earth’s surface.*

_{s}The trial, shown in Figure 7, sampled equivalent^{3} measurement noise in the range of 6 m to 1.8 km in 300 equal intervals. Each interval was evaluated using a random heading (uniform over 360°) and the star tracker boresight was randomly perturbed from zenith (normal distribution with standard deviation σ = 2°). At each interval, 500 star tracker images and inclinometer measurements were generated.

Centroid and inclinometer noise were added to each measurement according to the assumed instrument covariance values. A maximum of 10 star vectors were generated for each exposure, selected on the basis of star brightness. This process mimicked the star tracker’s selection logic^{4}. The two non-zero eigenvalues of **P**_{proj} were used to predict measurement covariance and the close agreement between synthetically generated data and analytical predictions demonstrates the accuracy of our covariance model.

The sensor represented in these simulations was oriented to point near the zenith, giving a near-isotropic uncorrelated error distribution in the easting and northing directions. When pointed at extreme angles near the horizon, there was more asymmetry and correlation in the error distributions.

The random variations in orientation allowed stars to enter and exit the field of view between intervals such that favorable star geometries were transient. This introduced variability in both the synthetically generated eigenvalues as well as the predicted ones. Figure 8 shows the same results as Figure 7, but with fixed orientations. We still wished to ensure the solution was heading-independent, so instead of fixing the heading to a single value, four pre-determined headings were simulated at each interval. This is in contrast to the previously discussed process using random headings and a 2° random variation in inclination. In doing so, the same stars would always be visible and the variability in predicted positioning error would vanish.

We note that the bias visible in Figure 8 is an artifact of the smallest observable eigenvalue, described in Section 2.5, having a small non-zero value. Its removal via projection imparted a small in-plane rotation between the remaining predicted and generated eigenvectors shown in Figure 9. This rotation was random and varied from interval to interval and so, too, did the visible bias in both magnitude and sign. Note, the ellipse in Figure 9 is a hand-picked sample which is representative of typical performance but these plots vary from showing near-perfect agreement between predicted and generated data, to as much as 15% disagreement in principal axes magnitudes. A bias resulting from this rotation implies that if one direction under-predicts, the other will tend to over-predict, which we see in other trials as well as in Figure 8. This random rotation also implies that the eigenvectors do not exactly represent east-west and north-south, but rather tend to be predominantly in those directions.

In practice, the best performance achievable with our formulation had an accuracy of 296 m between orientations using articulate mechanization. Using the star tracker manufacturer’s quoted performance estimate of 7.2 arc seconds and the inclinometer’s *x*- and *y*-axis repeatability of 0.001° (3.6 arc seconds), which we’ve confirmed across multiple data sets, a 500-sample simulation was conducted and the covariance was predicted to have an error of 241 m, a notable under-prediction. We have confirmed that this error was not the result of our simplification of .

We are aware that our inclinometer calibration resolution was inadequate for a range of motion of ± 0.02° (Jovanovic & Enright, 2020), creating the potential for systemic errors in addition to measurement noise. Fine-scale surveys of our inclinometer calibration’s residual space around measurement values of 0° revealed that unmitigated curvatures exist that could introduce up to 0.0005° of error over a ± 0.02° range of motion. This happens to be the amount of noise that was needed to bring our prediction up to 294 m, however, the manner in which we discovered this value was not rigorous and we will therefore refrain from drawing conclusions based on it.

For the strapdown articulation, our model predicted σ* _{b}* of 1.58 km versus the experimental 1.78 km. In this case, we used the inclinometer wide-range root-mean-square error (RMSE) validation accuracy of 0.0038° (Jovanovic & Enright, 2020). The over-prediction of performance was likely due to a fine-scale systemic errors in our inclinometer which exceeded our validation resolution.

Figure 10 shows performance predictions obtained by independently varying the noise of both sensors and averaging them over multiple trials. Since inclination away from the zenith has a predicted adverse effect on performance, Figure 10 was generated for the two cases captured by our field tests. The combined error, σ* _{geo}*, was the root-sum of the two non-zero eigenvalues of

**P**

*, converted to an equivalent linear distance error. Each data point represented the RMSE across 500 trials with randomized orientations. The quantity modeled by this figure need not necessarily be temporal noise. The listed inclinometer error was the per-axis error while the star tracker error was the individual centroid error.*

_{proj}The red dots in Figure 10 represent the predicted DSS performance for our system for each mechanization. The difference between predicted performance of the red dots is the realized improvement of using a smaller region of our inclinometer’s measurement space, exposing less of the inclinometer’s systemic errors. We also plotted on this figure a black line to indicate the threshold at which improving a single sensor would come with diminishing benefits to overall DSS performance. This threshold corresponds to:

63

for the strapdown (± 10°) case and:

64

for the articulate (± 0.02°) case.

Further improvement to our DSS system from the articulate mechanization, since it is close to this threshold, is only possible by improving both the star tracker and inclinometer.

## 6 CONCLUSION

This study presents a suite of mathematical tools based on a forward measurement model for conducting DSS calibration without the use of heading truth. Using SVD, we demonstrate that the accompanying cost function produces parameter estimates that are observable and uncoupled in a strapdown implementation. Presented data suggests that the set of parameters must be reduced for an articulate mechanization in order to maximize localization accuracy. For calibration, an articulate mechanization provides superior performance with a simpler calibration model. Our results did not achieve the same accuracy as other research in this field, but we believe this to be a limitation of our choice of hardware and not a deficiency of our approach.

As presented, this localization method relies on GPS truth for calibration. Near-term interplanetary missions cannot count on any such positioning system, so some modifications of the outlined procedures would be necessary to make this approach practical. Two observations can be used to guide designers. First, we noted that other sources of localization could be substituted for the GPS positioning without any procedural modifications. Thus, even a single-point observation—using radio ranging, orbital imagery, etc.—would be sufficient to replicate this calibration in situ. Second, we recognized that the primary purpose of the calibration was to compensate for inter-sensor geometry and sensor biases. We expect that careful design and validation (e.g., athermal mechanical design, temperature regulation, ground-based calibrations) would greatly reduce the need for post-launch recalibration. Some technical risks may remain but these would be present in any system dependent on autonomous localization.

The accompanying covariance model was shown to be effective by a Monte Carlo analysis. When used to predict field test data, the deficiency of our inclinometer calibration likely detracted from our ability to make accurate predictions. Based on this Monte Carlo study, we provided a framework for predicting DSS performance. Our results also identified the threshold at which further DSS performance improvement would be equally limited by star tracker centroid noise and inclinometer subsensor noise.

Using Figure 10, we can conclude that in order to optimally achieve a strapdown DSS localization accuracy of 100 m, a star tracker with a star vector accuracy of 2.1 arc seconds and an inclinometer with a subsensor accuracy of 1.4 arc seconds would be required. This level of performance lies within reach of contemporary technology. With careful calibration of an inclinometer, such a system would be practical. Improving the performance further to 25 m would require additional improvements to the star tracker and inclinometer. The necessary 0.55 arc second star tracker accuracy is within reach of current technology, but achieving 0.4 arc second performance will be a challenge for inclinometers without substantially limiting an inclinometer’s functional range of motion. Further innovations in calibration and measurement technologies would be required to advance DSS systems to this level of performance.

## HOW TO CITE THIS ARTICLE

Jovanovic, I., & Enright, J. (2022) Deterministic heading-independent celestial localization measurement model. *NAVIGATION, 69*(3). https://doi.org/10.33012/navi.529

## CONFLICT OF INTEREST

The authors declare no potential conflict of interests.

## ACKNOWLEDGMENTS

The authors would like acknowledge Sinclair Interplanetary for generously supplying an engineering model star tracker for this research.

## Footnotes

**Present address**Department of Aerospace Engineering, 350 Victoria St., Toronto, ON, CANADA, M5B 2K3↵1 We use the relations developed by Shuster (2003) for stellar aberration correction

↵2 In the two star case, there is a 2 : 1 solution ambiguity.

↵3 Sensor angular error converted to angular error along the Earth’s geoid

↵4 Adding additional stars can give some accuracy improvements but this quickly reaches the point of diminishing returns. Moreover, these star-rich areas are found in only some parts of the sky.

This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.