Trajectory Recovery and Terrain Reconstruction Based on Descent Images under Dual-Restrained Conditions: Tianwen-1

Abstract

Tianwen-1 is the first Mars probe launched by China and the first mission in the world to successfully complete the three steps of exploration (orbiting, landing, and roving) at the one time. Based on the unverifiable descent images which cover the full range of the landing area, trajectory recovery and fine terrain reconstruction are important parts of the planetary exploration process. In this paper, a novel trajectory recovery and terrain reconstruction (TR-TR) algorithm employing descent images is proposed for the dual-restrained conditions: restraints of the flat terrain resulting in an unstable solution of the descent trajectory and of the parabolic descent trajectory causing low accuracy of terrain reconstruction, respectively. A landing simulation experiment on a landing field with Mars-like landform was carried out to test the robustness and feasibility of the algorithm. The experiment result showed that the horizontal error of the recovered trajectory didn’t exceed 0.397 m, and the elevation error of the reconstructed terrain was no more than 0.462 m. The algorithm successfully recovered the descent trajectory and generated high-resolution terrain products using in-orbit data of Tianwen-1, which provided effective support for the mission planning of the Zhurong rover. The analysis of the results indicated that the descent trajectory has parabolic properties. In addition, the reconstructed terrain contains abundant information and the vertical root mean square error (RMSE) of ground control points is smaller than 1.612 m. Terrain accuracy obtained by in-orbit data is lower than that obtained by field experiment. The work in this paper has made important contributions to the surveying and mapping of Tianwen-1 and has great application value.

1. Introduction

Planetary exploration began in the 1950s, developed to launch various planetary probes in the late 1980s, and entered a new era of planetary exploration with the help of the computer and other technologies by the beginning of the 21st century [1]. Among all the planets, the moon and Mars are the primary targets for exploration [2]. Although China was a late starter in planetary exploration, the series of missions from Chang’e-1 to Chang’e-5 marked the success of the three-step mission of “around, down, and back” in the lunar exploration project [3,4,5]. In addition to the lunar exploration project, Mars exploration and sample return are also part of China’s planetary exploration programs in recent years [6].

Tianwen-1 is China’s first Mars probe to land successfully [7]. After the lander landed on the surface of Mars (Utopia Plain), the Zhurong rover was planned to complete a series of important tasks of patrol and scientific exploration [8]. A large number of images were acquired by sensors carried on the Tianwen-1 orbiter before landing. The Digital Elevation Model (DEM) and Digital Orthophoto Map (DOM) of the pre-selected landing area were produced based on these images [9,10]. However, the resolution of these terrain products is too low to identify small and easily overlooked ground objects (e.g., rocks, fine trenches, uneven slates) which can be obstacles to the rover’s tour [11]. Therefore, the detailed terrain information near the landing site is necessary for the path planning of the rover [12]. Fortunately, descent images have a higher resolution than orbiter images, which can be used to obtain more detailed digital terrain products [13]. Accurate recovery of the landing camera’s trajectory including position and pose is a prerequisite for generating accurate terrain products based on descent images. In addition, it helps us understand and study the influence of wind on the flight status of the lander [14]. Therefore, trajectory recovery and landing area terrain reconstruction using descent images are the important components of the Tianwen-1 mission, which would guarantee the achievement of scientific goals [15]. Moreover, this work is an essential part of near-surface planetary surveying and mapping.

Descent images have played an important role in the probe’s landing and subsequent inspections in planetary exploration missions [16]. The Descent Image Motion Estimation System (DIMES) in the Mars Exploration Rover (MER) mission used descent images to obtain the horizontal velocity of the lander during the entry, descent, and landing (EDL) [17,18]. Xiong et al. proposed to use inertial measurement unit (IMU) and laser altimeter measurements to constrain the descent trajectory solution [19,20]. Random sample consensus (RANSAC) [21] was used to estimate homography matrix (homography matrix is a mapping between two planes [22]) and removes mismatching points (matching points not corresponding to the same spatial point, where matching points means projections of the same spatial point on different images) to recover the Chang’e-3 lander’s descent trajectory and map the landing area [23]. Meng et al. improved Xiong et al.’s method [20] by utilizing the homography matrix in reconstructing the terrain of the landing area [24,25]. Di et al. reconstructed the terrain of the landing area using descent images of Chang’e-3, but no detailed algorithm is mentioned in the paper [15]. Xu et al. proposed an optimization algorithm to improve the surface normal vector by combining the reflectivity model [26]. For the Chang’e-4 mission, Liu et al. made use of the collinear equation combined with least squares to recover the lander trajectory [27], while Di et al. mapped the landing area [28,29,30,31].

Descent trajectory recovery requires feature points of the descent images. Regarding feature point extraction and matching, operators with scale, rotation, and affine invariance include Scale-invariant feature transform (SIFT) [32], Speeded Up Robust Features (SURF) [33] and Binary Robust Invariant Scalable Keypoints (BRISK) [34], among which SIFT is found to be the most accurate descriptor for feature detector [35]. Brute-Force (BF) Matcher and Fast Library for Approximate Nearest Neighbors (FLANN) based Matcher [36] are often applied for feature point matching and combined with k-Nearest-Neighbor (KNN) to filter bad matches [37]. RANSAC combined with a minimal pose solver is the most commonly used method to eliminate mismatching points [21]. However, these methods are not effective in eliminating mismatching points in the descent images of Tianwen-1.

Geometry model estimation based on matching points is an important step of trajectory recovery. The scene of the Tianwen-1 descent images is predominantly but not entirely planar, which is (quasi-) degenerate data for robust geometry model estimation. Rebert et al. proposed a parallax beam algorithm to process the degenerate planar scene [38]. RANSAC for the (Quasi-) Degenerate data (QDEGSAC) algorithm solves the degenerate plane problem by running the RANSAC algorithm multiple times to calculate the rank of the data matrix [39]. Decker et al. detect degeneracy by continuously testing motion models with decreasing degrees of freedom [40]. The Degenerate RANSAC (DEGENSAC) perceives the degeneracies within the RANSAC sample test [41]. Torr proposed the Geometric Robust Information Criterion (GRIC) model selection approach which selects the model with the lowest score [42]. However, the GRIC model selection approach would not obtain a right model for a dominant plane scene.

In the terrain reconstruction step, the plane-sweep method is proposed to be applied to generate dense point clouds for the specific baseline direction of descent images [43]. Xiong et al. realized depth map recovery with the plane-sweep method for unstructured scenes in descent images [19]. Meng et al. modified the plane-sweep method by using the Zero-normalized cross-correlation score (ZNCC) method instead of the Sum of Squared Differences (SSD) to find the image correlation after warping, followed by combining best seed first propagation (BSFP) strategy with neighborhood value propagation to accomplish depth recovery [24,25].

The above methods have several critical disadvantages when applied to descent images of Tianwen-1 due to the dual restraints: the first restraint of terrain conditions on the trajectory recovery and the second restraint of descent trajectory on the terrain reconstruction (Figure 1). For the first restraint, the previous methods (1) cannot extract enough feature points in the textureless region (flat region with no obvious features) of descent images, which reduces the accuracy of the recovered trajectory; (2) cannot effectively remove the mismatching points caused by a large number of similarly shaped dunes in the landing area which results in the failure of trajectory recovery; (3) sometimes give the wrong trajectory solution due to (quasi-) degenerated data caused by the flat terrain. For the second restraint, the previous methods result in low accuracy of terrain reconstruction based on descent trajectory with the baseline (line connecting cameras’ optical center) approximately parallel to the optical axis of landing camera. Therefore, the dual-restraints problem brought by the combination of terrain conditions and special descent trajectory creates difficulties for the TR-TR process.

