2. FEATURE LINE CORRESPONDENCES

In order to estimate motion parameters accurately, feature lines on the plane of the planetary surface need to be extracted from their Two-Dimensional (2D) images and matched to the same feature lines with known position in the database during the descent phase. This section focuses on feature line extraction and matching algorithms.

First of all, Akinlar and Topal's (Reference Akinlar and Topal2011) extraction algorithm, EDLine (Edge Drawing Line), is applied to our proposed algorithm to extract feature lines from landing images. As this algorithm is obviously faster than more traditional algorithms, it can have a better real-time performance (Burns et al., Reference Burns, Hanson and Riseman1986; Grompone von Gioi et al., Reference Grompone von Gioi, Jakubowicz, Morel and Randall2008; Etemadi, Reference Etemadi1992). As the method can directly extract the directions of feature lines and points on the feature lines, the outputs of the algorithm are consistent with the input of our visual navigation algorithm. Secondly, the algorithm presented in Zhang and Koch (Reference Zhang and Koch2013) is selected to match the extracted feature lines from the landing image with the feature lines in a database as it takes advantage of the Line Band Descriptor (LBD), which is robust to the rotation, scaling and brightness of image to finish feature line matching. In addition, because the extracted feature lines of EDLine are fitted by using a single pixel edge from an edge drawing, it helps that the main direction of Zhang and Koch's method is found exactly. Hence, Akinlar and Topal's (Reference Akinlar and Topal2011) extraction algorithm is combined with Zhang and Koch's (Reference Zhang and Koch2013) matching algorithm to finish feature line matching in this paper.

In this paper, the landing images taken by spacecraft Spirit at altitudes of approximately 1,986 m and 1,690 m and spacecraft Curiosity at altitudes of approximately 6,654 m and 6,221 m are extracted and matched by using the algorithm mentioned above (Resolution of landing images are 1024 × 1024 pixels). The matching results are shown in Figures 1 and 2. In the matching result of spacecraft Spirit, 106 and 194 feature lines from two landing images are extracted respectively in 920 ms and 1,277 ms, and 33 pairs of feature lines are matched in 97 ms. Meanwhile, the results matched by using landing images from spacecraft Curiosity show that 45 and 55 feature lines from two landing images are extracted respectively in 321 ms and 337 ms, and 20 pairs of feature lines are matched in 31 ms. In conclusion, because the numbers of extracted feature lines and matched feature lines are respectively more than 33 and 20, they provide a prerequisite for the development of the VN-FLC. As the time of feature line extraction and feature line matching is less than 1,277 ms and 97 ms respectively, the image processing has better real-time performance than Scale-Invariant Feature Transform (SIFT) in which the time of feature point extraction and feature line matching is more than 6 s and 3 s, respectively. (Matching conditions: CPU: Intel(R) Core(TM) i7-4558U CPU @ 2.80GHz, 4GB, Visual Studio 2013, Opencv2.4.10, ARPACK (ARnoldi PACKage), BIAS-2.8.0 (Basic Linear Algebra Subprograms), CLAPACK-3.1.1(C Language Interface of Linear Algebra PACKage), LAPACK (Linear Algebra PACKage) and SuperLU (Supernodal LU)).

Figure 1. Feature line extraction and matching based on the images of spacecraft Spirit.

Figure 2. Feature line extraction and matching based on the images of spacecraft Curiosity.

3. CONSTRAINT EQUATIONS

Here, it is assumed that the planet surface is a plane and all feature lines on the planet surface are in the same plane. The coordinate systems used in this paper are then introduced, and the nonlinear equations are deduced according to coordinate transformation and geometric constraints.

3.1. Coordinate Systems

The definition of coordinate systems is important for clearly explaining the coordinate transformation and geometric constraints. In this paper, the coordinate systems are built as Figures 3 and 4.

Figure 3. A geometry sketch map of several common coordinate systems.

Figure 4. A geometry sketch map of model coordinate systems.

3.1.1. The mass centre coordinate system O _w − X _wY _wZ _w

O _w is the mass centre of the targeted planet; O _wX _w directs to the intersection between the equatorial plane and the ecliptic plane; the Z _w-axis is perpendicular to the equatorial plane; the X _w-axis, Y _w-axis and Z _w-axis meet the right-hand corkscrew rule.

3.1.2. The cartographic coordinate system O _s − X _sY _sZ _s

O _s is the landing site; the direction of O _sZ _s is the direction of vector from the mass centre of the targeted planet to the landing site; O _sX _s is the vector along the tangent line of the meridian of the landing site to the south pole; X _s-axis, Y _s-axis and Z _s-axis meet the right-hand corkscrew rule.

3.1.3. The camera coordinate system O _c − X _cY _cZ _c

O _c is the optical centre of the navigation camera; Z _c-axis is defined as the optic axis and the X _c-axis and Y _c-axis are parallel to the horizontal axis and vertical axis of image respectively.

