2.1 Coordinate systems and LOS pointing model

The OpNav system involves five coordinate systems: the inertial coordinate system, ${O_w} - {X_w}{Y_w}{Z_w}$, the body coordinate system, ${O_b} - {X_b}{Y_b}{Z_b}$, the camera coordinate system, ${O_c} - {X_c}{Y_c}{Z_c}$, the image coordinate system, ${O_{u^{\prime}v^{\prime}}} - xy$, and the pixel coordinate system, ${O_{uv}} - uv$. As shown in Figure 1, the ${O_w} - {X_w}{Y_w}{Z_w}$ is used to describe the spatial position and attitude of both the celestial body and the spacecraft. The conversion relationship between the ${O_b} - {X_b}{Y_b}{Z_b}$ and the ${O_c} - {X_c}{Y_c}{Z_c}$ is determined by the camera's mounting position. The ${O_{u^{\prime}v^{\prime}}} - xy$ and the ${O_{uv}} - uv$ are in the same image plane.

Figure 1. Relationship between the five coordinate systems

Based on the monocular camera, the OpNav system extracts the LOS vector from the camera to the planetary centroid. To simplify the calculation, the ${O_{u^{\prime}v^{\prime}}} - xy$ is introduced into the model. The ${O_{u^{\prime}v^{\prime}}} - xy$ and the ${O_{uv}} - uv$ lie in the same plane, but their origins and coordinate units differ. The point $p({x,\;y} )$ in the ${O_{u^{\prime}v^{\prime}}} - xy$ can be expressed in pixel coordinates as

