Short-Arc Horizon-Based Optical Navigation by Total Least-Squares Estimation

Abstract

Horizon-based optical navigation (OPNAV) is an attractive solution for deep space exploration missions, with strong autonomy and high accuracy. In some scenarios, especially those with large variations in spacecraft distance from celestial bodies, the visible horizon arc could be very short. In this case, the traditional Christian–Robinson algorithm with least-squares (LS) estimation is inappropriate and would introduce a large mean residual that can be even larger than the standard deviation (STD). To solve this problem, a simplified measurement covariance model was proposed by analyzing the propagation of measurement errors. Then, an unbiased solution with the element-wise total least-squares (EW-TLS) algorithm was developed in which the measurement equation and the covariance of each measurement are fully considered. To further simplify this problem, an approximate generalized total least-squares algorithm (AG-TLS) was then proposed, which achieves a non-iterative solution by using approximate measurement covariances. The covariance analysis and numerical simulations show that the proposed algorithms have impressive advantages in the short-arc horizon scenario, for the mean residuals are always close to zero. Compared with the EW-TLS algorithm, the AG-TLS algorithm trades a negligible accuracy loss for a huge reduction in execution time and achieves a computing speed comparable to the traditional algorithm. Furthermore, a simulated navigation scenario reveals that a short-arc horizon can provide reliable position estimates for planetary exploration missions.

1. Introduction

Autonomous navigation is being widely valued with the bloom of deep space exploration missions [1,2]. It is an effective supplement and backup for ground-based tracking to reduce ground maintenance costs, overcome communication delays, and prevent communication loss [3]. Optical navigation is an image-based solution with strong autonomy, high accuracy, and good real-time performance [4]. It has been successfully applied to interplanetary exploration missions such as Apollo [5], Deep Space 1 [6], and SMART-1 [7] and will soon be applied to Orion missions [8].

OPNAV comes in many forms, including extracting the line-of-sight direction of natural or artificial targets [9,10,11], using the horizon of a celestial body [12,13,14,15], and tracking the landmarks on the body’s surface [16,17,18,19]. The most appealing solution for a planetary exploration mission is the horizon-based OPNAV, which guarantees a good balance of accuracy and complexity. The position of the spacecraft can be determined by using the planet’s centroid and apparent diameter in the image [12]. In most applications, the OPNAV camera’s field of view (FOV) is carefully designed to exactly cover the entire celestial body for better navigation accuracy. However, in some cases, the visible horizon arc could be very short, for example, due to the influence of lighting conditions, deviation of the pointing angle, and, most notably, change in altitude. In a typical Mars exploration mission, the distance between the spacecraft and Mars varies from a few hundred kilometers to tens of thousands of kilometers [1]. To ensure full-stage navigation accuracy, the FOV of the navigation camera is generally small such that only a small part of the horizon or scenery can be captured at the perigee. Although landmark-based OPNAV techniques can provide higher accuracy, they imply additional modes of navigation sensors and higher computational requirements. Thus, OPNAV, which only relies on a short-arc horizon, is essential.

The perspective projection of an ellipsoid planet or moon is an approximate ellipse, so many ellipse fitting algorithms have been used for the problem of horizon-based OPNAV [12,20]. In some applications, the celestial body with low ellipticity is directly approximated as a sphere, resulting in a simple and robust solution at the expense of some accuracy [21,22]. The traditional method of direct ellipse fitting in the image plane will incorrectly converge and introduce large errors in the case of a short-arc horizon. To solve this problem, Hu et al. [23] takes ephemeris and the spacecraft’s position as prior information and simplifies the ellipsoid to a sphere, thus obtaining the size and shape information of the ellipse in the image plane beforehand. The fitting accuracy of a short-arc horizon is greatly improved because there are only two unknown pieces of information, a scalar of the ellipse’s rotation and a two-dimensional vector of the ellipse’s translation on the image plane, to be estimated. This approach is inappropriate for our application due to the lack of position information, but it still reveals the difficulty of short-arc ellipse fitting and the importance of using prior information. The state-of-the-art solution is the Christian–Robinson algorithm [13,24,25], in which the ellipsoid celestial body is transferred into a sphere through Cholesky factorization, thus leading to a non-iterative solution. This is undoubtedly a novel and elegant solution that makes full use of prior information and greatly simplifies the problem. It has been employed in evaluating various OPNAV missions and has shown excellent performance [21,26].

In the Christian–Robinson algorithm, it was mentioned that the measurement equations should be solved using either the least-square (LS) estimation or the total least-square (TLS) estimation, but no detailed analysis was given [25], and LS estimation has widely been used for its simplicity [21,24,26]. The ordinary LS estimation assumes that the basis function $\mathbf{H}$ is error-free. However, in the measurement equation of the Christian–Robinson algorithm, the measurement errors are entirely embedded in the basis function $\mathbf{H}$. TLS estimation addresses errors in the $\mathbf{H}$ matrix, thus obtaining higher estimation accuracy and smaller estimation bias than the LS estimation [27]. TLS estimation is an asymptotic unbiased estimator whose error covariance achieves the Cramér–Rao lower bound [28,29]. It has been widely used in various measurement and estimation applications and demonstrates excellent performance [30,31,32,33,34]. However, compared with LS estimation, TLS estimation is more complex to apply and brings a non-negligible time consumption. In addition, TLS estimation is a general term for a class of algorithms with many variants [35,36,37]. It is necessary to further analyze why and when TLS estimation is needed and how it is implemented.

In this paper, we concentrate on improving the performance of horizon-based OPNAV based on a short-arc horizon. The main contributions of this paper lie in the following aspects: (1) Theoretical analysis reveals that LS estimation applied to the traditional Christian–Robinson algorithm is inappropriate and may introduce estimation bias. (2) A simplified measurement covariance model is proposed by further analyzing the propagation of measurement errors. The covariance model is no longer limited by the narrow FOV approximation and can be adapted to large FOV cameras. (3) An unbiased solution with the EW-TLS algorithm is developed. By considering the measurement equation and the covariance of each measurement, mean residuals close to zero can be achieved. (4) To further simplify this problem, an AG-TLS algorithm is proposed, in which a non-iterative solution is achieved by using the covariance of any one measurement as an approximation of all measurement covariances. Compared with the EW-TLS algorithm, the AG-TLS algorithm trades a negligible accuracy loss for a huge reduction in execution time. Its computing speed is comparable to the traditional algorithm and is suitable for onboard applications. (5) Covariance analysis and numerical results prove that the covariances of the LS estimation and the TLS estimation in this problem are identical, and the difference between them lies in the possible estimation bias. Due to the drawbacks of LS estimation, the Christian–Robinson method with ordinary LS estimation introduces large mean residuals under a short-arc horizon scenario, while the EW-TLS algorithm can achieve near zero mean residuals.

The remainder of this paper is organized as follows. Section 2 reviews the horizon-based OPNAV problem. Section 3 analyzes the problem in detail at the theoretical level and describes the proposed algorithms using TLS estimation. Section 4 presents the results of numerical simulations, and the main conclusions are summarized in Section 5.

2. Brief Review of Horizon-Based Optical Navigation

2.1. Geometry Fundamentals