3.1.4. The model coordinate system O _m − X _mY _mZ _m

In this paper, L_i = (V_i, P_i) denotes the 2D feature lines on the planet surface, in which V_i is defined as the unit vector of the direction of the feature lines and P_i is defined as a point on these feature lines. l_i = (s_i, e_i) denotes the projection of L_i on the 2D image plane, in which s_i and e_i are defined as the endpoints of l_i. Let O _m − X _m Y _m Z _m be the model coordinate system, in which V₀ is defined as the X _m-axis (the projection length of L₀ is the longest from {L_i}), the Z _m-axis is parallel to the Z _s-axis, and O _m coincides with the origin O _s of the cartographic coordinate system O _s − X _s Y _s Z _s.

3.2. The perspective-3-line constraint equation

In this paper, it is assumed that L_i = (V_i, P_i) and l_i = (s_i, e_i) are known. ${\bi R}_{{\bi s}}^{{\bi m}}$ denotes that the direction vector is rotated from the cartographic coordinate system to the model coordinate. Therefore, ${\bi V}_{0}^{{\bi m}} ={\bi R}_{{\bi s}}^{{\bi m}} {\bi V}_{0}^{{\bi s}}=\lsqb 1\comma \; 0\comma \; 0\rsqb ^{T}$. It can be shown that ${\bi R}_{{\bi s}}^{{\bi m}}$ can be written as $\lsqb \lpar {\bi V}_{0}^{{\bi s}} \rpar ^{T}\semicolon \; {\bf \varepsilon}_{1}^{T} \semicolon \; {\bf \varepsilon}_{2}^{T}\rsqb $ where ε₁ and ε₂ are mutually orthogonal and they are solutions from $\lpar {\bi V}_{0}^{{\bi s}}\rpar ^{T}{\bi X}=0$. Similarly, the directions of the remaining feature lines are transformed from the cartographic coordinate system to the model coordinate by ${\bi R}_{{\bi s}}^{{\bi m}} $, such that ${\bi V}_{{\bi i}}^{{\bi m}} ={\bi R}_{{\bi s}}^{{\bi m}} {\bi V}_{{\bi i}}^{{\bi s}}$.

The plane that passes through the origin O _c of the camera coordinate system and l_i is represented as Π_i. The normal of Π_i is defined as n_i. Therefore, the rotation matrix ${\bi R}_{{\bi m}}^{{\bi c}}$ by which ${\bi V}_{0}^{{\bi m}}$ in the cartographic coordinate system is transformed to ${\bi V}_{0}^{{\bi c}}$ in the camera coordinate system must meet the constraint that ${\bi n}_{0}^{T} {\bi R}_{{\bi m}}^{{\bi c}} {\bi V}_{0}^{{\bi m}} =0$. Due to this constraint, ${\bi R}_{{\bi m}}^{{\bi c}}$ can be decomposed as:

(1)$$\eqalign{{\bi R}_{{\bi m}}^{{\bi c}} &={\bi R}{\bi R}_{{\bi y}} \lpar \gamma \rpar {\bi R}_{{\bi x}} \lpar \beta \rpar \cr & = \left[{\matrix{ r_1 & r_2 & r_3 \cr r_4 & r_5 & r_6 \cr r_7 & r_8 & r_9}}\right]\left[{\matrix{\cos \gamma & 0 & -\sin \gamma \cr 0 & 1 & 0 \cr \sin \gamma & 0 & \cos \gamma \cr }}\right]\left[{\matrix{ 1 & 0 & 0 \cr 0 & \cos \beta & -\sin \beta \cr 0 & \sin \beta & \cos \beta \cr}} \right]} $$

in which R is an arbitrary orthogonal identity matrix whose second column equals n₀, R_y(λ) means a rotation whose angle around the Y-axis is equal to γ, and R_x(β) means a rotation whose angle around the X-axis is equal to β.

${\bi R}_{{\bi s}}^{{\bi c}}$ denotes that the direction vector is rotated from the cartographic coordinate system to the camera coordinate system. It is expressed as:

(2)$${\bi R}_{{\bi s}}^{{\bi c}} ={\bi R}_{{\bi m}}^{{\bi c}} {\bi R}_{{\bi s}}^{{\bi m}} $$

In this condition, ${\bi R}_{{\bi m}}^{{\bi c}}$ can be determined by two unknown variables λ and β. Owing to the geometrical constraint that the normal n_i of the plane Π_i can be perpendicular to all lines in this plane, it meets the constraint that:

(3)$${\bi n}_{{\bi i}}^T {\bi R}_{{\bi m}}^{{\bi c}} {\bi V}_{{\bi i}}^{{\bi m}} ={\bi n}_{{\bi i}}^T {\bi R}_{{\bi s}}^{{\bi c}} {\bi V}_{{\bi i}}^{{\bi s}} \equiv 0 $$