(1)\begin{equation}\left\{ \begin{array}{@{}c} x= (u - u^{\prime})/{dx}\\ y= (v - v^{\prime})/{dy}\end{array}\right.\end{equation}

where $({u,\;v} )$ are the pixel coordinates of the navigation object, $({u^{\prime},\;v^{\prime}} )$ are the coordinates of the image plane centre, and $dx$ and $dy$ denote the number of pixels per unit of physical size along the respective coordinate axes. Based on the geometric relationship of the small-aperture imaging model, the relationship between the point $P({{X_i},\;{Y_i},\;{Z_i}} )$ in the${O_c} - {X_c}{Y_c}{Z_c}$ and the point $p({{x_i},\;{y_i}} )$ in the ${O_{u^{\prime}v^{\prime}}} - xy$ is given by

(2)\begin{equation}\left\{ {\begin{array}{@{}c} {{x_i} = f\dfrac{{{X_i}}}{{{Z_i}}}}\\ {{y_i} = f\dfrac{{{Y_i}}}{{{Z_i}}}} \end{array}} \right.\end{equation}

where f denotes the focal length of the camera, which is considered constant. Figure 2 illustrates the concept for estimating the optical observation LOS using image coordinates $({x,\;y} )$. If the ${Z_c}$ axis of the ${O_c} - {X_c}{Y_c}{Z_c}$ is assumed to coincide with the navigation camera's optical axis, the relationship between the focal length, a measured point on the image plane, and the LOS vector to the image source is given by (Christian, Reference Christian2010) the following:

(3)\begin{equation}e_i^C= \frac{1}{{\sqrt {x_i^2+ y_i^2+ {f^2}} }}\left[\begin{array}{@{}c@{}} { - {x_i}}\\ { - {y_i}}\\ f \end{array}\right]\end{equation}

where $e_i^C$ is the LOS vector in the ${O_c} - {X_c}{Y_c}{Z_c}$. Further, the LOS vector in the ${O_w} - {X_w}{Y_w}{Z_w}$ can be obtained by the rotation operation:

(4)\begin{equation}e_i^I = T_B^IT_C^Be_i^C\end{equation}

where $T_C^B$ is the conversion matrix from the camera frame to the body frame, and $T_B^I$ is the conversion matrix from the body frame to the inertial frame.

Figure 2. Imaging process and LOS direction model

2.1 Principle of OpNav

OpNav uses the camera to determine the spacecraft's position relative to the planet. However, determining the location using image data is an abstract process. Therefore, this section discusses the projection of a triaxial ellipsoid onto the image plane.

For a two-dimensional (2D) conic curve, the implicit equation is given by

(5)\begin{equation}Ax_i^2 +B{x_i}{y_i} +Cy_i^2 +D{x_i} +E{y_i} +F\textrm{ = 0}\end{equation}

where $({{x_i},\;{y_i}} )$ is the planetary edge in the image plane. If $4AC\mathrm{\ > }{B^2}$, the curve is an ellipse. Based on the relationship between a point $({{x_i},\;{y_i}} )$ in the image plane and the corresponding point $({{X_i},\;{Y_i},\;{Z_i}} )$ in the camera coordinate system, as described by Equation (2), the three-dimensional form of the model is

(6)\begin{equation}A{\left( {f\frac{{{X_i}}}{{{Z_i}}}} \right)^2} +B\left( {f\frac{{{X_i}}}{{{Z_i}}}} \right)\left( {f\frac{{{Y_i}}}{{{Z_i}}}} \right) +C{\left( {f\frac{{{Y_i}}}{{{Z_i}}}} \right)^2} +D\left( {f\frac{{{X_i}}}{{{Z_i}}}} \right) +E\left( {f\frac{{{Y_i}}}{{{Z_i}}}} \right) +F = 0\end{equation}

Rewritten into matrix form as

(7)\begin{equation}\left[ {\begin{array}{@{}ccc@{}} {{X_i}}& {{Y_i}}& {{Z_i}} \end{array}} \right]\left[ {\begin{array}{@{}ccc@{}} {A{f^2}}& {B{f^2}/2}& {Df/2}\\ {B{f^2}/2}& {C{f^2}}& {Ef/2}\\ {Df/2}& {Ef/2}& F \end{array}} \right]\left[ {\begin{array}{@{}c@{}} {{X_i}}\\ {{Y_i}}\\ {{Z_i}} \end{array}} \right] = 0\end{equation}

The ${P_i}$ and K are defined as

(8)\begin{equation}{\boldsymbol{P}_i} = {\left[ {\begin{array}{@{}ccc@{}} {{X_i}}& {{Y_i}}& {{Z_i}} \end{array}} \right]^T},\quad \boldsymbol{K} = \left[ {\begin{array}{@{}ccc@{}} {A{f^2}}& {B{f^2}/2}& {Df/2}\\ {B{f^2}/2}& {C{f^2}}& {Ef/2}\\ {Df/2}& {Ef/2}& F \end{array}} \right]\end{equation}

Then,

(9)\begin{equation}\boldsymbol{P}_i^T\boldsymbol{KP} = 0\end{equation}

It can be found that Equation (9) represents a conic surface. This implies that the image plane intersects the conic surface, resulting in an ellipse, which corresponds to the planetary 2D model in the image plane (as shown in Figure 3).

Figure 3. The process of planetary projection into the image plane

After deriving the conic surface from Equation (9), the relationship with the triaxial ellipsoid representing the observed planet is established. In the planet's inertial frame, any point ${P_i}$ on its surface can be written as

(10)\begin{equation}P_i^T\left[ {\begin{array}{@{}ccc@{}} {1/{a^2}}& 0& 0\\ 0& {1/{b^2}}& 0\\ 0& 0& {1/{c^2}} \end{array}} \right]{P_i} = P_i^T{A_P}{P_i} = 1\end{equation}

Based on the rotation matrix $T_C^P$ from the camera to the planetary inertial frame, Equation (10) is expressed in the camera frame as

(11)\begin{equation}p_i^TT_C^P{A_P}{(T_C^P)^T}{p_i} = p_i^TA{p_i} = 1\end{equation}

where ${p_i}$ is the point in the camera frame corresponding to ${P_i}$. As shown in Figure 4, the vector from the origin of the camera frame to the planetary point, ${p_i}$, is given by

(12)\begin{equation}{s_i} = {p_i} - r = k{e_i}\end{equation}

where ${e_i}$ is the unit vector of ${s_i}$ and k is the corresponding length. Substituting Equation (12) into Equation (11) yields

(13)\begin{equation}(e_i^TA{e_i}){k^2} +2(e_i^TAr)k +({r^T}Ar - 1) = 0\end{equation}

Figure 4. Geometry of the planetary surface as observed from the camera frame

If k is the only real solution, it implies that the point ${p_i}$ lies on the horizon. Therefore, the determining equation of Equation (13) is obtained as

(14)\begin{equation}4{(e_i^TAr)^2} - 4(e_i^TA{e_i})({r^T}Ar - 1) = 0\end{equation}

Rewritten as

(15)\begin{equation}e_i^T[{Ar{r^T}A - ({r^T}Ar - 1)A} ]{e_i} = 0\end{equation}

If $\boldsymbol{M} = Ar{r^T}A - ({r^T}Ar - 1)A$, substitute it into Equation (15) and multiply ${k^2}$ both sides of the equation by

(16)\begin{equation}s_i^T\boldsymbol{M}{s_i}\textrm{ = 0}\end{equation}

It has been demonstrated that all possible ${s_i}$ form a cone, and this conic surface is consistent with the one described by Equation (9). Therefore, if the 2D ellipse parameters are known, the position relative to the observed planet can be calculated using the projection relation. Note that A consists of the rotation matrix, which can be accurately determined using the spacecraft's attitude and the planetary ephemeris.

4.1 Outlier elimination

The process of outlier elimination involves dividing the edge points into two categories: inliers, which provide useful data for fitting the ellipse model, and outliers. Proximity information between edge points is used to capture the characteristics of the data, offering a reliable basis for outlier elimination. The Euclidean distance is employed to quantify proximity and facilitate the removal of outliers. The procedure is as follows:

1. 1) Construct a proximity information space ${\mathbb Z}$ to store the proximity data, ${{\rm Z}_{ij}} = ||{({x_i},\;{y_i}),\;({x_j},{y_j})} ||$, between point $({x_i},\;{y_i})$ and point $({x_j},{y_j})$;

2. 2) Count the number, ${k_i}$, of neighbouring points around the point $({x_i},\;{y_i})$ based on the proximity information space, ${\mathbb Z}$;