The geometry of horizon-based OPNAV is shown in Figure 1. When the navigation camera takes a picture of an ellipsoid planet, the vectors from the camera’s optical center to the planet’s apparent horizon will form an approximate conical surface that tightly bounds the planet. Here is a brief review of the geometry fundamentals behind this.

To clarify the statement, several reference frames are introduced, which are defined as follows. The planetary frame (P) is established with the center of the planet ($O_P$) as the origin. The X-axis($X_P$), Y-axis($Y_P$), and Z-axis($Z_P$) of frame P are directions of the planet’s three principal axes, respectively. The origin of the camera frame (C) is the optical center of the camera ($O_c$). The X-axis ($X_c$) and Y-axis ($Y_c$) of frame C are directions of the row and column of the image sensor, respectively. The Z-axis ($Z_c$) of frame C is normal to the $X_c-Y_c$ plane. The image frame takes the imaging center of the image sensor ($O_s$) as the origin. The row and column of the image sensor are taken as the X-axis (x) and Y-axis (y) of the image frame, and the unit is a pixel. The image plane is right in the image frame, while its size is limited by the image sensor.

Note that the following analysis is performed in the camera frame C unless otherwise indicated. For an ideal navigation camera that perfectly follows the perspective imaging relationship [38], a detected planetary horizon point ${\mathbf{\psi }}_i={\left[x_i,y_i\right]}^T$ in the image frame corresponding to the following horizon vector ${\mathbf{s}}_i$:

$${\mathbf{s}}_i={\left[\begin{array}{ccc}{x}_i& y_i& f\end{array}\right]}^T$$

(1)

where f is the focal length of the navigation camera in pixels. The shape of most regular planets can be approximately described by a triaxial ellipsoid, then all the apparent horizons form a conical surface tightly bounding the triaxial ellipsoid during the imaging process [12]. For any horizon vector ${\mathbf{s}}_i$ belonging to the conical surface, it obeys the following constraints:

$${{\mathbf{s}}_i}^T\mathbf{M}{\mathbf{s}}_{\mathbf{i}}=0$$

(2)

where $\mathbf{M}$ is a full rank symmetric $3\times 3$ matrix. The expression of $\mathbf{M}$ is as follows:

$$\mathbf{M}=\mathbf{Ar}{\mathbf{r}}^T\mathbf{A}-\left({\mathbf{r}}^T\mathbf{Ar}-1\right)\mathbf{A}$$

(3)

where $\mathbf{r}$ is the position of camera center $O_c$ relative to planet center $O_P$, and $\mathbf{A}$ describes the shape and orientation of a triaxial ellipsoid. The expression of $\mathbf{A}$ is as follows:

$$\mathbf{A}={\mathbf{T}}_C^P\left[\begin{array}{ccc}1/a^2& 0& 0\\ 0& 1/b^2& 0\\ 0& 0& 1/c^2\end{array}\right]{\mathbf{T}}_P^C$$

(4)

where $\left\{a,b,c\right\}$ represents the length of the planet’s three principal axis, and ${\mathbf{T}}_C^P$ represents the rotation matrix from the planetary frame P to the camera frame C. With the development of high-precision ephemeris and star trackers, a high-precision ${\mathbf{T}}_C^P$ can be easily obtained, which will serve as prior information for the entire problem.

2.2. The Christian–Robinson Algorithm

The Christian–Robinson algorithm makes full use of prior ellipsoid shape and camera attitude information. The algorithm transforms the triaxial ellipsoid into a unit sphere through Cholesky factorization, thereby fully simplifying the triaxial ellipsoid projection constraint. Here is a brief review of the algorithm.

(1)

Convert the camera’s horizon measurement ${\mathbf{s}}_i$ to ${\overline{\mathbf{s}}}_i$ in Cholesky factorization space:

$${\overline{\mathbf{s}}}_i=\mathbf{U}{\mathbf{s}}_i$$

(5)

where $\mathbf{U}$ can be obtained from Equation (4) as:

$$\mathbf{U}=\mathbf{Q}{\mathbf{T}}_P^C$$

(6)

and the expression of $\mathbf{Q}$ is as follows:

$$\mathbf{Q}=\left[\begin{array}{ccc}1/a& 0& 0\\ 0& 1/b& 0\\ 0& 0& 1/c\end{array}\right]$$

(7)

(2)

Normalize ${\overline{\mathbf{s}}}_i$ to ${\mathbf{h}}_i$:

$${\mathbf{h}}_i=\frac{{\overline{\mathbf{s}}}_i}{∥{\overline{\mathbf{s}}}_i∥}$$

(8)

(3)

The problem is converted to solving for $\mathbf{n}$ in the following measurement equations:

$$\mathbf{Hn}=\mathbf{y}$$

(9)

where:

$$\mathbf{H}=\left[\begin{array}{c}{\mathbf{h}}_1^T\\ {\mathbf{h}}_2^T\\ ⋮\\ {\mathbf{h}}_n^T\end{array}\right]$$

(10)

$\mathbf{y}={\mathbf{1}}_{n\times 1}$, n is the total number of horizon measurement points.

(4)

The solution of $\mathbf{r}$ could be found by:

$$\mathbf{r}=-{\left({\mathbf{n}}^T\mathbf{n}-1\right)}^{-(1/2)}{\mathbf{U}}^{-1}\mathbf{n}$$

(11)

3. Total Least-Squares Estimation

3.1. Estimation Bias Caused by Least-Squares Estimation

With the help of the Christian–Robinson algorithm, the problem has been greatly simplified to solving for $\mathbf{n}$ in Equation (9). Although it was mentioned in [25] that the measurement equations should be solved using either the LS estimation or the TLS estimation, no detailed analysis was given. The LS estimation seems to bring the same results with much lower computational consumption in most cases, so it is more widely used [21,24,26]. However, the limitation of the LS estimation seriously affects the performance of the algorithm in short-arc horizon-based OPNAV. This will be demonstrated in the following derivations.

Let us revisit the problem described in Equation (9). If the problem is treated as a LS problem, then its loss function is:

$$J\left({\widehat{\mathbf{n}}}_{ls}\right)=\frac{1}{2}{\left(\tilde{\mathbf{y}}-\tilde{\mathbf{H}}{\widehat{\mathbf{n}}}_{ls}\right)}^T\left(\tilde{\mathbf{y}}-\tilde{\mathbf{H}}{\widehat{\mathbf{n}}}_{ls}\right)$$

(12)

where $\tilde{\mathbf{y}}$ is the measurement of $\mathbf{y}$, $\tilde{\mathbf{H}}$ is the measurement of $\mathbf{H}$, and ${\widehat{\mathbf{n}}}_{ls}$ is the LS solution. As shown in Section 2.2, $\tilde{\mathbf{H}}$ is converted from the camera’s horizon measurement by Equations (5) and (8), $\tilde{\mathbf{y}}$ is set as ${\mathbf{1}}_{n\times 1}$ directly from Equation (9), and the expression of ${\widehat{\mathbf{n}}}_{ls}$ is:

$${\widehat{\mathbf{n}}}_{ls}={\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T\tilde{\mathbf{y}}$$

(13)

If we treat ${\widehat{\mathbf{n}}}_{ls}$ as a function of the measurements, then the expectation of the LS solution is:

$$\begin{array}{cc}\hfill E\left({\widehat{\mathbf{n}}}_{ls}\right)& \hfill =E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T\tilde{\mathbf{y}}\right)\hfill \end{array}$$