In this paper, a novel trajectory recovery and terrain reconstruction (TR-TR) algorithm is proposed for the dual constraints to overcome the above drawbacks when descent trajectory recovery and terrain reconstruction interact with each other. Firstly, the TR-TR algorithm obtains adequate feature points in the textureless region of descent images by reducing the SIFT contrast threshold, which improves the accuracy of the recovered trajectory. Next, the algorithm eliminates mismatching points by using scale constraints and a priori terrain information, which improves the success rate of trajectory recovery. Then, the algorithm uses the epipoles (the intersection of baseline and image) constrained parallax beam method to obtain a robust trajectory from (quasi-) degenerated data caused by the flat terrain. Finally, the algorithm allows the virtual planes to have multiple directions in the plane-sweep method, which improves the accuracy of terrain reconstruction under the constraint of the special trajectory. The TR-TR algorithm was tested in field experiment with improvements made based on the results, and the accuracy was evaluated. The algorithm was subsequently applied to the descent images of Tianwen-1′s in-orbit data to recover the descent trajectory and reconstruct the terrain of Tianwen-1′s landing area, which provided effective product support for the mission.

This paper presents the work of obtaining the position and pose information of the landing camera and reconstructing the terrain of the landing area. Section 1 gives the background of the research. Section 2 introduces the Tianwen-1′s in-orbit data and the TR-TR algorithm proposed. Section 3 shows the validation of the novel algorithm through field experiment. Section 4 demonstrates Tianwen-1′s in-orbit data processing results. Section 5 provides a discussion, summarizes current work, and looks forward to future work.

2. Data and Methodology

2.1. Data

The descent images used were continuously captured during the EDL phase by the landing camera, which is a downward-facing optical camera (parameters are shown in

Table 1. Landing camera parameters.

Name	Value
Diagonal field of view (FOV)	43°
Square FOV	30° × 30°
Focal length	20 mm
Spectral range	500~800 nm
Entrance pupil diameter	≥4 mm
Entrance pupil position	21.9 mm
Integral time	0.03~64 ms
Image area resolution	2048 × 2048
Pixel size	5.5 μm

) mounted on the bottom of the landing platform.

During the whole descent process of Tianwen-1, the landing camera acquired a total of 84 descent images, among which the preceding images numbered 1~22 were sampled too early to be valuable, whereas the last images numbered 43~84 were acquired too late to be applied in the mission. Therefore, the 20 descent images numbered from 23 to 42 at a suitable altitude were chosen to recover the position and pose of the landing camera during the EDL phase as well as to map the landing area. The acquisition time of the next 19 images is as follows in

Table 2. Acquisition times of the descent images.

Image ID	Acquisition Time (s)	Image ID	Acquisition Time (s)
23	T23 + 0.000	33	T23 + 65.536
24	T23 + 28.672	34	T23 + 81.273
25	T23 + 32.768	35	T23 + 94.209
26	T23 + 36.864	36	T23 + 98.305
27	T23 + 40.960	37	T23 + 102.401
28	T23 + 45.056	38	T23 + 106.497
29	T23 + 49.152	39	T23 + 110.593
30	T23 + 53.248	40	T23 + 114.689
31	T23 + 57.344	41	T23 + 118.785
32	T23 + 61.440	42	T23 + 122.881

if the acquisition time of the first image (ID of 23) is recorded as T23.

2.2. Methodology

The four main steps of the TR-TR algorithm include feature extraction and matching, mismatching points removal, robust motion estimation, and terrain reconstruction. Details are as follows:

(I)

For the feature extraction and matching part, SIFT with a low contrast threshold is applied to extract and descript feature points, after which BF matcher with KNN (k = 2) is utilized to filter bad matches.

(II)

The TR-TR algorithm uses two constraints to remove mismatching points: the scale monotonicity constraint and the a priori terrain information constraint. The scale monotonicity constraint takes advantage of the fact that the feature point scale monotonically varies with the descent image order to remove mismatching points. The constraint of prior terrain information requires first finding the homography matrix (H) corresponding to the points in the scene plane by a contrario RANSAC (AC-RANSAC) method. Then, DEM of the landing area generated from the orbiter images is exploited as prior terrain information to find the upper boundary of the plane induced parallax (introduced in Section 2.2.2) corresponding to the homography matrix (H), and, finally, remove points beyond the boundary.

(III)

According to the two solutions obtained by the homography matrix (H) decomposition, two possible positions of epipoles are obtained as constraints to perform the epipoles constrained parallax beam (ECPB) algorithm. If the fundamental matrix is obtained, the determined motion solution can be obtained directly by decomposition of the fundamental. If the fundamental matrix cannot be obtained, the solution that does not conform to the parabolic descent trajectory constraint in the two solutions obtained by decomposition of the homography matrix (H) is eliminated to obtain the determined motion solution.

(IV)

After refining the motion by bundle adjustment, the accurate descent trajectory is obtained. A lot of sparse spatial points are also obtained when estimating motion between images in the previous step. Triangular mesh is established based on these sparse spatial points. Each triangle in the mesh is projected to all images that can see the triangle. Based on these projections, multi-image normalized cross-correlation score (NCC) curves are drawn and initial seed points are determined accordingly. Initial seed points are then propagated to obtain depth maps for generating dense point clouds. Finally, DEM and DOM of the landing area are generated based on dense point clouds.

The flow chart of the TR-TR algorithm proposed in this paper for processing the descent images is given below (Figure 2).

2.2.1. Feature Extraction and Matching

Compared with other operators, SIFT describes the position of points most accurately, which can improve the trajectory recovery accuracy. Therefore, SIFT is selected to extract and describe feature points in descent images. However, the uneven distribution of ground features in the landing area of Tianwen-1 leads to the emergence of large textureless areas in the descent images. Unfortunately, in these textureless areas, SIFT with default threshold in OpenCV (a famous open source computer vision library) is not able to extract sufficient feature points. In order to solve this problem, the contrast threshold of the SIFT operator should be set as small as possible.

The principle of feature point matching is to compare the Euclidean distance of two feature point descriptors as a similarity measure. For feature point $p_1$ in image $I_1$, to find its matching point $p_2$ in image $I_2$, the BF matcher combined with KNN ($\mathrm{k}=2$) is used to return the two best matches $p_{21}$ and $p_{22}$ (with similarity measure $s_1$ and $s_2$, respectively) from image $I_2$. If $s_1/s_2$ is smaller than a similarity ratio threshold (assuming $s_1<s_2$), $p_{21}$ is considered to be the correct matching point. The BF matcher is better than the FLANN matcher because although the FLANN matcher is faster, the processing time of the descent images is sufficient. Since the landing camera moves forward along the descent trajectory, there will be overlapping areas in successive descent images. Hence, it is not enough to match feature points only in adjacent images, but it should be done in multiple consecutive images to obtain feature tracks.

2.2.2. Mismatching Points Elimination

Meng et al. uses RANSAC to search the fundamental matrix among all matching points, and removes the matching points that do not conform to the calculated fundamental matrix as mismatching points [24]. This method is suitable for a general scene, but not suitable for the scene in the Tianwen-1 descent images. The scene in the Tianwen-1 descent images is a planet surface with flat terrain, resulting in (quasi-) degenerated data. (Quasi-) degenerated data may not provide sufficient constraints on the fundamental matrix due to rank deficiency; hence, an incorrect fundamental matrix may be obtained. RANSAC consistent with an incorrect fundamental matrix will remove a large number of correct matching points. Besides, Meng’s method does not take advantage of the characteristics of the descent images, such as the monotonicity of scale changes. Based on this problem, the TR-TR algorithm chooses to use the homography matrix to eliminate mismatching points instead of using the fundamental matrix. Since most of the points (or all points) in the (quasi-) degenerated data are located in a scene plane, these data may not form sufficient constraints on the fundamental matrix, but it must form sufficient constraints on the homography matrix. Therefore, a robust homography matrix can be obtained. In addition to robustness, the homography matrix has another advantage that it can make use of the prior terrain information, which the fundamental matrix cannot do. Apart from the consistency method, the TR-TR algorithm also makes use of the scale monotonicity of the descent images.

In the TR-TR algorithm, the scale monotonicity constraint is utilized first to eliminate part of the mismatching points. The landing camera’s forward motion causes the scale (obtained from the SIFT descriptor) of one feature track in the descent images to increase monotonically with the acquisition time. As a result, feature tracks not meeting this condition must contain mismatching points. For such feature tracks, the contained features whose scale do not vary monotonically are removed to maintain the scale monotonicity of the feature tracks.