3. 3) Set a threshold, $\Theta$, to determine whether a point $({x_i},\;{y_i})$ is an outlier, assuming ${k_i} < \Theta$.

The threshold $\Theta$ significantly affects the performance of the outlier elimination process. When $\Theta$ is too small, sample points tend to form strong connections, making it difficult to detect outliers. Conversely, when $\Theta$ is too large, inliers within the sample space, $\mathrm{\mathbb{S}}$_, are more likely to be considered outliers. Therefore, careful selection of the threshold value is crucial. Six sets of comparative trials were conducted to determine the optimal value of $\Theta$. In these tests, the range of integer values for $\Theta$ was $\Theta \in [1\;,6]$. The results, presented in Figure 5, indicated that the best performance for outlier rejection was achieved when $\Theta = 3$. Almost all outliers were accurately detected, with no inliers excluded. Thus, the optimal threshold for the outlier elimination process was $\Theta = 3$.

Figure 5. Results of outlier elimination trials for different thresholds $\Theta$ ($\Theta \in [1\;,6]$)

4.2 Finding the optimal model using a hybrid genetic algorithm

The genetic algorithm is a stochastic search method that simulates biological evolution. It aims to identify individuals with high fitness through selection, crossover, and mutation operations. In this problem, the chromosome in the hybrid genetic algorithm represents the ellipse model, with the points on the ellipse functioning similarly to genes in a chromosome. Given that an ellipse model can be fitted using multiple pixel points, this approach innovatively treats the pixel points as genes, facilitating the solution of the ellipse model within the hybrid genetic algorithm framework. This strategy eliminates the need for traditional gene encoding and decoding processes, significantly improving computational efficiency. The complete procedure of the hybrid genetic algorithm is illustrated in Figure 6.