(14)

According to Equation (9), we have:

$$\tilde{\mathbf{y}}={\mathbf{1}}_{n\times 1}=\mathbf{y}=\mathbf{Hn}$$

(15)

thus $\tilde{\mathbf{y}}$ can be rewritten as:

$$\tilde{\mathbf{y}}=\mathbf{Hn}=\left(\tilde{\mathbf{H}}+\mathbf{H}-\tilde{\mathbf{H}}\right)\mathbf{n}$$

(16)

Then, substituting Equation (16) into Equation (14) yields:

$$\begin{array}{cc}\hfill E\left({\widehat{\mathbf{n}}}_{ls}\right)& =E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T\tilde{\mathbf{y}}\right)\hfill \\ \hfill & =E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T(\tilde{\mathbf{H}}+\mathbf{H}-\tilde{\mathbf{H}})\mathbf{n}\right)\hfill \\ \hfill & =E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\mathbf{n}\right)+E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T(\mathbf{H}-\tilde{\mathbf{H}})\mathbf{n}\right)\hfill \\ \hfill & =\mathbf{n}+E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T(\mathbf{H}-\tilde{\mathbf{H}})\mathbf{n}\right)\hfill \end{array}$$

(17)

It can thereby be concluded that LS estimation is a biased estimation for this problem [27], and the estimation bias is:

$$E\left({\widehat{\mathbf{n}}}_{ls}\right)-\mathbf{n}=E\left({\left({\tilde{\mathbf{H}}}^T\tilde{\mathbf{H}}\right)}^{-1}{\tilde{\mathbf{H}}}^T(\mathbf{H}-\tilde{\mathbf{H}})\mathbf{n}\right)$$

(18)

Estimation bias is unavoidable in LS estimation due to the error contained in $\tilde{\mathbf{H}}$, and its existence is verified in subsequent simulations. The estimation bias will increase with the lack of measurement information due to the shortening of the horizon arc length. In practice, estimation biases are not easily detected without additional correction methods. Furthermore, it is undoubtedly catastrophic for subsequent Kalman filter or other estimation frameworks if the estimation bias is large enough compared with the estimation error.

3.2. Simplified Measurement Covariance Model

As shown in the last section, the key to the estimation bias is the coupling of the measurement error with the measurement model. To solve this problem, we need to find out the covariance of measurement errors beforehand, which means studying the propagation of the measurement error in the transition of horizon measurements in the image to $\tilde{\mathbf{H}}$ in the measurement model. The narrow FOV approximation was used in previous covariance derivations, which greatly limited its applicability [20,24,39,40]. In this paper, a more general and simple measurement covariance model is given by reducing the normalization step of the measurement vector.

At first, a detected horizon point ${\mathbf{\psi }}_i={\left[x_i,y_i\right]}^T$ in the image frame describes the horizon point’s line-of-sight direction. The measured error should be uniformly Gaussian distributed in the image plane, then the covariance of ${\mathbf{\psi }}_i$ can be expressed as:

$${\mathcal{R}}_{{\mathbf{\psi }}_i}=\left[\begin{array}{cc}{\sigma }^2& 0\\ 0& {\sigma }^2\end{array}\right]$$

(19)

where $\sigma$ is the STD of the horizon point location error, and the unit is one pixel.

Then, the covariance changes during the projection of the point ${\mathbf{\psi }}_i$ in the image frame to the vector ${\mathbf{s}}_i={\left[\begin{array}{ccc}{x}_i& y_i& f\end{array}\right]}^T$ in the camera frame C. The narrow FOV approximation is commonly used to estimate the covariances of ${\mathbf{s}}_i$ when the camera has a small FOV [20,24]. A small FOV indicates that the focal length f is much greater than the size of the image sensor, allowing for the assumption that the covariance of any normalized ${\mathbf{s}}_i$ within the narrow FOV remains constant [40]. However, we do not use the narrow FOV approximation in our covariance model because it would restrict the model from being applied to cameras with larger FOVs. Instead, we apply the normalization step only once at a later stage. This approach further simplifies the covariance calculation and increases the model’s versatility. The projection corresponds to the Jacobi matrix as follows:

$$\frac{\partial {\mathbf{s}}_i}{\partial {\mathbf{\psi }}_i}=\left[\begin{array}{cc}1& 0\\ 0& 1\\ 0& 0\end{array}\right]$$

(20)

Thus the covariance of ${\mathbf{s}}_i$ is:

$$\begin{array}{c}\hfill {\mathcal{R}}_{{\mathbf{s}}_i}=\frac{\partial {\mathbf{s}}_i}{\partial {\mathbf{\psi }}_i}{\mathcal{R}}_{{\mathbf{\psi }}_i}{\left(\frac{\partial {\mathbf{s}}_i}{\partial {\mathbf{\psi }}_i}\right)}^T=\left[\begin{array}{ccc}{\sigma }^2& 0& 0\\ 0& {\sigma }^2& 0\\ 0& 0& 0\end{array}\right]\end{array}$$

(21)

Due to the characteristics of camera measurements, the covariance of ${\mathbf{s}}_i$ is not uniformly distributed in three directions, and the errors mainly accumulate in $X_C$ and $Y_C$ directions. Note that a more general theory of camera measurement covariance could be found in [39], while the one used here is an adequate simplification.

Next, consider the transformation corresponding to Equation (5), in which ${\mathbf{s}}_i$ in the camera frame P is transferred into ${\overline{\mathbf{s}}}_i$ in the Cholesky factorization space. The transformation corresponds to the Jacobi matrix as follows:

$$\frac{\partial {\overline{\mathbf{s}}}_i}{\partial {\mathbf{s}}_i}=\mathbf{U}$$

(22)

In addition, as shown in Equation (8), the vector ${\overline{\mathbf{s}}}_i$ needs to be normalized to ${\mathbf{h}}_i$. The normalization corresponds to the Jacobi matrix as follows:

$$\frac{\partial {\mathbf{h}}_i}{\partial {\overline{\mathbf{s}}}_i}=\frac{1}{∥{\overline{\mathbf{s}}}_i∥}\left({\mathbf{I}}_{3\times 3}-{\mathbf{h}}_i{\mathbf{h}}_i^T\right)$$

(23)

Defining the $3\times 3$ matrix ${\mathsf{\Phi }}_i=\frac{\partial {\mathbf{h}}_i}{\partial {\overline{\mathbf{s}}}_i}$ thus obtains:

$$\frac{\partial {\mathbf{h}}_i}{\partial {\mathbf{s}}_i}={\mathsf{\Phi }}_i\mathbf{U}$$

(24)

Since the magnitude of every measurement is different, the covariance of each measurement varies after normalization.

Finally, the covariance of ${\mathbf{h}}_i$ is found to be:

$${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}={\mathsf{\Phi }}_i\mathbf{U}{\mathcal{R}}_{{\mathbf{s}}_i}{\mathbf{U}}^T{\mathsf{\Phi }}_i^T$$

(25)

3.3. Element-Wise Weighted Total Least-Squares Algorithm

As shown in the previous section, the measurement error propagation leads to different covariance for each measurement. In addition, the measurements $\mathbf{H}$ act as a model in the measurement equation, which means that this is an entirely error-in-variables problem. To better address this problem, the TLS estimation is adopted, which takes into account the error in the $\mathbf{H}$ matrix and weights the measurement equations based on the covariance of the measurements.