Next, mismatching points are eliminated by prior terrain information constraint. Data are degenerate when all points lie inside a plane, whereas those data are quasi-degenerate when most points lie inside a plane with the rest lying outside the plane. The terrain in the landing area of Tianwen-1 is flat and approximates a scene plane for the landing camera. Most of the ground points are inside the scene plane, with occasional undulating ground points (e.g., volcanoes, gullies, impact craters) located outside the scene plane. The maximum elevation difference of the landing area can be evaluated from the DEM generated by the orbiter images of Tianwen-1. Restricted by the maximum elevation difference, when the approximate acquisition height corresponding to each descent image is known (supplied by laser altimeter), the plane induced parallax (introduced in Figure 3 and denoted as $\rho$) of points outside the plane must have an upper boundary (denoted as ${\rho }_{max}$). The points beyond the boundary must be mismatching points to remove.

In order to obtain the plane induced parallax of corresponding points, the homography matrix needs to be solved first. The homography matrix computed in conjunction with the AC-RANSAC framework is robust to noise and mismatches [44]. The number of false alarms ($NFA$) is determined by AC-RANSAC as:

$$NFA\left(H,k\right)=\left(n-4\right)\left(\begin{array}{c}n\\ k\end{array}\right)\left(\begin{array}{c}k\\ 4\end{array}\right){\left({\epsilon }_k^2{\alpha }_0\right)}^{k-4}$$

(1)

where $H$ is the homography matrix to be tested, ${\epsilon }_k$ denotes the $k$th minimum residual, ${\alpha }_0$ indicates the proportion of the area of the circle with a radius of one pixel in the image, $n$ represents the number of total matching points.

The homography matrix ($H$) is considered valid if

$$NFA\left(H\right)=\underset{5\le k\le n}{min}NFA\left(H,k\right)\le 1$$

(2)

Let $\widehat{H}$ indicate the homography matrix obtained by AC-RANSAC and $\epsilon$ the corresponding maximum residuals; the correct matching points set $P$ can then be figured out according to $\widehat{H}$ and $\epsilon$. Introducing set $U$ as the total matching points, the rest of the matching point set is obtained as the complementary set of $P$ from $U$, which is represented by ${\complement }_{U}P$. ${\complement }_{U}P=\varnothing$ implies that all matching points lie inside the plane corresponding to $\widehat{H}$, which means that the data are degenerate. If not, matching points with too large plane induced parallaxes are considered to be mismatching points and are removed from ${\complement }_{U}P$, leaving the matching points set ${\complement }_U{P}_{\rho <{\rho }_{max}}$.

2.2.3. Robust Motion Estimation

Meng et al. uses RANSAC to search the fundamental matrix among all matching points, and decompose the fundamental matrix to obtain the motion [24]. In this way, there is an obvious disadvantage that the (quasi-) degenerated data caused by flat terrain will easily result in an incorrect fundamental matrix, and then get the wrong motion. By contrast, in the TR-TR algorithm, motion can be obtained robustly. As mentioned in the previous section, a robust homography matrix is obtained from (quasi-) degenerated data. Based on the obtained homography matrix, the corresponding fundamental matrix is estimated according to the relationship between the fundamental matrix and the homography matrix if the data provide sufficient constraints. Specifically, the robust fundamental matrix is estimated from the points (may including some mismatching points) located outside the plane corresponding to the obtained homography matrix using the epipoles constrained parallax beam method [38]. The robust motion is obtained by decomposing the fundamental matrix. If the data cannot provide sufficient constraints, the fundamental matrix will not be obtained. For such cases, the homography matrix is decomposed into two motion solutions. The TR-TR algorithm utilizes the parabolic property of the descent trajectory to eliminate the one that does not match the actual situation, and the remaining motion solution is robust.

Continuing from the previous section, ${\complement }_U{P}_{\rho <{\rho }_{max}}=\varnothing$ means that matching points, except those in the plane corresponding to $\widehat{H}$, are all mismatching points, indicating a case of degenerate data. Conversely, if ${\complement }_U{P}_{\rho <{\rho }_{max}}\ne \varnothing$, there may exist correct matching points that lie outside the plane which provides constraints on the epipolar geometry to estimate the fundamental matrix.

$\widehat{H}$ is decomposed to obtain two possible epipoles $e_1^{\prime }$ and $e_2^{\prime }$ on the second image. The $i$th matching point is imaged at $x_i$ and $x_i^{\prime }$ for the first and second view, respectively. Ideally, the line between $x_i^{\prime }$ and $\widehat{H}{x}_i$ would either pass through point $e_1^{\prime }$ or $e_2^{\prime }$, corresponding to the two possible fundamental matrices $\widehat{F_1}$ or $\widehat{F_2}$, respectively. However, in practice, it will not happen because both $x_i^{\prime }$ and $\widehat{H}{x}_i$ are corrupted with noise. For simplicity, let $x_i^{\prime }$ and $\widehat{H}{x}_i$ have the same point error that has a distribution area of the circle of radius $\epsilon$ with the point as the center since $\epsilon$ (figured out from AC-RANSAC) represents point error more accurately. The parallax beams [38] of $x_i^{\prime }$ and $\widehat{H}{x}_i$ are drawn to judge if they cover $e_1^{\prime }$ or $e_2^{\prime }$. Set the initial values of two counters $N_{\widehat{F_1}}$ and $N_{\widehat{F_2}}$ to be zero. Let $N_{\widehat{F_1}}=N_{\widehat{F_1}}+1$ if $e_1^{\prime }$ is covered. Similarly, let $N_{\widehat{F_2}}=N_{\widehat{F_2}}+1$ if $e_2^{\prime }$ is covered. The geometric model ($M$) of the image pair is determined as:

$$\{\begin{array}{c}M=\widehat{H}, if N_{\widehat{F_1}}=N_{\widehat{F_2}}\\ M=\widehat{F_1}, if N_{\widehat{F_1}}>N_{\widehat{F_2}} \\ M=\widehat{F_2}, if N_{\widehat{F_1}}<N_{\widehat{F_2}}\end{array}$$

A definite solution of motion comes from the decomposition of the fundamental matrix if $M=\widehat{F_1}$ or $M=\widehat{F_2}$ when knowing the camera matrix (K), whereas two possible solutions come from the decomposition of the homography matrix if $M=\widehat{H}$. In fact, the parabolic property of the descent trajectory can be used to remove the unrealistic solution.

For convenience, x and y are defined [45]:

$$x=\frac{R^{T}t}{\left|\right|t\left|\right|}$$

(3)

$$y=\left|\right|t\left|\right|n$$

(4)

where $R$ and $t$ are the rotation matrix and translation vector normalized with respect to the depth of the plane of which the normal vector is represented by $n$.

Equations (5) and (6) show the relationship of two solutions obtained from the decomposition of the homography matrix, which is distinguished by subscripts a and b:

$$y_b=\pm \frac{\left|\right|y_a\left|\right|}{\rho }\left(y_a+2x_a\right)$$

(5)

$$x_b=\frac{1}{\rho }\left(\frac{\nu }{\left|\right|y_a\left|\right|}{y}_a-||y_a||x_a\right)$$

(6)

where,

$$\rho =\sqrt{\left|\right|x_e\times y_e\left|\right|{}^2+{\left(x_e^T{y}_e+2\right)}^2}$$

(7)

$$\nu =2\left(x_e^T{y}_e+1\right)$$

(8)

where,

$$e=\left\{a,b\right\}$$

(9)

Replacing (3) and (4) into (5) and (6), we can figure out that

$$n_b=\pm \frac{\left|\right|t_a\left|\right|\left|\right|n_a\left|\right|}{\rho \left|\right|t_b\left|\right|}\left(||t_a||n_a+\frac{2R_a^T{t}_a}{\left|\right|t_a\left|\right|}\right)$$

(10)

$$R_b^T{t}_b=\frac{\left|\right|t_b\left|\right|}{\rho }\left(\frac{\nu n_a}{\left|\right|n_a\left|\right|}-||n_a||R_a^T{t}_a\right)$$

(11)

For one descent image pair of Tianwen-1, denoting $O_L$ as the camera center of the image taken at a lower altitude and $O_H$ as the camera center of the image taken at a higher altitude, there are two possible locations of $O_H$ (marked as $O_{Ha}$ and $O_{Hb}$) corresponding to two solutions if the camera $O_L$ is fixed with a view along the inverse direction of the y-axis of the camera coordinate system, as shown in Figure 4. It can be concluded from Figure 4 that the angle ${\theta }_a$ between $\left|\right|t_a\left|\right|n_a$ and $\overrightarrow{O_L{O}_{Ha}}$ satisfy Equation (12) when assuming that $O_{Ha}$ is in the quadrant $\left\{+x,-z\right\}$ ($O_{Ha}$ must be located in the quadrant $\left\{\pm x,-z\right\}$, since $O_{Ha}$ is higher than $O_L$):

$$\pi /2<{\theta }_a<\pi$$

(12)

$R_a^T{t}_a$ (equaling the vector $\overrightarrow{O_{Ha}{O}_L}$) has the same direction as $\frac{2R_a^T{t}_a}{\left|\right|t_a\left|\right|}$, which implies that $\frac{2R_a^T{t}_a}{\left|\right|t_a\left|\right|}$ is a vector pointing to the quadrant $\left\{-x,+z\right\}$. Equation (8) shows there are two possible solutions (in opposite directions of $\left\{-x,+z\right\}$ and $\left\{+x,-z\right\}$ according to the vector algorithm) of $n_b$, among which apparently only the one solution pointing $\left\{-x,+z\right\}$ fits the reality.

$\frac{\nu n_a}{\left|\right|n_a\left|\right|}$ ($\mathsf{\nu }>0$) has the same direction with $\left|\right|t_a\left|\right|n_a$, which points towards the plane perpendicularly. At the same time, $\left|\right|n_a\left|\right|R_a^T{t}_a$ has the same direction with $\frac{2R_a^T{t}_a}{\left|\right|t_a\left|\right|}$. It can be inferred from Equation (9) that $R_b^T{t}_b$ points to the quadrant $\left\{+x,\pm z\right\}$ since $\mathsf{\rho }>1$. Obviously, $\overrightarrow{O_{Hb}{O}_L}$ cannot point toward $\left\{+x,-z\right\}$ because $O_{Hb}$ is higher than $O_L$, which indicates that $R_b^T{t}_b$ actually points to the quadrant $\left\{+x,+z\right\}$, since $R_b^T{t}_b=\overrightarrow{O_{Hb}{O}_L}$. Therefore, the angle ${\theta }_b$ between $R_b^T{t}_b$ and $n_b$ satisfies the following equation:

$$0<{\theta }_b<\pi$$

(13)

where $0<{\theta }_b<\pi /2$ violates the visibility constraint of $I_H$, so:

$$\pi /2<{\theta }_b<\pi$$

(14)

In Figure 4, two surfaces are drawn for Solution a and Solution b separately, which does not correspond to the actual situation in which there is only one surface that can be used to analyze the relationship between two possible locations of $O_H$ when fixing $O_L$. Next, the surfaces of Solution a and Solution b are brought together to get the following schematic diagram:

According to Equations (12) and (14), as shown in Figure 5, it can be concluded that $O_{Ha}$ and $O_{Hb}$ must not lie on the same side of the surface normal line of the Mars surface at the point $O_L$.

The above discussion is for the case when $O_{Ha}$ lies in the $\left\{+x,-z\right\}$ quadrant. It still satisfies the above conclusion when $O_{Ha}$ lies in the $\left\{-x,-z\right\}$ quadrant due to the property of equivalent exchange of the two solutions [45].

For simplicity, a 2D case has been given above. It is not difficult to deduce that two solutions follow the same principle in the actual 3D case, except that the normal line is replaced by the normal plane which is perpendicular to the horizontal component of the velocity of the lander and passes through $O_L$. In practice, the solution in which $O_H$ lies on the opposite side of the direction of the horizontal component of the normal plane’s velocity is the correct solution.

2.2.4. Terrain Reconstruction

Due to the forward motion of the descent images, the epipole is located inside the image. The window correlation method used in binocular vision is not applicable in descent images. If we resample one image along the epipolar line to make their resolution equivalent as in stereo, then the images will be oversampled near the epipole while undersampled at the image boundaries. Hence, Meng proposed using the plane-sweep method to reconstruct the terrain [24]. The specific method is to use several parallel virtual planes to cut the terrain. For one point on image 1, the projection window of its correlating window on each virtual plane on image 2 can be obtained by homography projection. Consequently, several candidate windows are obtained on image 2 about the correlation window of that point on image 1. The ZNCC of the correlation window on image 1 and the candidate windows on image 2 are calculated as the correlation measure, and the depth value of the virtual plane corresponding to the candidate window with the greatest correlation is taken as the depth value of that point on image 1. The depth map of image 1 is generated by this method. Subsequently, Meng studied the shape of ZNCC curves and formulated rules to filter out ZNCC curves with poor conditions (such as too flat or with two similar peaks). A point whose ZNCC curve is qualified is selected as the initial seed point. Then, the initial seed points are propagated using BSFP strategy to generate a depth map [25].

Meng’s method has several disadvantages. Firstly, the method only discussed the depth recovery of two images, but did not utilize the other images. Secondly, the virtual planes used by this method in the plane-sweep algorithm has only one direction; thus, this method has low depth accuracy for the recovery of forward-tilting terrain, which is not good at the recovery of detailed terrain. Thirdly, this method does not use the sparse spatial points obtained in the previous motion estimation steps. These sparse spatial points are very accurate prior information which can be used to improve the terrain accuracy. In contrast, the TR-TR algorithm proposes the plane-sweep method using multiple images. The TR-TR algorithm uses virtual planes that do not have only one direction, so that they can better fit the small-scale undulating terrain of the Tianwen-1 landing area. TR-TR algorithm takes advantage of sparse spatial points to increase the accuracy of reconstructed terrain. Specific methods are as follows.

After camera motion refining and bundle adjustment, a series of sparse spatial points are obtained, which are employed to constitute a mesh containing many spatial triangles. The image $I_L$ with the lowest acquisition height is moved out from the image set obtained by searching images containing some triangle $\mathrm{P}$, leaving an image set denoted as $\mathsf{\Gamma }=\left\{I_i|\mathrm{i}=1\dots \mathrm{N}\right\}$. $\mathrm{M}$ virtual planes parallel to $\mathrm{P}$ are created around $\mathrm{P}$ and are marked as $\left\{V_k | k=1\dots \mathrm{M}\right\}$. The homography ${\mathrm{H}}_k^i$ induced by $V_k$ between $I_L$ and $I_i$ in $\mathsf{\Gamma }$ is:

$${\mathrm{H}}_k^i=K\left(R_i{R}_L^{-1}+\frac{R_i\left(C_L-C_i\right)n_k^T}{n_k^T{X}_k}\right)K^{-1}$$

(15)

where, K is the camera calibration matrix, $X_k$ represents a point passed through by $V_k$, $n_k^T$ is defined as the normal of $V_k$, $R_i$ is the rotation matrix of $I_i$, $R_L$ is the rotation matrix of $I_L$, $C_i$ is the location of $I_i$ and $C_L$ is the location of $I_L$.

Marking $P_L$ and $P_i$ as projections of P on $I_L$ and $I_i\left(I_i\in \mathsf{\Gamma }\right)$, respectively, the NCC of pixel $\left(u,v\right)\in P_L$ between $I_L$ and $I_i\left(I_i\in \mathsf{\Gamma }\right)$ is given:

$$C_k^i\left(u,v\right)={\displaystyle \sum }_{\left(u,v\right)\in W_c}w\left(u,v\right)·\tau \left(I_L^{ik}\left(u,v\right),I_i\left(u,v\right)\right)$$

(16)

where $W_c$ is the correlation window, $w\left(u,v\right)$ indicates the weight of pixel $\left(u,v\right)$ and $\tau \left(I_L^{ik}\left(u,v\right), I_i\left(u,v\right)\right)$ computes the dissimilarity of $I_L^{ik}$ and $I_i$ with regard to $\left(u,v\right)$. $w\left(u,v\right)$ is obtained by:

$$w\left(u,v\right)=\{\begin{array}{c}1,\left(u,v\right)\in P_i\\ e^{-d},\left(u,v\right)\notin P_i\end{array}$$

(17)

where $d$ is the distance between $\left(u,v\right)$ and $P_i$, implying less weight corresponding to further distance.

$\tau \left(I_L^{ik}\left(u,v\right),I_i\left(u,v\right)\right)$ is defined as:

$$\tau \left(I_L^{ik}\left(u,v\right),I_i\left(u,v\right)\right)=\gamma \left|\right|I_L^{ik}\left(u,v\right)-I_i\left(u,v\right)\left|\right|+\left(1-\gamma \right)\left|\right|\nabla I_L^{ik}\left(u,v\right)-\nabla I_i\left(u,v\right)\left|\right|$$

(18)

where $\left|\right|I_L^{ik}\left(u,v\right)-I_i\left(u,v\right)\left|\right|$ computes the L1-distance of $I_L^{ik}\left(u,v\right)$ and $I_i\left(u,v\right)$’s grey value space, $\left|\right|\nabla I_L^{ik}\left(u,v\right)-\nabla I_i\left(u,v\right)\left|\right|$ computes the absolute difference of the gray-value gradients of $I_L^{ik}\left(u,v\right)$ and $I_i\left(u,v\right)$, and $\gamma$ is a user-defined parameter.

The final score $C_k\left(u,v\right)$ of $\left(u,v\right)$ with respect to $V_k$ is obtained by summing up the $C_k^i\left(u,v\right)$ calculated from all images in $\mathsf{\Gamma }$:

$$\tau C_k\left(u,v\right)={{\displaystyle \sum }}_{i=1}^N{C}_k^i\left(u,v\right)$$

(19)

The points at which the $C_k\left(u,v\right)$ curve fluctuated with respect to $k$ fits the condition are selected as seed points, as well as sparse spatial points. The seed points are propagated within the triangle $P_L$ to generate the depth map of $I_L$ of triangle $P$, thus obtaining the point clouds of $P$. In the same way, the point clouds of other triangles are calculated and fused to obtain the final point clouds of the whole terrain. The terrain point clouds are projected to the horizontal plane in the vertical direction according to the Mars geographic coordinate system to generate a dense Delaunay triangle network. For each triangle, the elevation and grayscale values of the raster points within the triangular surface slice can be obtained by interpolation based on the elevation and grayscale values of the three vertices. Using this method, the elevation and grayscale values of all raster points in the landing area can be obtained to generate the DEM and DOM of the landing area.

3. Field Experiment

3.1. Experiment Overview

In order to verify the method proposed in this paper, a corresponding field experiment in the Red Cliff area with a Mars-like terrain was designed. The experimental data were captured in November 2019 by the unmanned aerial vehicles (UAVs), DJI Phantom 3 and Phantom 4, during descent when simulating the lander of Tianwen-1, through which far more than 2000 images were acquired. The Red Cliff area is located in the Dachaidan region of Qinghai Province, China. It is the first Mars simulation basement in China and has topographic of great value for studying planetary landforms [46]. Before aerial photography, a number of checkerboard grid plastic boards were deployed as ground control points (Figure 6). The coordinates of the ground control points were measured using the Topcon total station. The true value of the position and pose of the camera corresponding to each descent image were evaluated from the coordinates of the ground control points. In addition to the simulated landing, the UAV flew along multiple airstrips and acquired a large amount of image data which were processed to obtain the true values of the DEM of the experimental area (Figure 7). It can be seen that there are three similarities between the landforms of the experimental area and the landing area of Tianwen-1 (refer to the Tianwen-1′s orbiter image): (1) the terrain is relatively flat (Utopia Plain of Mars); (2) both contain striated feature (sand ridges on Mars); and (3) both contain small areas of obvious topographic relief (impact craters on Mars).

3.2. Matching Accuracy and Initial Motion

The combination of the mismatching points elimination methods proposed in this paper achieves higher matching accuracy. In addition, the proposed initial motion estimation algorithm in this paper achieves more robust and accurate motion results. In order to verify this, the result of Meng’s method [24] (searching fundamental matrix directly by RANSAC) is compared in the experiment (Figure 8). The comparison of matching accuracy and the comparison of initial motion error are put together in this section because the matching result directly affects the accuracy and robustness of initial motion estimation.

3.2.1. Matching Accuracy

A combination of scale constraints and homography matrix constraints are proposed to eliminate mismatching points in this paper. In the experiment, the true values of the corresponding fundamental matrices can be evaluated based on the true values of the translations and poses of the descent images. It is assumed that the fundamental matrix of an image pair $\left\{I_1,I_2\right\}$ is denoted as $F$. A matching point is marked as $x$ on image $I_1$ and $x\prime$ on image $I_2$. The epipolar line through $x\prime$ is represented by $l$. Taking the distance from $x\prime$ to $l$ as a measure, a matching point with such distance less than one pixel is considered to be a correct matching point. The matching accuracy is obtained by dividing the number of correct matching points by the total number of matching points. Twelve out of a total of forty-five pairs were sampled to compare matching accuracy. It is obvious that the matching accuracy of four image pairs obtained by Meng’s method is too low. This is due to the elimination of correct matching points since the incorrect fundamental matrices are evaluated from (quasi-)degenerate data. These phenomena appear randomly among image pairs, whereas the geometric models obtained by the method of this paper avoid this randomness and thus can robustly eliminate the mismatching points. Apart from this, the matching accuracy of the Meng’s method is still not as high as that obtained by the method of this paper. The reason for this is that Meng’s method does not consider using the scale constraint of descent images to eliminate the mismatching points.

3.2.2. Initial Motion

The ECPB algorithm was used to estimate an accurate geometry model of an image pair, followed by the utilization of parabolic trajectory to remove the ambiguity of solutions derived from the decomposition of the homography matrix. In order to verify that the proposed algorithm can obtain more robust motion results, the relative pose and translation of twelve image pairs selected from total forty-five image pairs by equal interval are compared in the experiment with Meng’s method [24]. Compared to the vertical translation, the horizontal translation has more influence on the accuracy of terrain reconstruction. Thus, only the horizontal translation error (obtained by differentiation from the true value) is compared. The rotational angle error calculated by Meng’s method is sometimes too large (gross error). In order to reflect the details of the comparison, the maximum value of the angle error is set to 0.5°, for which reason all angle errors over 0.5° are displayed as 0.5° (Figure 9).

From Figure 9, it can be seen that almost half of the motion solutions obtained through Meng’s method [24] have gross errors in angle, which means they are completely wrong solutions although they do not have gross errors in horizontal motion error. These wrong motion solutions are considered to be results of decomposition of incorrect fundamental matrices obtained from the (quasi-)degenerate data. In fact, this result corresponds to the matching accuracy result in Figure 9 (the incorrect fundamental matrices lead to very low matching accuracy and wrong pose solutions).The results of this paper’s algorithm do not have gross errors, which provides stable initial values for the subsequent steps of bundle adjustment.

It can also be seen from Figure 9 that the accuracy of other correct pose solutions obtained by Meng’s method is still not as high as that obtained by this paper’s method, which is mainly because this paper’s method has higher matching accuracy in other image pairs than the method of Meng.

3.3. Feature Points Extraction and Refined Motion

The feature points extraction method proposed in this paper obtains denser and more uniform feature points, which achieves more accurate refined motion results after bundle adjustment. In order to verify it, the feature points extracted by SIFT with default threshold in OpenCV is compared in the experiment.

3.3.1. Feature Points Extraction

In OpenCV, the default contrast threshold for SIFT is 0.04. In the method of this paper, this value is set to 0.01; thus, more matching points in the textureless area of the descent images can be obtained; but at the same time, there will be more points in other texture-rich areas, resulting in an increase in computation. Hence, the method of this paper also sets the maximum number of extracted feature points by region to limit the number of feature points that are too dense in some areas of descent images.

In order to visually compare the extraction of feature points, a descent image containing a textureless area is selected to show the results of the two algorithms (Figure 10).

It is shown in Figure 10 that evenly distributed feature points are obtained by the method of this paper. These feature points are matched across all descent images according to the method described in Section 2.2.1 to obtain a large number of feature tracks. The comparison of the number of feature tracks obtained by the SIFT and the method of this paper is displayed in

Table 3. Number of feature tracks obtained by SIFT and the method of this paper.

Length of Feature Tracks	SIFT	The Method of This Paper
2	1105	2176
3	772	3680
4	383	4262
5	34	2997
6	7	1821
7	0	782
8	0	215
9	0	16

. Feature tracks of different lengths (number of images spanned by a feature track) are shown in the table.

3.3.2. Refined Motion

Employing reliable initial values, the uniform and dense feature points extracted by this paper’s algorithm and the SIFT feature points extracted by OpenCV with default threshold are used as the observed data of the bundle adjustment. After making a difference with true value, camera pose and translation errors are displayed for comparison (Figure 11).

It can be seen that uniformly and densely extracted and matched feature points help to obtain a more accurate motion. This is the result of control points involved adjustment. The error is small at the beginning of the descent, followed by a slow accumulation with descending height. The calculation results using feature points extracted by OpenCV fluctuate more, and the results of this paper have higher stability.

3.4. Terrain Reconstruction

Regarding terrain reconstruction, the method of this paper adopts the plane sweep algorithm with multiple images in triangular units, whereas the previous method is the plane sweep algorithm with just image pairs in image units. In order to compare the number of initial seed points, a descent image is selected in the experiment, with all the initial seed points obtained by the method of this paper in the area covered by this image projected onto it for comparison with the initial seed points of this image calculated by Meng’s method (Figure 12) [24,25].

Clearly, in Figure 12, there are missing initial seed points in the circular area around the center of the image. This region, in which the variations in parallaxes are not sensitive to changes in virtual plane depth, is actually the region near the epipole of the image, resulting in a score curve that is too flat for pixels to meet the seed point selection condition. Overall, it seems that the proposed method in this paper yields a much larger number of seed points, even in the region near the epipoles.

In order to compare the terrain reconstruction error, the true terrain point clouds are projected to a selected descent image to obtain the ground-truth (GT) depth map. The depth maps of two methods are generated by projecting point clouds reconstructed by the methods of this paper and Meng’s to the selected descent image. Then, the depth error maps of two methods are evaluated through making difference with the GT depth map respectively. (Figure 13).

As shown in Figure 13, the depth error maps obtained by two methods both show a contour-like distribution due to the plane sweep algorithm, in which the error variation estimated by the method in this paper is smoother. In the area near the epipole, the error obtained by Meng’s method reaches more than 0.40 m, whereas the error evaluated by the method in this paper does not exceed 0.25 m. Apart from the area near the epipole, the method in this paper also calculates a lower error, which indicates that the accuracy of terrain reconstruction has been significantly improved in this paper.

The average value of depth error maps of twelve images selected by equal interval are obtained for comparison between the methods of this paper and Meng’s (Figure 14) [24,25].

It can be seen from Figure 14 that the error of the depth map obtained by both methods has a tendency of first increasing and then decreasing with the landing process. The reason for the increase is probably the gradual decrease in the tilt angle of the photographic baseline, resulting in a small swing of the ray that creates a large variation in the vertical direction, whereas the later decrease is caused by the increasing ground resolution of the image when approaching the surface. During the whole landing process, it is illustrated that the advantage of the method proposed in this paper increases with a lower acquisition height of images. This is because at a low acquisition height, the method of this paper exploits multiple images to score the correlation instead of just two images used by Meng’s method. The other images acquired from higher altitude have larger tilt angles which alleviate the insensitivity of depth changes to parallax variation.

4. In-Orbit Data Processing

In the Tianwen-1 mission, using the TR-TR algorithm proposed in this paper, the corresponding positions and poses of the Tianwen-1′s landing camera during the EDL were recovered by processing the descent images, based on which the descent velocity was analyzed. Finally, the DEM and DOM of the landing area of Tianwen-1 were generated.

4.1. Footprint Map

Based on the descent images of Tianwen-1, a footprint map was created. The footprint map shows the projection field of each descent image (the broader the field, the higher the acquisition height), which is more intuitive for us to understand the sampling method of the descent images. It is evident from Figure 15 that in the first half of the descent, there was a large swing of the landing camera which swung in the changing direction of mainly east–west and slightly north–south. A more violent swing of the landing camera denotes more deformation of the image and a greater distance from the landing site to the intersection of the camera’s optical axis with the surface of Mars. Besides, it also leads to a larger area far from the landing site covered by an image and a smaller overlap area with other images. In the second half of the descent, the pose of the landing camera gradually stabilized with a smaller swing in a mainly east–west and slightly north–south direction, along with the gradual disappearance of the deformation of the image. As the center of the image gradually approached the landing site, the overlapping area of adjacent images gradually becomes larger, with the north direction of the camera gradually shifting from the northwest to northeast, indicating that a clockwise rotation in the horizontal plane has occurred in the lander. It was found that descent images cover mainly a large area to the east, west, and south of the landing site, as well as a small area to the north of the landing site. The last few images taken before landing demonstrated that the camera was no longer swinging and no longer rotating horizontally during the final stage of descent. The upper borders of the descent images were more densely distributed, whereas the lower borders were more sparsely distributed, indicating that the horizontal moving component was directed from southwest to northeast at the final stage of descent.

The descent trajectory and landing area terrain reconstructed from the descent images of Tianwen-1 are in the camera coordinate system corresponding to the first descent image, in which case the space relationships between the cameras are based on a hypothetical scale. Therefore, the result needs to be registered under the Mars geographic coordinate system, for which purpose a total of 24 ground control points were selected, as shown in Figure 16.

Figure 16 shows that the distribution range of ground control points gradually tightens from the southern region of the landing site toward the north, because the coverage of the descent image gradually shrinks with the landing process from a large area south of the landing site to the small area around the landing site. Most of the ground control points are mainly located in the southern region of the landing site, with few in the northern region, leading to an asymmetrical distribution from north to south but a basically symmetrical distribution in the east–west direction. The reason for this is the smaller imaging range of the landing camera covering the northern region of the landing site (refer to the footprint in Figure 15). In addition, the topographic information of the landing area can be obtained in advance based on the DEM generated from the images acquired by the Tianwen-1′s orbiter. Compared with the rugged topography of the northern region of the landing site, the southern region is flatter, which can reduce the risk of obstacle interference when the Zhurong rover carries out inspection and exploration work. Hence, the accurate acquisition of topographic information in the southern region of the landing site is more conducive to decision-making departments on the Earth to provide data support for the subsequent Zhurong rover’s path planning, which is also an important reason why the ground control points are mainly located in the southern region of the landing site.

4.2. Recovery of Trajectory and Velocity

The positions and poses of the recovered landing camera during the EDL in the Mars geographic coordinate system are shown in Figure 17.

In Figure 17, the landing camera is represented by a polyhedron, so that not only the descent trajectory of the landing camera can be seen, but also the camera poses corresponding to each descent image. In order to facilitate the display, the simplified models of the landing camera corresponding to the lowest eight descent images are reduced to a certain size. This is due to the fact that the interval heights of those eight images became significantly smaller because of decreasing landing velocity. As indicated in Figure 17, there are three adjacent descent image positions with large intervals, namely, View1 (V1, and the same next) and V2, V11 and V12, and V12 and V13. The reason is that the sampling intervals of the descent images at the other positions are 4.096 s, whereas the sampling intervals of the camera became longer when these few images were taken; the sampling interval is 28.672 s between V1 and V2, 15.737 s between V11 and V12, and 12.936 s between V12 and V13 (refer to Section 2.1). Therefore, this phenomenon does not illustrate that the descent velocity has changed dramatically at the positions corresponding to these few images. The analysis about the footprint can be verified by Figure 17, the lander swings more in the east–west direction and less in the north–south direction during the EDL. The swing becomes negligible during the descent process corresponding to the last eight images, especially the swing in the east–west direction. During the time frame of the 20 descent images in total, the lander descended over 5000 m, spanning a distance of about 800 m from south to north and about 1100 m from west to east, for a total of about 1360 m in the horizontal direction, with the horizontal movement component pointing from eastward to northeastward.

The quadratic curves in two directions are fitted respectively. In the east-facing View, the equation of the quadratic curve is (4-bit valid digits retained with integer bits less than 4):

$$z=-0.002350y^2+7701y-\mathrm{5,715,107,590}$$

(20)

In the north-facing View, the equation of the quadratic curve is (4-bit valid digits retained with integer bits less than 4):

$$z=-0.002950x^2-34.53x-\mathrm{98,627}$$

(21)

In order to verify that the descent trajectory has the property of parabola, the residual error in the XY-plane between the fitted curves and the actual trajectory value is obtained as in

Table 4. Residual error in the X-axis direction and Y-axis direction.

Z (m)	Res in X-axis (m)	Res in Y-axis (m)	Z (m)	Res in X-axis (m)	Res in Y-axis (m)
2508.252	8.213	3.935	−1482.290	5.689	8.766
828.214	−8.655	−0.465	−2542.272	6.976	0.182
591.026	−8.452	−0.334	−3094.474	0.944	−1.198
301.692	−1.730	−3.466	−3210.715	0.092	−5.491
130.631	−3.613	−8.914	−3324.941	−0.003	−1.711
−161.540	−0.994	−6.517	−3437.747	−0.582	−2.156
−404.045	−2.488	−2.847	−3545.003	−0.247	−0.095
−657.444	4.802	3.305	−3647.823	−1.170	2.323
−896.076	3.575	5.838	−3738.504	−2.928	0.689
−1213.155	4.299	7.291	−3822.342	−3.728	0.865

.

Statistically, the root mean square error (RMSE) and the coefficient of determination of the recovered trajectory are obtained, respectively (

Table 5. RMSE and coefficient of determination in the X-axis direction and Y-axis direction.

Statistic Variate	X-axis	Y-axis
RMSE (m)	4.458	4.351
Coefficient of Determination	0.903	0.967

).

In the vertically downward direction, the descent trajectory can be seen as shown in Figure 18.

As illustrated by Figure 18, the overall landing trajectory is relatively smooth, except for slight fluctuations in the middle section. During the first half of the descent, the horizontal component of the velocity was mainly pointing from west to east, whereas in the second half of the descent, a slight left turn occurred, after which the direction of the horizontal component pointed from the southwest to northeast. The horizontal interval variation of the camera corresponds to the vertical interval variation of the camera in Figure 17. The V1 and the V12 far away from adjacent images, which is also the same situation in Figure 17. In the final stage of descent, the green dots are densely distributed, no longer turning and oscillating when facing northeast, which corresponds exactly to the way the last few images in Figure 15 overlap.

According to the translation of the landing camera, combined with the acquisition time of each descent image, the velocity components along the three directions of the coordinate axis were calculated, from which the velocity components in the XY plane (horizontal plane) were obtained, as depicted in Figure 19.

As displayed in the Figure 19, the velocity along the x-axis underwent a large fluctuation in the middle section with two extreme maximum values and two extreme minimum values, and gradually decreased to almost zero in the second half with a smooth trend. There were two extreme maximums and an extreme minimum of the velocity along the y-axis which experienced three significant decreases and two rebounds. Besides, the initial velocity along the y-axis was not as large as along the x-axis. Small fluctuations were observed at the end of the landing, presumably due to the thrust that was mostly along the y-axis. The thrust was provided by the thrusters during the auto-tuning process. From Figure 19, it can also be known that the horizontal velocity in the first half was mainly determined by the velocity along the x-axis, whereas at the end of the landing it was greatly influenced by the velocity along the y-axis. The initial velocity along the vertical direction was larger and relatively stable, reaching more than twice the initial velocity in the horizontal direction. After experiencing fluctuations in the first half, the vertical velocity reached its maximum at an altitude of about 3000 m above the surface. Then, the velocity started a rapid descent until it was about 30 m/s, during which the lander descended by an altitude of about 2000 m. At the end of the descent trajectory, a slight change in acceleration occurred as the thrusters were turned on about 1000 m above the surface, which corresponds to a velocity oscillation in the y-axis direction before landing.

4.3. Terrain Reconstruction

The DEM (Figure 20) and DOM (Figure 21) were reconstructed from the descent images of Tianwen-1.

This DEM of the landing area spans about 1600 m in the east–west direction, corresponding to a longitude range of about 2′, and about 1600 m in the north–south direction, corresponding to a latitude range of about 2′, which mainly contains an area of about 1200 m south of the landing site, an area of about 400 m north of the landing site, and an area of 800 m to both the west and east of the landing site. It is demonstrated that the topographic height difference in the area is within 20 m, which implies that the overall topography is relatively flat with lower terrain in the northwest and southeast areas and higher terrain in the southwest and northeast areas. The morphology of the impact crater is mainly represented by a high circular crater lip and a relatively low crater center, according to which several larger impact craters (five located in the western part with one in the eastern part as marked in Figure 20) in the landing area are clearly distinguished in the DEM. In order to evaluate the accuracy of DEM, the RMSE in the vertical direction of the control points was calculated to be 1.612 m.

To facilitate comparison with the DEM, the DOM is produced to the same extent as the DEM (Figure 21). Since the sun was located to the east of the lander during the descent process, the east side of the descent image was brighter and the west side was darker, which was most acute at the edge of the image. In addition, affected by the combination of the vignetting of the landing camera’s lens and light propagation attenuation, the descent image shows a phenomenon of brightness in the middle and darkness around, which exacerbates the problem of low illumination on the western edge of the descent image. For this reason, two templates were used for the illumination filtering of the image. Both of them are dark in the middle and bright at the periphery with the relation between light intensity and distance from the center of the template conforming to the quadratic functions. The first template is centered on the projection of the sun in the image, with the second template being centered on the center of the image. This filtering greatly attenuates but does not completely eliminate the difference in the edge’s brightness in the point cloud fusion process. As a result, some of the edges of the descent images are still visible in the DOM, which is not yet enough to interfere in the recognition of the feature’s information on Mars’ surface. As described by the DOM, there are a large number of sand ridges that mainly exist in the flat area and rarely lie inside the impact crater or across the crater lip in the landing area. The sand ridges are mainly arc-shaped, pointing from northwest to southeast with the convex direction facing north. Most of the sand ridges exist independently, whereas some of them appear to be connected. In addition, impact craters of different sizes and depths in the landing area can be clearly seen, which together with sand ridges pose a potential safety hazard to the roving rover. On this account, their effective identification will assist teleoperation departments on Earth in route planning to avoid danger.

5. Discussion and Conclusions

5.1. Discussion

In this paper, the motivation of proposing the TR-TR algorithm for dual restrained conditions comes from the mission support for Tianwen-1. When using the descent images to recover the descent trajectory and reconstruct the terrain of the landing area, various shortcomings of previous algorithms were found, thus failing to robustly recover the trajectory or reconstruct the accurate terrain. In view of this problem, this paper has made relevant research and found that there are dual restraints between the descent trajectory and the terrain of the landing area. The terrain condition in the landing area leads to the inability to calculate the motion of the landing camera robustly, and the forward motion of the landing camera along the parabolic trajectory leads to the inability to reconstruct the accurate terrain. Aiming at the problem of dual restraints, this paper studied the key technologies that can solve it, and proposed the TR-TR algorithm. The TR-TR algorithm has been verified and improved by field experiment and then applied to in-orbit data of Tianwen-1.

The advantage of the TR-TR algorithm is that it can solve the dual-restraints problem at the same time, while previous work has not been able to propose effective solutions to the dual-restraints problem, or only partial solutions to one of them. For example, data obtained by sensors were utilized to overcome the unrobustness of the trajectory recovery due to the flat terrain [19]. The plane-sweep algorithm with BSFP strategy was proposed to solve the problem of low terrain reconstruction accuracy caused by the forward motion of the landing camera [25]. Based on these works, the TR-TR algorithm made improvements for each restraint, and finally integrated them together to become an algorithm that can overcome both restraints. Thus, the TR-TR algorithm outperforms previous algorithms even when dealing with a single restraint. Some of the other previous algorithms were also tested in the field experiment. Under the condition of a single restraint in the TR-TR algorithm, the maximum error of the recovered trajectory is 0.712 m, and the accuracy is better than that of the original algorithm, but it is not as high as that of the TR-TR algorithm (0.397 m) for dual restraints [19]. The reason may be that TR-TR uses denser feature points and more accurate matching points.

The TR-TR algorithm was mainly proposed for the mission of Tianwen-1; hence, it is good at processing similar data with Tianwen-1 descent images. These kinds of descent images are sequence images, but the sampling frequency is low, which leads to the sparse images that cannot be processed according to the idea of video processing. By contrast, descent images acquired by other planet exploration missions do not necessarily have such conditions. For instance, the sensors of Chang’e-4 acquired dense descent images with high sampling frequency [27]. In addition, the descent images of Chang’e-4 have high contrast illumination since there is no air on the lunar [31]. The applicability of the TR-TR algorithm in the descent images of other planets will be better. When encountering densely sampled descent images or videos, the TR-TR algorithm will achieve higher accuracy, but also increase the computational complexity. Hence, the generality of the TR-TR algorithm still needs improvement.

The results of the field experiment showed that the TR-TR algorithm did outperform the previous methods. The TR-TR algorithm obtained refined motion with higher accuracy compared to the SIFT algorithm. This is because the TR-TR algorithm obtained uniform and dense feature points by reducing the contrast threshold. Of course, it must also have something to do with the high matching accuracy. The accuracy of refined motion is lower than that of initial motion because the cumulative error is added to it. Compared with Meng’s method [24], the TR-TR algorithm obtained more stable and higher matching accuracy after removing mismatching points. The reason for the stability is mainly because the TR-TR algorithm exploited more robust information (homography matrix and maximum height difference of the terrain), while Meng’s method exploited unstable information (fundamental matrix). More specifically, it is possible to obtain the incorrect fundamental matrix, but it is not possible to obtain the incorrect homography matrix using (quasi-) degenerated data. This is exploited by this algorithm to search for the robust homography matrix first. In fact, the homography matrix has also been combined with RANSAC to remove mismatching points [23]. However, the inlier threshold of RANSAC is not given. The threshold is explicitly given in the TR-TR algorithm based on prior terrain information. By using the DEM of the landing area, the robust maximum height difference of the landing area is known. Then, the robust homography matrix is combined with the maximum height difference to get the robust maximum plane induced parallax. Finally, matching points are checked and eliminated according to the robust maximum plane induced parallax. In addition to stability, the TR-TR algorithm also has higher matching accuracy, because the TR-TR algorithm uses the monotonicity of feature scale variation in the descent images to eliminate some other mismatching points, which is not used by Meng’s method. The field experiment of robust motion estimation mainly includes two parts: the estimation of initial motion and the adjusted motion. As for the initial value of motion, it can be seen that Meng’s method will randomly get the error solutions with the gross error when the (quasi-)degenerate data are used [24]. The error solutions come from decomposition of the incorrect fundamental matrix which is obtained randomly. In comparison, the TR-TR algorithm is able to obtain a more stable motion solution for two reasons. The first one is that the TR-TR algorithm is able to obtain two motion solutions based on the robust homography matrix decomposition and uses the parabolic property of the descent trajectory to eliminate the one that does not match the real situation, leaving the robust motion solution. The second is that according to the robust homography matrix, the TR-TR algorithm makes full use of the points outside the scene plane, and uses the epipoles constrained parallax beam method to search for a robust fundamental matrix among the matching points that may contain mismatching points. For the adjusted motion, the TR-TR algorithm is compared with the SIFT method with default threshold in OpenCV, and higher trajectory recovery accuracy is obtained. This is because the TR-TR algorithm obtains more and uniformly distributed matching points by reducing the contrast threshold of SIFT, which provides more sufficient constraints, especially for the (quasi-) degenerated data. For example, for (quasi-) degenerated data, there are points outside the scene plane, but if the number of matching points is small, it is likely to be all in the plane. However, if the number of matching points is large enough and the distribution is uniform, it is more likely to contain points outside the plane. Such a point as an observed value participating in the adjustment can provide the constraint in the vertical direction and obtain a more accurate motion solution. Compared with Meng’s algorithm [25], the TR-TR algorithm obtained higher terrain accuracy. The most direct reason is that the TR-TR algorithm found more initial seed points as exact depth constraints. The better fit of the multi-directional virtual planes to the small undulations of the terrain resulted in more initial seed points and made more ZNCC curves eligible to be the initial seed points. Meng’s method is proposed for two images, while the TR-TR algorithm is for multiple images. In the case of two images, one of the images has only one epipole, and the insensitivity of depth to parallax near the epipole cannot be resolved by more images. In contrast, in the case of multiple images, one of the images can have more than one epipole, and the insensitivity of depth recovery near a certain epipole can be solved by the sensitivity of other images in this region, thus improving the accuracy of terrain reconstruction near the epipole.

The in-orbital descent images data of Tianwen-1 were used by the TR-TR algorithm to generate the descent trajectory and reconstruct the terrain of the landing area. Some similar algorithms have been used to process in-orbit descent images of Chang’e-3 and Chang’e-4 [15,27,47]. They also successfully recovered the descent trajectory and reconstructed the terrain of landing area. The difference is that their descent trajectory is not parabolic, hence, their algorithm does not necessarily apply to Tianwen-1.

5.2. Conclusions

This paper presents a detailed description of the descent images acquired by the Tianwen-1′s landing camera and points out the difficulties encountered in the descent trajectory recovery and reconstructing terrain of the landing area using descent images in the dual-restraints conditions. In response to these difficulties, a novel trajectory recovery and terrain reconstruction (TR-TR) algorithm was proposed, which solved the problems of uneven distribution of feature points, high numbers of mismatching points owing to similar features, unrobust motion estimated from (quasi-) degenerate data, and low accuracy of terrain reconstruction in baselines along the optical axis. Robust motion was estimated by using parallax beams with epipoles constraints and utilizing parabolic characteristics of trajectories to eliminate uncertainty. More accurate trajectory was recovered by obtaining evenly distributed feature points and removing more mismatching points. More accurate and fine terrain was reconstructed by using a multi-image and multi-direction plane-sweep method with sparse spatial points. Verified by the simulated landing experiment in the Red Cliff area, the TR-TR algorithm obtained more accurate results than previous methods, including more accurate trajectories and more accurate terrain. The maximum error of recovered trajectory was reduced from 0.681 m of previous algorithm to 0.397 m. The maximum error of terrain reconstruction was reduced from the 0.692 m of previous algorithm to 0.462 m of the TR-TR algorithm. In practical applications, the TR-TR algorithm was applied to the in-orbit data of Tianwen-1 by accurately recovering the landing trajectory of the lander and mapping the high-precision terrain of the landing area. The recovered trajectory was fitted to parabola in X-axis and Y-axis directions, and the corresponding parameters and correlation coefficients were obtained. The fitting residuals were also calculated. The vertical RMSE of the reconstructed terrain in the landing area of Tianwen-1 was calculated as 1.612 m. The terrain accuracy obtained by in-orbit data is lower than that obtained by field experiment, because the accuracy of the control points in the in-orbit data is relatively lower, and the selection error of control points on the in-orbit descent image is relatively larger. The results met the accuracy requirements of the mission for terrain reconstruction and supplied a strong product to support the follow-up inspection and detection of the Zhurong rover and provided valuable data feedback for the future planetary landing mission.

The innovations of this article are as follows:

(1) A novel algorithm (TR-TR) for descent trajectory recovery and landing area terrain reconstruction based on Tianwen-1′s descent images was developed in a dual-restraints context, which has higher precision, accuracy, and robustness than previous methods;

(2) The landing simulation experiment of the Mars lander was carried out for the first time in the Mars-like landform area to verify the descent trajectory recovery and terrain reconstruction algorithms based on descent images;

(3) For the first time, the descent images of Tianwen-1 were employed to realize the recovery of the descent trajectory of the lander and the terrain reconstruction of the landing area.

5.3. Future Work

Although the TR-TR algorithm proposed in this paper has achieved good results, there are still limitations, such as the requirement of the internal parameters of the descent camera, the failure to utilize the shape and size of the ground objects in the landing area to provide scale constraints for the reconstructed terrain, and the lack of the establishment of the multi-scale precision evaluation method of terrain mapping based on images with different resolutions. Moreover, as mentioned in the discussion, the generality of the TR-TR algorithm needs to be enhanced. Further research is needed for China’s future Mars landing missions that still face enormous challenges. The terrain reconstruction algorithm based on descent images is an important prerequisite and component of point cloud based terrain reconstruction. However, dual restraints are still an inevitable problem for point cloud terrain reconstruction. Therefore, the TR-TR algorithm in this paper can also provide reference for point cloud reconstruction in planetary mapping.