Figure 6. Overview of the method: orange arrows represent inputs, and green arrows represent outputs

In the current application, the initialisation operation randomly selects 180 edge points from the sample space, $\mathrm{\mathbb{S}}$, to construct the model population, M. Using these edge points, 30 chromosomes (ellipse models ${M_i} \in M$) are generated as the initial population, with each chromosome consisting of 6 genes (edge points $({x_i},\;{y_i}) \in \mathrm{\mathbb{S}}$). It is important to note that the outlier elimination procedure is performed before this step.

After the outlier elimination operation, it is assumed that the sample space $\mathrm{\mathbb{S}^{\prime}} \in \mathrm{\mathbb{S}}$ consists of normal data. A model, ${M_i}$, that includes a higher number of inliers typically corresponds to a lower error. Therefore, the performance of an ellipse model ${M_i}$ can be evaluated based on the number of inliers, ${N_{in}}({M_i})$. However, it is also important to consider the distribution of outliers around the model. The fitness evaluation function, $\rho$, is defined as follows:

(25)\begin{align} \begin{aligned} \rho & = \varDelta \ast \chi\\ & = \dfrac{{\left|{{N_{in}}({M_i}) - \dfrac{1}{n}\sum\limits_{i \le n} {{S_{in}}({M_i})} } \right|}}{{|\mathrm{\mathbb{S}} |}}exp \left( { - \dfrac{1}{2}{{\left[ {\dfrac{{\mathrm{{\cal F}}({M_i},\;\mathrm{\mathbb{S}})}}{\tau }} \right]}^2}} \right)\end{aligned}\end{align}

where $\varDelta$ is the reward function used to evaluate the number of inliers, and $\chi$ is the distance scale function employed to assess the distribution of outliers around the model. The maximum value of the $\chi$ function is 1, indicating that the current model represents an optimal solution.

During the selection operation, the probability of selecting an individual is proportional to its fitness. This strategy has the advantage of automatically identifying better-performing individuals while ensuring that no individuals in the population are discarded. The probability, P, of selecting an individual is given by

(26)\begin{equation}{P_i} = \frac{{|{a \cdot \rho_\varDelta^i + b \cdot \rho_\chi^i} |}}{{\sum_{j = 1}^n {\sqrt {{{({a \cdot \rho_\varDelta^j + b \cdot \rho_\chi^j} )}^2}} } }}\end{equation}

where a and b denote the factors of the reward function, ${\rho _\varDelta }$, and the distance scale function, ${\rho _\chi }$, respectively. The properties of the selection operation can be controlled by adjusting these two factors.

The crossover of two chromosomes is designed to evolve better individuals. Since edge points are treated as the genes of a chromosome, crossover is achieved by swapping points between individuals. Parental chromosomes are randomly mated to produce two offspring chromosomes. A crossover factor, $\alpha$, is introduced to control the crossover point between genes, with its value set in the range of $\alpha \in (0,\;1)$. In each crossover event, the parental chromosomes mate only once.

The original mutation method is a one-way, random operation, which often causes the algorithm to converge to locally optimal results. To improve the performance of the mutation operation, an annealing mutation operator is introduced, which operates in a two-way, probabilistic manner. The mutation factor, $\beta$, is used to control the mutation probability, with its value set in the range of $\beta \in (0,\;1)$. Mutation results are either accepted or rejected based on the improved Metropolis criterion. If ${\rho _{i + 1}} > \;{\rho _i}$, the probability of accepting the new variant result is 1. Otherwise, the acceptance probability is determined using the distance scale function, $\chi$, of the outliers. The improved Metropolis criterion is defined as