For this problem, define the following matrix:

$$\mathbf{D}=\left[\begin{array}{cc}\mathbf{H}& {\mathbf{1}}_{n\times 1}\end{array}\right]$$

(26)

$$\mathbf{z}={\left[\begin{array}{cc}{\mathbf{n}}^T& -1\end{array}\right]}^T$$

(27)

The constraint condition Equation (9) can then be rewritten as:

$$\mathbf{Dz}=\mathbf{0}$$

(28)

TLS solution ${\widehat{\mathbf{n}}}_{tls}$ for this problem is the optimal estimation of $\mathbf{n}$ that minimizes [28]:

$$\begin{array}{c}J\left({\widehat{\mathbf{n}}}_{tls}\right)=\frac{1}{2}{\mathrm{vec}}^T\left({\tilde{\mathbf{D}}}^T-{\widehat{\mathbf{D}}}^T\right){\mathcal{R}}^{-1}\mathrm{vec}\left({\tilde{\mathbf{D}}}^T-{\widehat{\mathbf{D}}}^T\right),\hfill \\ \mathrm{subject}\mathrm{to}\widehat{\mathbf{D}}\widehat{\mathbf{z}}=\mathbf{0}\hfill \end{array}$$

(29)

where $\widehat{\mathbf{z}}={\left[\begin{array}{cc}\widehat{\mathbf{n}}{}_{tls}^T& -1\end{array}\right]}^T$ denotes the estimation of $\mathbf{z}$, $\tilde{\mathbf{D}}=\left[\begin{array}{cc}\tilde{\mathbf{H}}& {\mathbf{1}}_{n\times 1}\end{array}\right]$ denotes the measurement of $\mathbf{D}$, $\widehat{\mathbf{D}}=\left[\begin{array}{cc}\widehat{\mathbf{H}}& {\mathbf{1}}_{n\times 1}\end{array}\right]$ denotes the estimation of $\mathbf{D}$, $\widehat{\mathbf{H}}$ denotes the estimation of $\mathbf{H}$, $\mathcal{R}$ is the covariance matrix of $\tilde{\mathbf{D}}$, and vec denotes a vector formed by stacking the consecutive columns of the associated matrix. $\mathcal{R}$ is given by following block diagonal matrix:

$$\mathcal{R}=\left[\begin{array}{cccc}{\mathcal{R}}_1& {\mathbf{0}}_{4\times 4}& \cdots & {\mathbf{0}}_{4\times 4}\\ {\mathbf{0}}_{4\times 4}& {\mathcal{R}}_2& \cdots & {\mathbf{0}}_{4\times 4}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathbf{0}}_{4\times 4}& {\mathbf{0}}_{4\times 4}& \cdots & {\mathcal{R}}_n\end{array}\right]$$

(30)

where ${\mathcal{R}}_i$ is covariance matrix of ${\mathbf{h}}_i$, and it is given by:

$${\mathcal{R}}_i=\left[\begin{array}{cc}\hfill {\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}\hfill & {\mathcal{R}}_{\mathbf{h}{\mathbf{y}}_i}\\ \hfill {\mathcal{R}}_{\mathbf{h}{\mathbf{y}}_i}^T\hfill & {\mathcal{R}}_{\mathbf{y}{\mathbf{y}}_i}\end{array}\right]$$

(31)

As defined in Equation (15), there is:

$${\mathcal{R}}_{\mathbf{h}{\mathbf{y}}_i}={\mathbf{0}}_{3\times 1}\mathrm{and}{\mathcal{R}}_{\mathbf{y}{\mathbf{y}}_i}={\mathbf{0}}_{1\times 1}$$

(32)

In this case, since the errors are element-wise correlated and non-stationary, each measurement element needs to be weighted according to its covariance. It implies that this is an element-wise weighted TLS (EW-TLS) problem and can be solved by using an iterative procedure [27,36]. In addition, it is necessary to impose strict time limits on the algorithm to ensure implementation onboard. The specific steps of the EW-TLS algorithm are as follows:

(1)

Calculate the covariance ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ of each set of measurements with Equation (25).

(2)

Take the LS estimation ${\widehat{\mathbf{n}}}_{ls}$ as the initial estimate ${\widehat{\mathbf{n}}}^{\left(1\right)}$ with Equation (13).

(3)

The estimate ${\widehat{\mathbf{n}}}^{(j+1)}$ at the ${(j+1)}^{th}$ iteration is given by [36]:

$${\widehat{\mathbf{n}}}^{(j+1)}={\left[\sum _{i=1}^n\frac{{\tilde{\mathbf{h}}}_i{\tilde{\mathbf{h}}}_i^T}{{\gamma }_i\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right)}-\frac{e_i^2\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right){\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}}{{\gamma }_i^2\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right)}\right]}^{-1}\left[\sum _{i=1}^n\frac{{\tilde{y}}_i{\tilde{\mathbf{h}}}_i}{{\gamma }_i\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right)}\right]$$

(33)

where:

$${\gamma }_i\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right)={\widehat{\mathbf{n}}}^{\left(j\right)T}{\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}{\widehat{\mathbf{n}}}^{\left(j\right)}$$

(34)

$$e_i\left({\widehat{\mathbf{n}}}^{\left(j\right)}\right)={\tilde{\mathbf{h}}}_i^T{\widehat{\mathbf{n}}}^{\left(j\right)}-{\tilde{\mathbf{y}}}_i$$

(35)

(4)

The iteration is stopped when the final change is less than the tolerance value, or when the number of iterations is greater than the maximum number of iterations (MaxIterations). The iteration stopping conditions are specified as:

$${∥{\widehat{\mathbf{n}}}^{(j+1)}-{\widehat{\mathbf{n}}}^{\left(j\right)}∥}_2\le \mathrm{tolerance}\mathrm{value}\mathrm{or}j>\mathrm{MaxIterations}$$

(36)

The tolerance value is set according to the calculation precision, and it is set to ${10}^{-10}$ in this paper. MaxIterations is limited to five to avoid excessive time consumption. Although it is certainly satisfactory that the algorithm usually converges in only one to three iteration cycles, it is definitely safer to impose a limit on the number of iterations.

3.4. Approximate Generalized Total Least-Squares Algorithm

Although the accuracy of the EW-TLS algorithm is satisfactory, its computation time is still too long. Considering the limitations of onboard computing capabilities, it is usually worthwhile to trade some precision loss for faster computation. To further simplify this problem, an AG-TLS algorithm is proposed in this paper, which uses the covariance of any one measurement as an approximation of all covariances, thus achieving a non-iterative solution.

Recall the analysis of measurement error propagation in Section 3.2. For a planet with a small oblateness, its three major axes a, b, c are very close, and Equation (6) can be approximated as:

$$\mathbf{U}=\mathbf{Q}{\mathbf{T}}_P^C\approx \frac{1}{R_e}{\mathbf{I}}_{3\times 3}{\mathbf{T}}_P^C=\frac{1}{R_e}{\mathbf{T}}_P^C$$

(37)

where $R_e$ is the mean radius of the target planet. Thus, ${\overline{\mathbf{s}}}_i$ may be rewritten as:

$${\overline{\mathbf{s}}}_i=\frac{1}{R_e}{\mathbf{T}}_P^C{\mathbf{s}}_{\mathbf{i}}$$

(38)

The attitude transformation matrix ${\mathbf{T}}_P^C$ does not change the norm of ${\mathbf{s}}_{\mathbf{i}}$, which lead to:

$$∥{\overline{\mathbf{s}}}_i∥=\frac{\sqrt{x_i^2+y_i^2+f^2}}{R_e}$$

(39)

Then, substituting Equations (38) and (39) into Equation (24) yields:

$$\frac{\partial {\mathbf{h}}_i}{\partial {\mathbf{s}}_i}\approx \frac{\sqrt{x_i^2+y_i^2+f^2}}{R_e^2}\left({\mathbf{I}}_{3\times 3}-{\mathbf{h}}_i{\mathbf{h}}_i^T\right){\mathbf{T}}_P^C$$

(40)

For a narrow FOV camera, the focal length f is much larger than $\sqrt{x_i^2+y_i^2}$ and is approximately:

$$\sqrt{x_i^2+y_i^2+f^2}\approx f$$

(41)

Further substituting Equations (40) and (41) into Equation (25) yields:

$${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}\approx \frac{f^2}{R_e^4}\left({\mathbf{I}}_{3\times 3}-{\mathbf{h}}_i{\mathbf{h}}_i^T\right){\mathbf{T}}_P^C{\mathcal{R}}_{{\mathbf{s}}_i}{\mathbf{T}}_C^P{\left({\mathbf{I}}_{3\times 3}-{\mathbf{h}}_i{\mathbf{h}}_i^T\right)}^T$$

(42)

As shown in Equation (21), ${\mathcal{R}}_{{\mathbf{s}}_i}$ is constant and mainly accumulated in $X_C$ and $Y_C$ directions. Then, as shown in Equation (42), when horizon measurements ${\mathbf{s}}_i$ in the image transitions to ${\mathbf{h}}_i$ in the measurement model, ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ are redistributed to three directions with $\left({\mathbf{I}}_{3\times 3}-{\mathbf{h}}_i{\mathbf{h}}_i^T\right){\mathbf{T}}_P^C$. The variation of ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ is dominated by the variation of ${\mathbf{h}}_i{\mathbf{h}}_i^T$, while ${\mathbf{h}}_i{\mathbf{h}}_i^T$ changes relatively little since all measurements ${\mathbf{h}}_i$ are gathered in a small area for a narrow FOV camera.

According to the above analysis, it can be concluded that when both camera’s FOV and planetary oblateness are small, ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ are relatively close for all measurements ${\mathbf{h}}_i$. Therefore, the problem can be further simplified with a small loss of precision by taking the covariance of any one measurement as an approximation of all covariances. Considering the error distribution and correlation among the elements, the simplified problem is a generalized TLS estimation problem, which has a closed-form solution [27,35,36]. The specific steps of the AG-TLS algorithm are as follows:

(1)

Calculate the covariance ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ of any one measurement ${\mathbf{h}}_i$ with Equation (25), then construct ${\mathcal{R}}_i$ with Equation (31).

(2)

Prepare for Cholesky factorization of matrix ${\mathcal{R}}_i$. According to Equation (32), ${\mathcal{R}}_{\mathbf{y}{\mathbf{y}}_i}$ should be equal to zero, while this will lead to a rank deficient of ${\mathcal{R}}_i$. In addition, the covariance matrix ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ is a positive semi-definite matrix. With the change in spacecraft position and attitude, and considering the machine computational accuracy, ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$ will sometimes have an eigenvalue close to zero or even negative. These factors block the Cholesky decomposition of the matrix ${\mathcal{R}}_i$ and affect the implementation of the algorithm. Thus, the approximate covariance matrix $\overline{\mathcal{R}}$ is obtained by making a simple but effective adjustment to ${\mathcal{R}}_i$.

$$\overline{\mathcal{R}}={\mathcal{R}}_i+\epsilon {\mathbf{I}}_{4\times 4}$$

(43)

where $\epsilon$ is a very small value, e.g., ${10}^{-15}$, and is much smaller than every element in the matrix ${\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}$. Such an approximation hardly affects the distribution of covariances and ensures that the matrix $\overline{\mathcal{R}}$ is positive definite.

(3)

Perform a Cholesky factorization on the matrix $\overline{\mathcal{R}}$:

$$\overline{\mathcal{R}}={\mathbf{C}}^T\mathbf{C}$$

(44)

Then partition the inverse of $\mathbf{C}$ as:

$${\mathbf{C}}^{-\mathbf{1}}=\left[\begin{array}{cc}\hfill {\mathbf{C}}_{11}\hfill & \hfill \mathbf{c}\hfill \\ \hfill {\mathbf{0}}_{1\times 3}\hfill & \hfill c_{22}\hfill \end{array}\right]$$

(45)

where ${\mathbf{C}}_{11}$ is an $3\times 3$ matrix, $\mathbf{c}$ is an $3\times 1$ vector, and $c_{22}$ is a scalar.

(4)

Perform a singular value decomposition on the following matrix:

$$\widehat{\mathbf{D}}{\mathbf{C}}^{-\mathbf{1}}=\mathbf{N}\mathbf{S}{\mathbf{V}}^T$$

(46)

Then partition $\mathbf{V}$ as:

$$\mathbf{V}=\left[\begin{array}{cc}\hfill {\mathbf{V}}_{11}\hfill & \hfill \mathbf{v}\hfill \\ \hfill {\mathbf{w}}^{\mathbf{T}}\hfill & \hfill v_{22}\hfill \end{array}\right]$$

(47)

where $\mathbf{N}$, $\mathbf{S}$ and $\mathbf{V}$ are $4\times 4$ matrixes, ${\mathbf{V}}_{11}$ is an $3\times 3$ matrix, $v_{22}$ is a scalar, and $\mathbf{w}$ and $\mathbf{v}$ are both $3\times 1$ vectors.

(5)

The final solution is given by [27]:

$${\widehat{\mathbf{n}}}_{agtls}=\frac{1}{c_{22}}\left(-\frac{1}{v_{22}}{\mathbf{C}}_{11}\mathbf{v}-\mathbf{c}\right)$$

(48)

3.5. Analytic Solution Covariance

In this section, the analytic solution covariance of the TLS estimation and the LS estimation are evaluated to compare their performance under this problem.

According to [28], the analytic solution covariance of $\mathbf{n}$ with TLS estimation is:

$${}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{n}}={\left(\sum _{i=1}^n\frac{{\mathbf{h}}_i{\mathbf{h}}_i^T}{{\mathbf{z}}^T{\mathcal{R}}_i\mathbf{z}}\right)}^{-1}$$

(49)

According to Equation (32), it can be simplified as:

$${}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{n}}={\left(\sum _{i=1}^n\frac{{\mathbf{h}}_i{\mathbf{h}}_i^T}{{\mathbf{n}}^T{\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}\mathbf{n}}\right)}^{-1}$$

(50)

Take the variation of Equation (11) with $\mathbf{n}$:

$$\frac{\partial \mathbf{r}}{\partial \mathbf{n}}=\mathbf{F}$$

(51)

where:

$$\mathbf{F}=-{\left({\mathbf{n}}^T\mathbf{n}-1\right)}^{-(1/2)}{\mathbf{T}}_C^P{\mathbf{Q}}^{-1}\left({\mathbf{I}}_{3\times 3}-\frac{\mathbf{n}{\mathbf{n}}^T}{{\mathbf{n}}^T\mathbf{n}-1}\right)$$

(52)

So the analytic solution covariance of $\mathbf{r}$ is:

$${}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{r}}=\mathbf{F}{}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{n}}{\mathbf{F}}^T$$

(53)

The analytic covariance of $\mathbf{n}$ with LS estimation is:

$${}^{\mathrm{ls}}{\mathbf{P}}_{\mathbf{n}}={\left(\sum _{i=1}^n\frac{{\mathbf{h}}_i{\mathbf{h}}_i^T}{{}^{\mathrm{ls}}{\mathcal{R}}_{\mathbf{y}{\mathbf{y}}_i}}\right)}^{-1}$$

(54)

Since ${\mathbf{h}}_i$ is assumed to be unperturbed in the ordinary LS estimation, there are:

$${}^{\mathrm{ls}}{\mathcal{R}}_{\mathbf{y}{\mathbf{y}}_i}={\mathbf{n}}^T{\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}\mathbf{n}$$

(55)

Therefore, it can be obtained:

$${}^{\mathrm{ls}}{\mathbf{P}}_{\mathbf{n}}={\left(\sum _{i=1}^n\frac{{\mathbf{h}}_i{\mathbf{h}}_i^T}{{\mathbf{n}}^T{\mathcal{R}}_{\mathbf{h}{\mathbf{h}}_i}\mathbf{n}}\right)}^{-1}={}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{n}}$$

(56)

That is to say:

$${}^{\mathrm{ls}}{\mathbf{P}}_{\mathbf{r}}=\mathbf{F}{}^{\mathrm{ls}}{\mathbf{P}}_{\mathbf{n}}{\mathbf{F}}^T={}^{\mathrm{tls}}{\mathbf{P}}_{\mathbf{r}}$$

(57)

It can be seen that the analytic covariances of the LS estimation and the TLS estimation are completely consistent, then the difference between the two kinds of estimations lies in the possible estimation bias. As demonstrated in [28], the TLS estimation is asymptotically unbiased, which is the main advantage over the ordinary LS estimation in this problem, so it is certainly attractive to adopt the TLS estimation without any loss of precision. In the actual performance of the EW-TLS algorithm, the analytic solution covariance will be gradually approximated as the iteration number grows, while in the AG-TLS algorithm, there may be a gap between the actual performance and the analytics solution covariance due to the presence of approximations. The actual performance of the proposed algorithms is verified in the following numerical simulations.

4. Numerical Results

4.1. Performance under a Short-Arc Horizon

The simulated navigation camera assumed a FOV of 8x8 degrees and a resolution of 1024 × 1024 pixels. The exact numerical solution of the apparent horizon in the image plane could be found using Equations (1)–(3). One horizon point is sampled from each pixel where the horizon appears, and Gaussian distributed errors with $\sigma =0.3$ pixel are added to the horizon points, thus getting closer to the actual horizon detection situation. The above simulation conditions of navigation camera are applied to all subsequent simulations. The simulated navigation camera is pointing to the center of Mars at a range of 65,000 km. At this point, the horizon of Mars occupies almost the entire FOV of the navigation camera. However, due to the possible deviation of the pointing angle and the influence of lighting conditions, only a small part of the horizon may be detected. Simulated images of the horizon ellipse with different arc lengths in the image frame are shown in Figure 2.

As an example scenario, assume that a 15-degree of horizon arc is detected, and the actual number of valid horizon points is about 100. For convenience, the Christian–Robinson algorithm with LS estimation is called the LS algorithm, while the algorithms proposed in this paper are called the EW-TLS algorithm and AG-TLS algorithm. A total of 5000 Monte Carlo simulations are performed with three algorithms. Various horizon measurements are produced by adding randomly generated Gaussian noise to the same scenario, and the results are shown in Figure 3. Errors of OPNAV are presented in the camera frame C. Gray points represent measurement residuals of Monte Carlo simulations, the black cross represents mean residuals of simulations, the red dashed line ellipse represents 3$\sigma$ bounds of analytic covariance, and the black solid line ellipse represents 3$\sigma$ bounds of numeric covariance from simulations. Note that the covariance in Figure 3 is superimposed on the mean residuals. In the simulations of all three algorithms, the analytical covariance and the numeric covariance are in good agreement, and the correctness of the analytic covariance is fully verified.

The detailed results of Monte Carlo simulations are listed in

Table 1. Monte Carlo results with a 15-degree arc horizon.

Algorithm	Direction	Mean Residuals, km	STD, km	MSTDR	RMSE, km
	$X_C$ axis	−296.82	95.25	311.63%	311.73
LS	$Y_C$ axis	−39.71	13.18	301.23%	41.84
	$Z_C$ axis	5717.96	1834.61	311.67%	6005.01
	$X_C$ axis	0.84	95.77	0.88%	95.76
EW−TLS	$Y_C$ axis	0.04	13.22	0.34%	13.21
	$Z_C$ axis	−16.24	1845.42	0.88%	1845.30
	$X_C$ axis	−1.89	95.87	1.97%	95.88
AG−TLS	$Y_C$ axis	−0.37	13.23	2.78%	13.24
	$Z_C$ axis	36.54	1847.53	1.97%	1847.71

. The magnitude of the solution covariances is the same for all three algorithms. However, as shown in Figure 3a, the mean residuals of the LS algorithm are clearly not zero. The mean to the standard deviation ratio (MSTDR) is introduced here to evaluate the influence of mean residuals. MSTDR is set to be a positive value. MSTDR of the LS algorithm in $X_C$, $Y_C$, and $Z_C$ directions are 311.63%, 301.23%, and 311.67%, respectively. This means that the mean residuals are much larger than the STD, which is unacceptable in practical applications. In contrast, as shown in Figure 3b,c, the MSTDR of the EW-TLS algorithm and AG-TLS algorithm are all near zero. The root mean squared error (RMSE) is introduced to evaluate the estimation effectiveness of all three algorithms with respect to mean residuals, and two TLS algorithms both show a significant accuracy advantage of three times over the LS algorithm. The mean residuals of the AG-TLS algorithm are slightly larger than the average error of the EW-TLS algorithm. From the perspective of MSTDR, the performance difference between the two algorithms is within 2%, while it is basically no different from the perspective of RMSE, which indicates that the approximation and simplification of the AG-TLS algorithm are effective under a short-arc horizon.

Although there is no circle fitting or ellipse fitting of the projected horizon in the image frame in the Christian–Robinson algorithm, it is still more intuitive to use the projected horizon in the image frame to show the solution results. Substituting the true value of $\mathbf{r}$, LS algorithm’s estimation of $\mathbf{r}$, EW-TLS algorithm’s estimation of $\mathbf{r}$, and AG-TLS algorithm’s estimation of $\mathbf{r}$ into Equations (1)–(3) yields the actual projected horizon ellipses, LS-projected horizon ellipses, EW-TLS-projected horizons, and AG-TLS-projected horizon ellipses in the image frame, respectively. These ellipses are shown in Figure 4 together with the measured horizon arc. Compared with the actual projected horizon ellipse, LS-projected horizon ellipses are generally offset and relatively small, which means that the distance from the camera to the planet estimated by the LS algorithm is larger than the actual value, and there will be a direction error. On the contrary, EW-TLS-projected horizon ellipses and AG-TLS-projected horizon ellipses can both basically fit the actual projected horizon ellipse. This is consistent with our previous analysis that the LS estimation introduces a non-negligible bias, while the TLS estimation eliminates the effect of the bias.