Hence, the other two constraints can be expressed as:

(4)$$\left\{\matrix{{\bi n}_1 \cdot {\bi V}_1^{{\bi c}} ={\bi n}_1^T {\bi R}_{{\bi m}}^{{\bi c}} {\bi V}_1^{{\bi m}} =0 \cr {\bi n}_2 \cdot {\bi V}_2^{{\bi c}} ={\bi n}_2^T {\bi R}_{{\bi m}}^{{\bi c}} {\bi V}_2^{{\bi m}} =0\cr}\right. $$

Let ${\bi n}_{{\bi i}}^{T} =\lsqb {x}^{\prime}_{i} \comma \; {y}^{\prime}_{i} \comma \; {z}^{\prime}_{i} \rsqb $ and ${\bi V}_{{\bi i}}^{{\bi m}} =\lsqb a_{i}^{m} \comma \; b_{j}^{m} \comma \; 0\rsqb ^{T}$. By substituting Equation (1) into Equation (4), the constraints can be deduced as:

(5)$$\left\{\matrix{\sigma_1 \cos \gamma +\sigma_2 \sin \gamma +\sigma_3 =0 \cr \sigma_4 \cos \gamma +\sigma_5 \sin \gamma +\sigma_6 =0 \cr}\right. $$

in which

$$\left\{\matrix{ \sigma_1 =a_1^m {x}^{\prime}_1 +b_1^m {z}^{\prime}_1 \sin \beta \hfill \cr \sigma_2 =b_1^m {x}^{\prime}_1 -a_1^m {z}^{\prime}_1 \hfill \cr \sigma_3 =b_1^m {y}^{\prime}_1 \cos \beta \hfill \cr \sigma_4 =a_2^m {x}^{\prime}_2 +b_2^m {z}^{\prime}_2 \sin \beta \hfill \cr \sigma_5 =b_2^m {x}^{\prime}_2 -a_2^m {z}^{\prime}_2 \hfill \cr \sigma_6 =b_2^m {y}^{\prime}_2 \cos \beta \hfill \cr}\right.$$

By solving Equation (5), cos γ and sin γ are written as:

(6)$$\cos \gamma =\displaystyle{{\sigma_2 \sigma_6 -\sigma_3 \sigma_5 }\over{\sigma_1 \sigma_5 -\sigma_2 \sigma_4 }}\comma \; \sin \gamma =\displaystyle{{\sigma_3 \sigma_4 -\sigma_1 \sigma_6 }\over{\sigma_1 \sigma_5 -\sigma_2 \sigma_4 }} $$

By the geometrical constraint that cos ²γ + sin²γ = 1, we can get:

(7)$$\lpar \sigma_2 \sigma_6 -\sigma_3 \sigma_5 \rpar ^2+\lpar \sigma_3 \sigma_4 -\sigma _1 \sigma_6 \rpar ^2=\lpar \sigma_1 \sigma_5 -\sigma_2 \sigma_4 \rpar ^2 $$

Similarly, due to cos ²β + sin ²β = 1, we obtain a perspective-3-line constraint equation:

(8)$$f\lpar r\rpar =a_4 r^4+a_3 r^3+a_2 r^2+a_1 r+a_0 =\sum\limits_{k=0}^4 {a_k r^k} =0 $$

in which r denotes sin β.

The perspective-3-line constraint equation can be applied to estimate the four candidate solutions of sin β, when three arbitrary feature lines on the same plane are matched with their projections.

3.3 The Perspective-n-Line (PnL) constraint equations

As longer feature line projections are less affected by noise, the feature line whose projection is the second longest is selected as a guideline, and the feature line whose projection is the longest has been selected as the X _m-axis. We can construct n-2 constraint equations together with the rest of the feature lines such as Equation (8):

(9)$$\left\{\matrix{f_1 \lpar r\rpar =\sum_{k=0}^4 {a_{1k} r^4} =0 \hfill \cr f_2 \lpar r\rpar =\sum_{k=0}^4 {a_{2k} r^4} =0 \hfill \cr \vdots \cr f_{n-2} \lpar r\rpar =\sum_{k=0}^4 {a_{\lpar n-2\rpar k} r^4} =0 \hfill \cr}\right. $$

In order to get the four candidate solutions of sin β, the Perspective-n-Line (PnL) must be solved. However, the calculation process is difficult. Aiming at solving Equation (9) easily, an objective function is proposed in Section 4.

4. MOTION ESTIMATION

In this section, the four candidate solutions of sin β are solved based on the nonlinear Equation (9). Then, the four candidate solutions for the lander's pose are estimated according to the four candidate solutions of sin β. Finally, the unique position and attitude of the lander are selected.

4.1. Candidate solutions

Firstly, for solving the nonlinear Equation (9), an objective function Equation (10) is constructed to pick the optimal solutions by least square methods.

(10)$$F=\sum_{i=1}^{n-2} {f_i^2 \lpar r\rpar } $$

In order to obtain the candidate solutions, the local minima of Equation (10) is computed. Its derivative is calculated as:

(11)$$F^{\prime}=\sum_{i=1}^{n-2} {f_i \lpar r\rpar {f}^{\prime}_i \lpar r\rpar } $$

The minima is determined by computing the roots of its derivative Equation (11). Four roots are counted easily by the eigenvalue method and they are candidate solutions of sin β.

Then, let ${\bi P}_{{\bi i}}^{{\bi s}}$ be a point on the L_i in the cartographic coordinate system and t be the coordinate of O _c in the cartographic coordinate system. Hence, we can get:

(12)$${\bi n}_{{\bi i}}^T {\bi R}_{{\bi s}}^{{\bi c}} \lpar {\bi P}_{{\bi i}}^{{\bi s}} -{\bi t}\rpar =0 $$

Substituting Equations (1), (2) and candidate solutions of sin β into Equation (4) and Equation (12), the following equation is obtained:

(13)$$\left\{\matrix{{\bi n}_{{\bi i}}^T {\bi R} {\bi R}_{{\bi y}} \lpar \gamma \rpar {\bi R}_{{\bi x}} \lpar \beta_j\rpar {\bi R}_{{\bi s}}^{{\bi m}} {\bi V}_{{\bi i}}^{{\bi s}} =0 \hfill \cr {\bi n}_{{\bi i}}^T \lsqb {\bi R} {\bi R}_{{\bi y}} \lpar \gamma \rpar {\bi R}_{{\bi x}} \lpar \beta_j\rpar {\bi R}_{{\bi s}}^{{\bi m}} {\bi P}_{{\bi i}}^{{\bi s}} -\bar{{\bi t}}\rsqb =0 \hfill \cr}\right. $$

in which $\bar{{\bi t}}={\bi R}_{{\bi s}}^{{\bi c}} {\bi t}=\lsqb \matrix{{x^{c}} & {y^{c}} & {z^{c}}}\rsqb ^{T}$.

Let ${\bi n}_{{\bi i}}^{T} {\bi R}=\lsqb \bar{x}_{i} \comma \; \bar{y}_{i} \comma \; \bar{z}_{i}\rsqb $, ${\bi R}_{{\bi x}} \lpar \beta_{j}\rpar {\bi R}_{{\bi s}}^{{\bi m}} {\bi V}_{{\bi i}}^{{\bi s}} =\lsqb \matrix{{a_{ij} } & {b_{ij} } & {c_{ij} }}\rsqb ^{T}$ and ${\bi R}_{{\bi x}} \lpar \beta_{j}\rpar {\bi R}_{{\bi s}}^{{\bi m}} {\bi P}_{{\bi i}}^{{\bi s}} =\lsqb \matrix{{x_{ij}^{m} } & {y_{ij}^{m} } & {z_{ij}^{m} }}\rsqb ^{T}$. A linear equation is obtained:

(14)$${\bi M}_{{\bi j}} {\bi Q}=0 $$

in which

$$\eqalign{{\bi M}_{{\bi j}} & =\left[{\matrix{{\bar {x}_i a_{ij} +\bar {z}_i c_{ij} } & {\bar {z}_i a_{ij} -\bar {x}_i c_{ij} } & 0 & 0 & 0 & {\bar {y}_i b_{ij}} \cr {\bar {x}_i x_{ij}^m +\bar {z}_i z_{ij}^m } & {\bar {z}_i x_{ij}^m -\bar {x}_i z_{ij}^m } & {-{x}^{\prime}} & {-{y}^{\prime}} & {-{z}^{\prime}}& {\bar {y}_i y_{ij}^m}}} \right]\cr {\bi Q}&=\left[{\matrix{{\cos \lambda } & {\sin \lambda } & {x^c} & {y^c} & {z^c} & 1 }}\right]^T}$$

By solving Equation (14) with SVD, the candidate solutions of the rotation angle γ and the rotated translation vector $\bar{{\bi t}}$ are obtained. Then, the candidate solution of the translation vector is calculated as ${\bi t}=\lpar {\bi R}_{{\bi s}}^{{\bi c}} \rpar ^{T}\bar{{\bi t}}$.

4.2. Determination of lander's attitude and position

The candidate solutions of ${\bi R}_{{\bi s}}^{{\bi c}}$ and t are calculated from the linear system Equation (14). The solutions which are affected by noise in the data extracted from images may be not a normalisation. Hence, the solutions should be normalised by a standard Three-Dimensional (3D) alignment scheme (Umeyama, Reference Umeyama1991).

After the normalisation, the orthogonal error E _er of each candidate solution is computed as:

(15)$$E_{er} =\sum_{i=1}^n \lpar {\bi n}_{{\bi i}}^T {\bi R}_{{\bi s}}^{{\bi c}} {\bi V}_{{\bi i}}\rpar ^2 $$

As the results of multiple experiments show that their orthogonal errors are less than10⁻⁴, the orthogonal errors of candidate solutions which are larger than 10⁻² are deleted. The other m ₁ candidate solutions are retained.

Then, ${\bi P}_{{\bi ni}}^{{\bi s}}$ is defined as the closest point on the feature line L_i to the origin of the cartographic coordinate. $P_{ni}^{s}$ can be shown as:

(16)$${\bi P}_{{\bi ni}}^{{\bi s}} ={\bi P}_{{\bi i}}^{{\bi s}} -\lpar {\bi P}_{{\bi i}}^{{\bi s}} \cdot {\bi V}_{{\bi i}}^{{\bi s}}\rpar {\bi V}_{{\bi i}}^{{\bi s}} $$

${\bi P}_{{\bi ni}}^{{\bi c}}$ and ${\bi P}_{{\bi i}}^{{\bi c}}$ on the feature line L_i in the camera coordinate are obtained by using the candidate solution of the rotation matrix and the translation vector from Equation (14), and are projected onto the interpretation plane of lines $\lcub \bar{{\bi P}}_{{\bi ni}}^{{\bi c}}\rcub $ and $\lcub \bar{{\bi P}}_{{\bi i}}^{{\bi c}}\rcub $ as follows:

(17)$$\left\{\matrix{{\bi P}_{{\bi i}}^{{\bi c}} ={\bi R}_{{\bi s}}^{{\bi c}} \lpar {\bi P}_{{\bi i}}^{{\bi s}}-{\bi t}\rpar \hfill \cr \bar{{\bi P}}_{{\bi i}}^{{\bi c}} ={\bi P}_{{\bi i}}^{{\bi c}} -\lpar {\bi P}_{{\bi i}}^{{\bi c}} \cdot {\bi n}_{{\bi i}}\rpar {\bi n}_{{\bi i}}}\right. $$

(18)$$\left\{\matrix{{\bi P}_{{\bi ni}}^{{\bi c}} ={\bi R}_{{\bi s}}^{{\bi c}} \lpar {\bi P}_{{\bi ni}}^{{\bi s}} -{\bi t}\rpar \hfill \cr \bar{{\bi P}}_{{\bi ni}}^{{\bi c}}={\bi P}_{{\bi ni}}^{{\bi c}}-\lpar {\bi P}_{{\bi ni}}^{{\bi c}} \cdot {\bi n}_{{\bi i}}\rpar {\bi n}_{{\bi i}}}\right. $$

Following Figure 5, we define a re-projection residual E _re as follows:

Figure 5. Illustration of observed feature line model and re-projected feature line model on the projection plane.

(19)$$E_{re} =\vint_0^{\ell_{i}} D_i^2 \lpar s\rpar d\lpar s\rpar =\sum_{i=1}^n \displaystyle{{\ell_i }\over{3}}\lpar D_{is}^2 +D_{is} D_{ie} +D_{ie}^2\rpar $$

in which ℓ_i is the length of the projection l_i on the image plane and D _is and D _ie denote the distance between the observed feature line endpoints and the re-projected feature line endpoints.

Substituting Equations (17) and (18) into Equation (19), we can obtain the re-projection residual E _re of each candidate solution and select the solution with the smallest re-projection residual as the lander's attitude and position.

Lastly, in order to improve accuracy, the lander's attitude and position are calibrated by using Kumar and Hanson's (Reference Kumar and Hanson1994) algorithm.

5. APPLICATION OF EXTENDED KALMAN FILTER

In this section, in order to improve the navigation accuracy, the position and attitude of the lander are estimated by an extended Kalman filter.

5.1. Dynamic equations

Following Figure 3, the dynamic model of the spacecraft's landing orbit is derived as follows. Firstly, the universal dynamic model of the lander's position in the cartographic coordinate system O _s − X _s Y _s Z _s is established (Li and Cui, Reference Li and Cui2008):

(20)$$\eqalign{\dot{{\bi r}} & ={\bi v} \cr \dot{{\bi v}} & =-2{\bf \omega}_{{\bi a}} \times \dot{{\bi r}}+{\bf \omega}_{{\bi a}} \times \lpar {\bf \omega}_{{\bi a}} \times {\bi r}\rpar +{\bi U}+\Delta {\bi F}} $$

in which r is defined as the vector from the centroid of the planet to the mass centre of the lander, ΔF is defined as the solar radiation pressure, the third body gravity and so on, v denotes the velocity of the lander, ω_a represents the angular velocity of the cartographic coordinate system relative to the mass centre coordinate system and U means the gravitational acceleration of the planet.

Similarly, the dynamic model of the lander's attitude at landing can be developed by the angular velocity ω of the lander relative to the cartographic coordinate system and the angular velocity ω_a:

(21)$$\left[\matrix{{\dot {\beta }} \cr {\dot {\gamma }} \cr {\dot {\alpha}}}\right]={\bi R}_{{\bi b}} \left[\matrix{ {\omega_x } \cr {\omega_y } \cr {\omega_z }}\right]-{\bi R}_{{\bi a}} \left[\matrix{{\omega_{ax} } \cr {\omega_{ay} } \cr {\omega_{az} }}\right]$$

in which

$${\bi R}_{{\bi a}} =\left[\matrix{ {\displaystyle{{\sin \alpha }\over{\cos \gamma}}} & {\displaystyle{{\cos \alpha }\over{\cos\gamma }}} & 0 \cr {\cos \alpha } & {-\sin \alpha } & 0 \cr {\displaystyle{{\sin \gamma \sin \alpha }\over{\cos \gamma }}} & {\displaystyle{{\sin \gamma\cos \alpha }\over{\cos \gamma }}} & {-1} \cr} \right]\comma \; \quad {\bi R}_{{\bi b}} =\left[\matrix{ 1 & {\displaystyle{{\sin \gamma \sin \beta }\over{\cos \gamma }}} & {\displaystyle{{\sin \gamma \cos \beta }\over{\cos \gamma }}} \cr 0 & {\cos \beta } & {-\sin \beta } \cr 0 & {\displaystyle{{\sin \beta }\over{\cos \gamma }}} & {\displaystyle{{\cos \beta }\over{\cos \gamma }}}} \right]$$

5.2. State equation

In order to correct the attitude and position, these motion parameters are input into the state vector X:

(22)$${\bi X}=\left[\matrix{ {x_s } \cr {y_s } \cr {z_s } \cr {v_x } \cr {v_y } \cr {v_z } \cr \alpha \cr \beta \cr \gamma } \right]$$

By taking the first derivative of X with respect to time T, the following state equation is deduced

(23)$$\dot{{\bi X}}=f\lpar {\bi X}\rpar $$

in which f(·) is derived from the dynamic Equation (20) and Equation (21).

5.3. Observation equation

According to the result of 1,000 independent simulations of our algorithm without an Extended Kalman Filter (EKF), shown in Figure 6, the errors of the pose estimation of the lander can be defined as white noise whose averages are 0. The lander's parameters are shown in Table 1. The lander's motion parameters estimated in Section 4 can be input into the observation vector Y:

Figure 6. Simulation without the extended Kalman filter for detecting noise: (a) lander's attitude errors; (b) lander's position errors.

Table 1. Representative parameters of the simulation.

(24)$${\bi Y}=\left[\matrix{{x_s } \cr {y_s } \cr {z_s } \cr \alpha \cr \beta \cr \gamma}\right]$$

Next, a linear observation equation is obtained on the basis of the relationship between the state vector X and the observation vector Y

(25)$${\bi Y}={\bi H}\lpar {\bi X}\rpar $$

5.4. Extended Kalman Filter (EKF)

Here, through state Equation (23) and observation Equation (25), an EKF is created.

Firstly, the discrete state model and the discrete observation model are deduced by using Equations (23) and (25):

(26)$${\bi X}_{{\bi k}} =f\lpar {\bi X}_{{\bi k}-1} \rpar +{\bi w}_{{\bi k}-1} $$

(27)$${\bi Y}_{{\bi k}} ={\bi H}\lpar {\bi X}_{{\bi k}} \rpar +{\bi v}_{{\bi k}} $$

where w and v respectively refer to white noise in the system and in the measurements, and

$$\eqalign{{\bi K} & ={\bi E}\lpar {\bi ww}^T\rpar \cr {\bi L} & ={\bi E}\lpar {\bi vv}^T\rpar }$$

in which K and L are defined as the covariance matrix of state and the covariance matrix of observation and E( · ) denotes the expectation.

Then, we can obtain a first order Taylor expansion of the discrete state model:

(28)$${\bi X}_{{\bi k}} ={\bi F}_{{\bi k}-1} {\bi X}_{{\bi k}-1} +f\lpar \hat{{\bi X}}_{{\bi k}-1}\rpar -{\bi F}{{\bi k}-1} \hat{{\bi X}}_{{\bi k}-1} +{\bi w}_{{\bi k}-1} $$

in which Jacobi matrices are obtained by computing the first derivative of f(X) with respect to X as follows:

(29)$${\bi F}=\displaystyle{{\partial f\lpar {\bi X}\rpar }\over{\partial {\bi X}}} $$

Lastly, we can obtain the EKF:

(30)$$\hat{{\bi X}}_{{\bi k}\vert {\bi k}-1} =f\lpar \hat{{\bi X}}_{{\bi k}-1}\rpar $$

(31)$${\bi P}_{{\bi k}\vert {\bi k}-1} = {\bi F}_{{\bi k}-1} {\bi P}_{{\bi k}-1} {\bi F}_{{\bi k}-1}^T +{\bi K}_{{\bi k}-1} $$

(32)$${\bi G}_{{\bi k}} ={\bi P}_{{\bi k}\vert {\bi k}-1} {\bi H}_{{\bi k}}^T \lpar {\bi H}_{{\bi k}} {\bi P}_{{\bi k}\vert {\bi k}-1} {\bi H}_{{\bi k}}^T +{\bi L}_{{\bi k}}\rpar ^{-1} $$

(33)$$\hat{{\bi X}}_{{\bi k}} =\hat{{\bi X}}_{{\bi k}\vert {\bi k}-1} +{\bi G}_{{\bi k}} \lpar {\bi Y}_{{\bi k}} -{\bi H}\lpar \hat{{\bi X}}_{{\bi k}\vert {\bi k}-1}\rpar \rpar $$

(34)$${\bi P}_{{\bi k}} =\lpar {\bi I}-{\bi G}_{{\bi k}} {\bi H}_{{\bi k}}\rpar {\bi P}_{{\bi k}\vert {\bi k}-1} $$

where $\hat{{\bi X}}_{0}$ and P₀ are initial values of the state vector and the covariance matrix of estimation error:

(35)$$\hat{{\bi X}}_0={\bi E}\lsqb {\bi X}_0\rsqb $$

(36)$${\bi P}_0 ={\bi E}\lsqb \lpar {\bi X}_0-\hat{{\bi X}}_0\rpar \lpar {\bi X}_0-\hat{{\bi X}}_0\rpar ^T\rsqb $$

6. SIMULATION RESULTS AND ANALYSIS

In this section, the accuracy and robustness of the proposed algorithm will be tested according to the simulation results of the lander's landing process under different conditions.

In this paper, the measurement noise of landing images is assumed to be white noise. The method for adding the measurement noise to observed feature lines on landing images is shown in Figure 7. The white noise of measurement {(Δs_i, Δe_i)} is added to the feature line endpoints {(s_i, e_i)} on the matched 2D projections in our simulations.

Figure 7. The extracted errors of a feature line's 2D image.

Firstly, the proposed algorithm is simulated without the EKF under different white noise intensity. The lander's attitude and position are estimated when the endpoints of observed feature lines are added to white noise intensity between 1 and 10. The rest of the simulation parameters are shown in Table 1 and the surface altitude of the planet are not considered. In this condition, we simulate the operation of our algorithm 1,000 times for each noise level. The motion errors along with the noise level variation are expressed in Figure 8 after eliminating gross errors by using the 3σ rule.

Figure 8. Simulation without the EKF for testing the impact of noise intensity of the endpoints on accuracy: (a) lander's attitude errors; (b) lander's position errors.

As shown in Figure 8, the lander's attitude errors increase along with the growth of white noise intensity, but they are less than 1·5^° as long as the noise intensity is less than 10 pixels.

The trend of the lander's position error variation is similar to the trend of attitude errors, and the error of the Z-axis which is less than 35 m within the white noise intensity from 1 to 10 pixels is more sensitive than the others which are less than 5 m. In addition, the errors of attitude and position and the corresponding measurement noise are approximately linear, and the errors are not very sensitive with the variation of the noise intensity.

The simulation without the EKF based on the variation of the number of feature lines matched with their 2D projections is discussed to testify to the accuracy and robustness of our algorithm. As the force is very weak when the number of feature lines is equal to three, gross errors in the lander's pose parameters occur easily and the lander's pose parameters can be easily affected by the white noise of feature line endpoints. By using the simulation parameters as shown in Table 1, the average errors of the lander's attitude are approximately 2·5^°, and the average error of the Z-axis is equal to approximately 120 m. Hence, in order to visually display the variation of motion parameters’ error along with the variation of the number of matched feature lines in the simulation figures, it is assumed that the number of matched feature lines increases from four to 14 in the 2D image taken by the on board navigation camera, and the feature line endpoints {(s_i, e_i)} are corrupted with two pixels noise {(Δs_i, Δe_i)}. The rest of parameters of this simulation are kept as per Table 1. The motion errors under the different number of matched feature lines are shown in Figure 9 after eliminating gross errors by using the 3σ rule.

Figure 9. The simulation without the EKF for testing the impact of the number of matched feature lines on accuracy: (a) the lander's attitude errors; (b) the lander's position errors.

As depicted in Figure 9(a), the more the number of matched feature lines, the smaller the lander's attitude errors. When the link number of matched feature lines is more than ten, the errors begin to stabilise within 1^°. At the same time, the error of the X-axis and the error of the Y-axis are less than 1 m, and the error of the Z-axis is less than 10 m in Figure 9(b). The reasons for these are simple. With the increase of the number of matched feature lines, the constraint equation PnL contains more constraint equations, and its force becomes stronger, but when the number of feature lines reaches ten, the increase of its force is minimal. Hence, the errors are close to stable.

It is assumed that the visual navigation method meets the requirements of a precise landing mission when the position errors are less than 20 m and the attitude errors are less than 1·5^° at an altitude of 2,000 m. The ranges of the visual navigation algorithm's applications are detected without the EKF. As the planet surface is assumed to be a plane which includes all feature lines in our algorithm, the relation between the surface altitudes of the planet and the pose errors is explored and simulated. In this simulation, we let the surface altitudes of the planet vary from 0 m to 50 m. The rest of the simulation parameters are kept as per Table 1.

The results of simulation of the surface altitudes of the planet are shown in Figure 10 after eliminating gross errors by using the 3σ rule. The errors of the lander's motion parameters are growing and the error of the Z-axis is easily affected with the increase of the surface altitudes of the planet. When the surface altitudes of the planet are larger than 40 m, the error of the Z-axis is more than 20 m and the errors of β and γ are more than 1·5^°. Hence, in order to meet the requirements of a precise landing mission, the range of our visual navigation algorithm's application for the surface altitudes of the planet is 0 to 40 m at an altitude of 2,000 m.

Figure 10. The simulation without the EKF for testing the impact of the surface altitudes of planet on accuracy: (a) the lander's attitude errors; (b) the lander's position errors.

When the white noise which is added to the endpoints of the feature lines do not change, the different lengths of the feature lines have different influences on our algorithm. Hence, the range of our method's application with regard to the length of the longest feature line is considered. For this, the simulation of our visual navigation algorithm is carried out by using the parameters in Table 1 without the EKF when the length of the longest feature line varies from 100 m to 1,100 m at an altitude of 2,000 m.

In our simulation, the lander's pose errors with regard to the different lengths of the longest feature lines are displayed in Figure 11. When the length of the longest feature line becomes longer, the effects of white noise which are added to the endpoints of the longest feature line are, in theory, reduced. In this simulation, the theory is consistent with the results of the lander's pose estimation. When the length of the longest feature line is more than 200 m, the attitude errors are less than 1^° and the position errors are less than 20 m at an altitude of 200 m. Hence, the range of our visual navigation algorithm's application with regard to the length of the longest feature line is 200 to 1,100 m due to the length of the longest feature line which can be taken at an altitude of 2,000 m being approximately 1,100 m. From Figure 1, the length of the longest feature line which is matched is more than 200 m at an altitude of 2,000 m. Therefore, our algorithm can estimate the lander's attitude and position by using landing images.

Figure 11. The simulation without the extended Kalman filter for testing the impact of the lengths of feature lines on accuracy: (a) the lander's attitude errors; (b) the lander's position errors.

Lastly, the landing process of the NASA spacecraft Curiosity is simulated by using our method with the EKF (Steltzner et al., Reference Steltzner, Miguel San Martin, Rivellini, Chen and Kipp2014). The simulation result can confirm whether the proposed algorithm meets the requirements of a pin-point landing. We hypothesise that the lander's initial height is 8·4 km above the surface after the heat shield separated and the radar-based solution converged. We assume the lander takes 103 seconds to land at an altitude of 247·9 m. The feature line endpoints {(s_i, e_i)} are the same as for the previous simulation. The covariance matrices of the error vectors of this simulation are obtained from the white noise variances of the 1,000 repetitions of the simulation experiments at each height, and the covariance matrices of error vectors are input into the algorithm of Chang et al. (Reference Chang, Xu and Wang2017a; Reference Chang, Xu, Wang and Liu2017b). The covariance matrices of measurement errors are obtained and applied to this simulation of the EKF. Then, the covariance matrices of the state errors are obtained from the system's white noise at each height. The rest of the parameters are kept, and the motion errors along with the landing time are illustrated in Figure 12.

Figure 12. The simulation for landing on Mars: (a) the lander's attitude errors; (b) the lander's position errors.

From Figure 12, the conclusion can be drawn that the result of the simulation by using the EKF is clearly better than the previous solution and the visual navigation algorithm based on crater matching of Shao et al. (Reference Shao, Gao, Xi, Leng and Gu2016). After the extended Kalman filtering, the lander's attitude errors are stable within 0·5^° and the lander's position errors decrease with decreasing altitude to within 1 m at a height of about 247·9 m. This is because 2D images become sharper with the decreasing altitude of the lander. Hence the effect of noise is smaller.

In summary, all simulation of our proposed algorithm proves that the new algorithm can estimate the lander's pose parameters during the descent phase.