(27)\begin{equation} \textrm{Improved}\;\textrm{Metropolis}:= \left\{ \begin{array}{@{}cc@{}} {1,}& {{\rho_{i + 1}}> {\rho_i}}\\ {exp \left( { - \dfrac{1}{2}{{\left( {\dfrac{{{\rho_{i + 1}}}}{{{\rho_i}}}} \right)}^2}} \right),}& {{\rho_{i + 1}} < {\rho_i}} \end{array} \right.\end{equation}

5.1 $\alpha ,\beta$ test

Before evaluating the performance of the algorithm, it is necessary to determine the optimal crossover parameter, $\alpha$, and the annealing mutation parameter, $\beta$. The crossover parameter, $\alpha$, controls the mating operations between two chromosomes, enabling the population to evolve into better individuals. The annealing mutation parameter, $\beta$, is introduced to prevent the algorithm from converging to a local optimum.

The general performance analysis of the method is based on the following terms. The inlier function, G, calculated using the parameters $\alpha ,\beta$, is denoted as

(28)\begin{equation}{G_t}(\alpha ,\beta ):\; = N_{in}^t({M_{\alpha ,\beta }})\end{equation}

The average performance curve, L, is defined as

(29)\begin{equation}{L_t}(\alpha ,\beta ):\; = \frac{{{G_t}(\alpha ,\beta )}}{m} = \frac{1}{m}\sum\limits_{i = 1}^m {N_{in}^t({M_i})}\end{equation}

${L_t}(\alpha ,\beta )$ represents the average number of inliers in the evolutionary process of the algorithm at time t. It reflects the average performance of the algorithm during its evolution. The parameters $\alpha$ and $\beta$ in the algorithm exhibit a degree of coupling, and it is important to consider the impact of this pair of parameters on the algorithm's performance. Therefore, the relative goodness curve, $\mathrm{{\cal L}}$, is defined to represent the performance graph:

(30)\begin{equation}\mathrm{{\cal L}}(\alpha ,\beta ):\; = \frac{{L(\alpha ,\beta ) - \mathop {\min }\limits_{\forall \alpha ^{\prime},\beta ^{\prime}} G(\alpha ^{\prime},\beta ^{\prime})}}{{\mathop {\max }\limits_{\forall \alpha ^{\prime},\beta ^{\prime}} G(\alpha ^{\prime},\beta ^{\prime}) - \mathop {\min }\limits_{\forall \alpha ^{\prime},\beta ^{\prime}} G(\alpha ^{\prime},\beta ^{\prime})}}\end{equation}

In the test, the parameter ranges for $\alpha$ and $\beta$ are set as $(0,0 \cdot 95]$. The performance of $\mathrm{{\cal L}}$ (as shown in Figure 7) improves significantly after several iterations. Specifically, the maximum value of $\mathrm{{\cal L}}$ increases from 0 ⋅ 3824 in the fifteenth generation to a peak value of 0 ⋅ 6384. As a result, the optimal values for the parameters were found to be $\alpha = 0 \cdot 55$ and $\beta \textrm{ = 0} \cdot \textrm{15}$.

Figure 7. Relative goodness tests are performed to determine a pair of coupled parameters: the crossover parameter, $\alpha$, and the annealing mutation parameter, $\beta$

5.2 Ellipse fitting test

The edge points of the ellipse are used to evaluate the performance of the algorithm. To test its accuracy and robustness, various levels of Gaussian noise are introduced into the data. The standard deviations, $\varDelta$, of the Gaussian noise are set to 0 ⋅ 3, 1, 1 ⋅ 5, and 2, respectively. After the iterative process, the results are shown in Figures 8–11. The experimental results demonstrated that the proposed method achieved superior accuracy and robustness. As noted previously, RANSAC tends to converge to a local optimum as noise levels increase. In contrast, the proposed method consistently and accurately extracted the centroid of the ellipse model.

Figure 8. RANSAC and the proposed method utilise edge points to fit the ellipse model $({\varDelta = 0 \cdot 3} )$

Figure 9. RANSAC and the proposed method utilise edge points to fit the ellipse model $({\varDelta \textrm{ = 1}} )$

Figure 10. RANSAC and the proposed method utilise edge points to fit the ellipse model $({\varDelta \textrm{ = 1} \cdot \textrm{5}} )$

Figure 11. RANSAC and the proposed method utilise edge points to fit the ellipse model $({\varDelta = 2} )$

The number of inliers involved in the ellipse model is shown in Figure 12, where the total number of edge points is 506. The results indicated that the outlier elimination operation effectively mitigated the impact of noise. Furthermore, the ellipse models computed by the proposed method incorporated nearly all valuable edge points. Consequently, the precision and robustness of this approach outperformed those of RANSAC.

Figure 12. Comparison of the number of inliers between RANSAC and the proposed method under different levels of Gaussian noise

The statistical results of the comparison test are presented in Table 1. The minimum fitting error for RANSAC was 0 ⋅ 671582 pixels, while the maximum error reached 2 ⋅ 462185 pixels. In contrast, the maximum fitting error of the proposed method was 0 ⋅ 190068 pixels, demonstrating that the proposed approach for deep space OpNav yielded superior performance in celestial centroid extraction. However, the accuracy of centroid extraction does not directly reflect the overall performance of the OpNav system. To assess this, the LOS direction for OpNav must be computed according to Equation (3).

Table 1. Errors of the planetary centroid

For this evaluation, the camera parameters are assumed to be a focal length of $f = 200$ mm, a field of view of $\vartheta = 3 \cdot {5^\circ }$, and scale factors of $dx = 78 \cdot 4$ pixel/mm and $dy = 68 \cdot 7$ pixel/mm. The LOS pointing accuracy for OpNav is shown in Table 2. The results indicated that the LOS direction errors of the proposed method were less than $0 \cdot 902536 \times {10^{ - 5}}$ rad under noise $\varDelta = 0 \cdot 3$, and less than $1 \cdot 66826 \times {10^{ - 5}}$ rad under noise $\varDelta = 2$, both of which are smaller than the corresponding errors of the RANSAC method.

Table 2. Errors of the LOS direction

Next, the centroid extraction errors of the two methods were compared under varying levels of noise. Both approaches used edge points with different noise levels to fit ellipse models, and the results are presented in Figure 13. As the standard deviation, $\varDelta$_, of the noise increased, the fitting errors for both methods also increased. However, it was evident that the fitting error of the proposed method gradually decreased to a small value as the number of iterations grew. Due to the diverse population in the algorithm, the proposed method consistently demonstrated smaller fitting errors compared to RANSAC. In conclusion, the proposed method not only accurately extracted the planetary centroid but also effectively mitigated the impact of noise.

Figure 13. Errors in the ${x_0}$ and ${y_0}$ coordinates of centroid extraction under different noise levels

To simulate realistic OpNav scenarios, celestial images were used to evaluate the performance of the two algorithms. These images were generated using Celestia, an authoritative software developed by NASA. Three images of a planet from different scenes were considered, with varying levels of Gaussian noise ($\varDelta = 5,\;10,\;25$) were added to the images. The results of the two methods are shown in Figures 14–16. The position of the celestial centroid is located at the centre of the image, marked by the white spot. It can be observed that the proposed method accurately extracted the celestial centroid, even under different noise levels. Therefore, the effectiveness of the proposed method was demonstrated, making it suitable for application in deep space OpNav.

Figure 14. Results of RANSAC (green) and the proposed method (red) for fitting planetary images. The white spot is the precise centroid $({\varDelta = 5} )$

Figure 15. Results of RANSAC (green) and the proposed method (red) for fitting planetary images. The white spot is the precise centroid $({\varDelta = 10} )$

Figure 16. Results of RANSAC (green) and the proposed method (red) for fitting planetary images. The white spot is the precise centroid $({\varDelta = 25} )$