4.2. Performance under Different Horizon Arc Lengths

More Monte Carlo simulations are conducted to evaluate the algorithms’ performance under different lengths of detected horizon arc. The actual number of valid horizon points increases linearly with increasing arc length. A total of 5000 sets of independent simulations are conducted for each experimental condition. Various horizon measurements are produced by adding randomly generated Gaussian noise to the same scenario, and the results shown are the average values of these measurements. Errors of OPNAV are presented in the camera frame C. RMSE and MSTDR are used to evaluate the performance of three algorithms.

The results of the simulations are shown in Figure 5 and Figure 6. The RMSE of the three algorithms gradually decreased as the arc length was increased, which is consistent with the previous study of horizon-based OPNAV covariance [41]. At the same time, the gap between the RMSE and MSTDR between the LS algorithm and the TLS algorithms also decreased. As seen in Figure 5, both the EW-TLS algorithm and AG-TLS algorithm show great advantages over the LS algorithm in terms of RMSE when the horizon arc is less than 50 degrees, while the advantage of the two TLS algorithms is no longer obvious as the horizon arc length grows further. As shown in Figure 6, the mean residuals of the LS algorithm are larger than the STD when the horizon arc is less than 35 degrees. Even when the horizon arc is 95 degrees, the MSTDR of the LS algorithm is still larger than 20%. This will cause problems for the subsequent Kalman filter or other estimation frameworks, since they all assume that the measurements are zero mean error. In contrast, the two TLS algorithms maintain almost zero mean residuals as the horizon arc degree is changed. When the horizon arc length is long enough, the LS algorithm has almost no performance loss when compared with the two TLS algorithms. This indicates that as the measurement information is increased, the bias introduced by the LS algorithm becomes progressively smaller. As long as the captured horizon arc is longer than a certain value (e.g., 120 degrees), the LS algorithm can be used to obtain a faster calculation speed without any loss of accuracy; otherwise, it is recommended that the TLS algorithms be used to avoid mean residuals. As the horizon arc length increases, the MSTDR of the EW-TLS algorithm never exceeds 4%, while the MSTDR of the AG-TLS algorithm never exceeds 9%. From the perspective of RMSE, there is almost no significant difference in the performance of the two algorithms, and such a performance gap is negligible in practice.

Additionally, it is necessary to evaluate the execution times of the algorithms. Simulations are running multiple times at different horizon arcs and the average value is obtained as execution times. They are running in a computer with an intel${}^{\circledR }$ i7-10700 CPU and MATLAB${}^{\circledR }$ 2022a. The execution times include only the navigation algorithms, not the time for image processing, etc. The number of horizon points gradually increases with the horizon arc, and the relationship between the number of horizon points and execution time is shown in Figure 7. It can be seen that the number of measurement points dominates the execution time of the algorithms. The execution time of the EW-TLS algorithm increases from 1.3 to 8.1 ms with the increase of horizon points. Most execution time of the EW-TLS algorithm is consumed during the iterations, and a small portion is used to calculate the covariance of each measurement. The execution time of the LS algorithm never exceeds 0.23 ms, while the execution time of the AG-TLS algorithm never exceeds 0.43 ms, and such extra computation consumption is often acceptable. Compared with the EW-TLS algorithm, the AG-TLS algorithm brings a huge speed improvement with a slight loss of precision, which is very suitable for onboard applications with limited computing resources.

4.3. Opnav Performance in a Highly Elliptical Orbit

To fully evaluate the performance of the short-arc horizon-based OPNAV, a navigation scenario in a highly elliptical orbit is tested. The scientific exploration orbit is usually a highly elliptical orbit to meet the requirements of communication and observation at the same time. Elements of the orbit around Mars are listed in

Table 2. Elements of a highly elliptical orbit around Mars.

Parameter	Value
Semi-major axis	9500 km
Eccentricity	0.61
Inclination	87${}^{\circ }$
Argument of perige	320${}^{\circ }$
RAAN	64${}^{\circ }$
True anomaly	0${}^{\circ }$

. The test starts from the epoch 1 January 2022 07:49:00 (UTCG). In the frame P, the spacecraft’s initial position is (1.3212 $\times {10}^3$, 2.4107 $\times {10}^3$, −2.4834 $\times {10}^3$) km, and its initial velocity is (1.0999, 2.6389, 3.2311) km/s. The test covers 500 min, with perigee at the 0 and 470th minutes and the apogee at the 236th minute. The sampling interval for OPNAV is one minute. The parameters of the navigation camera are the same as those introduced earlier. The optical axis of the navigation camera is always aligned with the edge of Mars during the simulation, and the Mars horizon is on the diagonal of the image plane to ensure the best navigation performance. The influence of lighting conditions is not considered here to show the ideal situation, which may also be excluded by using an infrared camera, and the visible limb portion is always the entire one fitting within the camera’s FOV. The actual number of valid horizon points is always between 1700 and 1950. A total of 5000 sets of independent simulations are performed for each position, and the results shown are average values. Errors of OPNAV are presented in the camera frame C.

The Monte Carlo results of the LS algorithm, EW-TLS algorithm, and AG-TLS algorithm are shown in Figure 8. Gray points represent measurement residuals of Monte Carlo simulations, the blue solid line represents mean residuals of simulations, the red dash line represents 3$\sigma$ bounds of analytic covariance, and the black solid line represents 3$\sigma$ bounds of numeric covariance from simulations. Note that the covariance in Figure 8 is superimposed on the mean residuals. The three algorithms’ performances with the change in orbital altitude during the simulation are shown in Figure 9. The orbital altitude is raised from 300 to 12,000 km from perigee to apogee. The visible horizon arc length in the image plane gradually increases from 2 to 18 degrees as the orbital altitude is increased while, at the same time, the accuracy of each edge point gradually decreases. The horizon arc length and the edge point accuracy both affect the navigation accuracy and have opposite effects.

As shown in Figure 8b,c and Figure 9, there is a significant improvement in navigation accuracy with the increase of orbital altitude in the vicinity of the perigee. This shows that the increase in horizon arc length plays a dominant role in navigation accuracy at a low altitude. The navigation accuracy slowly decreases as the orbit altitude is further increased, which indicates that the decrease in the edge point accuracy at this time will affect the navigation accuracy to a greater extent. The STD of the navigation error ranges from 0.040 to 0.054 km in the $X_C$ and $Y_C$ axes and from 2.374 to 3.188 km in the $Z_C$ axis. The phase angle of the spacecraft does not show much effect in our simulation, which should be related to the small oblateness of Mars. As seen in Figure 8a and Figure 9, the LS algorithm shows a significant mean residual in the navigation results for all three axes at low orbital altitudes. The mean residuals of all three axes gradually decreases as the orbital altitude is further increased, and the influence in the $X_C$ and $Y_C$ axes can be neglected. However, a certain mean residual has been maintained in the $Z_C$ axis, and even at the apogee, there is still a mean residual of 1.16 km. While the mean residuals appear to be concentrated primarily along the $Z_C$ axis, it is important to note that the impact of residuals becomes more intricate when analyzed within the planetary frame. The $Z_C$ axis, which corresponds to the camera’s optical axis, is aligned with the edge of Mars. The attitude of the camera within the planetary frame varies as the spacecraft’s position changes. Consequently, the direction of the mean residuals will also be variable within the planetary frame. This means that the effect of the mean residuals should not be overlooked. On the contrary, the EW-TLS algorithm and AG-TLS algorithm both showed a better performance with near-zero mean residuals in different orbital altitudes. In general, the short-arc horizon-based OPNAV exhibits good accuracy throughout the entire orbital period, and it provides unbiased position estimation with the help of the TLS estimation.

In addition to orbital altitude, the performance of short-arc horizon-based OPNAV is affected by pointing errors. Ideally, the optical axis of the navigation camera should point to Mars’ horizon, and the horizon should remain on the diagonal of the image plane to ensure the best navigation performance, while it is difficult to achieve such precise pointing in practice. Taking the perigee as an example scenario, the performance of horizon-based OPNAV with different pointing errors on the $X_C$, $Y_C$, and $Z_C$ axes of the navigation camera are shown in Figure 10, Figure 11, and Figure 12, respectively.

As the pointing error increases, the horizon arc taken by the navigation camera becomes shorter and the OPNAV performance gradually decreases. The performance degradation of the LS algorithm is particularly significant compared to the EW-TLS algorithm and AG-TLS algorithm. When the pointing error on $X_C$ or $Y_C$ axis is -4 degrees, the LS algorithm’s positioning accuracy is worse than 10 km for both $X_C$ and $Y_C$ axes and 200km for $Z_C$ axis, while the EW-TLS and AG-TLS algorithms are better than 1 km for both $X_C$ and $Y_C$ axes and better than 25 km for $Z_C$ axis. With the pointing error around the $Z_C$ axis, the horizon shifts from the diagonal of the image plane to the horizontal or vertical direction. For the EW-TLS and AG-TLS algorithms, the positioning accuracy of $X_C$ and $Y_C$ axes decreases from 0.054 km to the worst, 0.108 km, and the positioning accuracy of $Z_C$ axis decreases from 3.188 km to the worst, 6.000 km. Meanwhile, for the LS algorithm, the positioning accuracy of $X_C$ and $Y_C$ axes ranges from 0.050 to 0.272 km, and the positioning accuracy of the $Z_C$ axis ranges from 10.214 to 25.052 km. The pointing error in the $Z_C$ axis has less impact on the optical navigation performance than the other two axes because the reduction of the horizon arc is smaller.

Pointing errors lead to the degradation of horizon-based OPNAV performance, and the EW-TLS algorithm and AG-TLS algorithm demonstrate a higher tolerance for pointing errors. The huge performance gap between the LS algorithm and TLS algorithms reveals the necessity of TLS estimation. Higher requirements should be placed in the pointing accuracy on $X_C$ and $Y_C$ axes than the $Z_C$ axis in the horizon-based OPNAV.

4.4. OPNAV Performance with Synthetic Images

To further validate the algorithm, the navigation performance after image processing is tested. Blender [42], an open source 3D computer graphics software program, is used for the synthesis of navigation images. Blender’s realism and convenience have made it popular in the field of OPNAV [43,44]. The specific steps of image synthesis and processing are as follows:

(1)

First, the basic simulation environment is built in Blender, including building Mars using its parameters and surface textures [45], setting up the navigation camera, and adding a light source to model the Sun. The parameters of the navigation camera are consistent with the previous section.

(2)

Next, automated rendering is implemented using the Python interface provided by Blender. Specifically, Python scripts are used to set the position and attitude of the sun and camera for each simulation scene, and the cycles engine is used to render photorealistic images. The camera is always pointed with its back to the Sun and toward the horizon of Mars to achieve good imaging.

(3)

Finally, MATLAB${}^{\circledR }$ is used for image processing of rendered images. The images are convoluted by a Gaussian kernel with an STD of 1.5 pixels to simulate the effect of defocusing, then the algorithm proposed in [13] is employed to achieve subpixel edge localization.

The example synthetic images for OPNAV around Mars are shown in Figure 13. It can be seen that as the orbit altitude decreases, the horizon no longer bends and gradually becomes straighter. OPNAV performance with synthetic images under different altitudes is shown in Figure 14. Due to the small number of simulations, the performance gap between the LS algorithm, EW-TLS algorithm, and AG-TLS algorithm is not significant. The navigation errors on the $X_C$ and $Y_C$ axes are between ± 0.1 km, while the navigation error on the $Z_C$ axis did not exceed ± 15 km. OPNAV performance with synthetic images at an orbital altitude of 300 km under different latitudes is shown in Figure 15. The navigation performance of the three algorithms seems to be similar on the $X_C$ and $Y_C$ axes. However, there is a noticeable difference in their performance on the $Z_C$ axis. The RMSE on the $Z_C$ axis of the LS algorithm is 10.2 km, while the RMSE on the $Z_C$ axis of the EW-TLS and AG-TLS algorithms are both 5.0 km. It can be seen that the LS algorithm showed certain mean residuals, which are well eliminated by both the EW-TLS algorithm and the AG-TLS algorithm. The performance of the proposed algorithms with synthetic images is not as good as in numerical simulation due to the limitation of simulation accuracy and edge positioning accuracy. However, the proposed algorithms can still bring considerable improvement to the performance of horizon-based OPNAV under a short-arc horizon scenario compared with the traditional algorithm.

5. Conclusions

To improve the performance of horizon-based OPNAV based on a short-arc horizon, a simplified measurement covariance model is proposed by further analyzing the propagation of measurement errors. Then, an unbiased solution with EW-TLS estimation is given by considering every element’s covariance of measurement equations. To further simplify this problem, an AG-TLS algorithm is proposed, which achieves a non-iterative solution by using the covariance of any one measurement as an approximation of all measurement covariances. Covariance analysis demonstrates that the covariances of the traditional Christian–Robinson algorithm with ordinary LS estimation and proposed algorithms are identical, and the difference between them lies in the possible estimation bias.

The numerical simulations under different horizon arc lengths show that when the horizon arc length is sufficiently long, both the traditional algorithm and the proposed algorithms achieve near-zero mean residual with consistent accuracy. With the shortening of the horizon arc length, the traditional algorithm will introduce a large mean residual that may even be larger than the STD, while the proposed algorithms both showed impressive advantages of near-zero mean residuals in the short-arc horizon scenario. Compared with the EW-TLS algorithm, the AG-TLS algorithm trades a negligible accuracy loss for a huge reduction in execution time and achieves a computing speed comparable to the traditional algorithm. Thus, the AG-TLS algorithm is very suitable for onboard applications with limited computing resources. Furthermore, as verified by the simulated navigation scenario with a highly elliptical orbit around Mars, a short-arc horizon can provide reliable position estimates for deep space exploration missions. Synthetic images with image processing further validate the effectiveness of the OPNAV method. The proposed algorithms extend the applicability of the horizon-based OPNAV method.