Landing site positioning and descent trajectory reconstruction of Tianwen-1 on Mars

Tianwen-1 (TW-1) is the first Chinese interplanetary mission to have accomplished orbiting, landing, and patrolling in a single exploration of Mars. After safe landing, it is essential to reconstruct the descent trajectory and determine the landing site of the lander. For this purpose, we processed descent images of the TW-1 optical obstacle-avoidance sensor (OOAS) and digital orthophoto map (DOM) of the landing area using our proposed hybrid-matching method, in which the landing process is divided into two parts. In the first, crater matching is used to obtain the geometric transformations between the OOAS images and DOM to calculate the position of the lander. In the second, feature matching is applied to compute the position of the lander. We calculated the landing site of TW-1 to be 109.9259° E, 25.0659° N with a positional accuracy of 1.56 m and reconstructed the landing trajectory with a horizontal root mean squared error of 1.79 m. These results will facilitate the analyses of the obstacle-avoidance system and optimize the control strategy in the follow-up planetary-exploration missions.


Introduction
Tianwen-1 (TW-1) is the first Chinese Mars-exploration mission to have realized orbiting, landing, and patrolling simultaneously. Its success of TW-1 made China the third country after the former Soviet Union and the United States to achieve a Mars-landing exploration [1][2][3]. The probe was launched on July 23, 2020 and the spacecraft was captured by Mars on February 10, 2021 after seven months of transit. For the next three months, the probe conducted reconnaissance for choosing the target landing site. On May 15, 2021, the lander separated from the orbiter and landed in the preselected landing zone in the southern part of the Martian Utopian Plain. Landing process is the most critical aspect of a Mars exploration. For a better understanding of the autonomous-control results and analysis of the obstacle-avoidance strategy used, it is essential to reconstruct the descent trajectory and locate the landing site of lander after safe landing.
Several studies have focused on the descent trajectory reconstruction and landing site positioning of lander. For example, radio measurements [4][5][6], including radio ranging and very long baseline interferometry (VLBI), have been widely used to reconstruct landing trajectories. However, owing to the complexity and uncertainty of the Martian environment, the positioning accuracy with this method is only of the order of few kilometers. Moreover, this method cannot be applied without Earth-based radio tracking or telemetry data. Therefore, an image-based method was developed that utilizes images captured by the onboard camera along with known georeferenced data of the planetary surface. This method enables more precise calculations of the landing site and descent trajectory and is not dependent on particular exploration environments. For the Chang'e-3 (CE-3) and Chang'e-4 (CE-4) missions, a landing camera with a high frame rate of 10 frames/s, which obtained thousands of continuous images, was employed during the entire descent process with an image overlap rate of more than 94%. A number of image-based methods have been proposed to make use of these images. For instance, using the photogrammetric method, Liu et al. [7] selected 180 images in 1 s intervals from the landing images and used CE2TMap2015 as the georeferenced data for reconstructing the powered-descent trajectory. This method needs to automatically extract and match a large number of feature points throughout the overlaps of the landing images. In addition, the matched points could be input into the photogrammetric bundle adjustment to reconstruct the descent trajectory of the CE-4 lander. Taking engineering requirements into account, Wang et al. [8] determined the landing point of the CE-4 lander using high-precision image-matching and geometric-transformation methods. The landing point was initially located with near-real-time high-compression-ratio descent sequence images. Then, the estimate was refined using the replayed low-compression-ratio descent images. This method can calculate the landing point fast. In addition, in Ref. [9], control points were acquired by matching the descent images with the digital orthophoto map (DOM) generated from CE-2 images, and the descent trajectory of the CE-3 lander was obtained using the SIFT feature-matching and bundle-adjustment methods.
In this study, an image-based method was employed to fulfill the precise positioning requirements of TW-1. However, the limited transmission rate and high time interval for transmission between the Earth and Mars makes it unrealistic to download the data obtained by the landing camera. Consequently, OOAS image data were used in this study. The optical obstacle-avoidance sensor (OOAS) is one of the payloads installed on the TW-1 lander, which recorded the entire landing process in a small number of images. Note that there are two disadvantages to using OOAS images. First, the images were downloaded with compression, leading to loss of information. Second, the high sun-elevation angle during the landing of TW-1 resulted in low image contrast and poor image quality.
Considering the properties of the OOAS images, we propose a hybrid-matching approach which executes descent-trajectory reconstruction in two parts. First, for the case of the lander at a high altitude, the crater is clearly recognizable in the OOAS images, which enables us to employ the crater-matching method. In contrast, in the case of the lander at a low altitude, the number of craters in the OOAS images is limited, and the feature-matching method needs to be applied. Based on hybrid-matching method, the landing site and trajectory were determined and reconstructed smoothly, which does not rely on Earth-based telemetry data and multi-landing images. In conclusion, for the first time, we located the landing site of TW-1 as 109.9259 • E, 25.0659 • N with an accuracy of 1.56 m and reconstructed the landing trajectory with a horizontal root mean squared (RMS) error of 1.79 m. These results can aid the analysis of the obstacle-avoidance system and control-strategy optimization for the follow-up planetary-exploration missions.
The rest of this paper is organized as follows. In Section 2, we describe the dataset including the OOAS images and DOM data and define the coordinate systems. In Section 3, the technical details of the hybrid-matching method and related experimental processes are presented. Sections 4 and 5 present the procedures and results, respectively, of the landing site determination and landing trajectory reconstruction. Finally, a summary and conclusion are presented in Section 6.
2 Dataset description and coordinatesystem definition

Dataset description
To reconstruct the landing trajectory and determine the landing site, we used the images captured by the onboard OOAS and high-resolution DOM data of the Martian surface. The position of the lander at the time of OOAS image capture was calculated by detecting corresponding features between the OOAS image and DOM data. The details of the two datasets are presented in the following.

OOAS images
The OOAS is one of the payloads on the TW-1 lander responsible for ensuring safe landing. The OOAS records the entire powered descent starting from when the lander releases the heat shield (approximately 10 km from the Martian surface). The frame rate of OOAS was approximately 1 frame per 4 s, and its main parameters are listed in Table 1. The OOAS images contain position information for TW-1, which can be used to locate the landing site and reconstruct the descent trajectory. In this study, 36 OOAS images were selected for use in landing-trajectory reconstruction, as shown in the Appendix. The capture time of the image is presented in Fig. 1: The capture time of the first image is the origin, and the horizontal and vertical axes represent OOAS images 1-36 and their capture time, respectively.

DOM data
The DOM data were used as geographical reference data   for landing site positioning and trajectory reconstruction. In this study, the DOM data (shown in Fig. 2) of the landing zone were generated from the images captured by the onboard high-resolution camera [10,11]. The DOM data cover the landing area from 109 • 48 ′ E to 109 • 58 ′ E and from 24 • 58 ′ N to 25 • 6 ′ N with a spatial resolution of 0.7 m, which is available at the Data Publishing and Information Service System of China's Lunar Exploration Program [12].

Coordinate-system definition
For trajectory reconstruction, three coordinate systems need to be defined: a georeferenced coordinate system X m , an image coordinate system X i , and a camera coordinate system X c . These were defined as follows (Fig. 3). First, the origin of X m was placed at the center of the map, and the axes x m , y m , and z m defined as the east, north, and vertical directions, respectively. Second, X i was defined as a 2D coordinate system with the x i and y i axes parallel to the column and row of the OOAS image, respectively. Finally, for X c , the origin was placed at the optical center of the camera with the axes x c and y c pointing in the same directions as those in X i . The z c axis points along the main optical axis and forms a right-handed set with x c and y c .

Hybrid-matching method
Descent trajectory and landing site calculation requires a geometric transformation H nD between the lander and Martian surface, where H nD is computed using the matching-point pairs between the OOAS images and DOM data. In this paper, a hybrid-matching method based on crater matching and feature matching is proposed for this, as shown in Fig. 4  because of the low contrast and poor quality of the OOAS images, image enhancement was performed prior to the matching process. After image processing, the following procedures were applied to the OOAS images. First, for the lander at a high altitude (images 1-24), the H nD (n = 1, 2, · · · , 24) was estimated by the crater-matching method, which matches the craters detected in the OOAS images with the known crater database built from the DOM data. Subsequently, for the case of the low-altitude lander (images 25-36), the matched craters are insufficient. Thus, in this case, the H nD (n = 25, 26, · · · , 36) was calculated using the feature-matching method as follows. First, we computed the corresponding relationship T n (n = 25, 26, · · · , 36) between the OOAS images from the 25th to the 36th and the 24th OOAS image through inter-frame feature matching. With the obtained geometric transformation H 24D at high altitudes, T n (n = 25, 26, · · · , 36) was converted to the georeferenced coordinate system, and therefore, the last part of H nD (n = 25, 26, · · · , 36) was acquired.
Being one of the most common features on the surfaces of extraterrestrial bodies, craters are suitable for use in autonomous navigation and obstacle identification by spacecraft. In the crater-matching process, the corresponding points can be obtained by matching the craters extracted from the OOAS images with those in the pre-established crater database [13][14][15]. In this study, crater matching was mainly applied to the OOAS images taken at altitudes above 3 km from the Martian surface, and the matching procedure is outlined in Algorithm 1.
Step 2: Extract craters from the OOAS images and generate the list of crater triangles and triangular invariants.
Step 3: Match the list obtained in Step 2 with the crater database using triangle matching method and determine the crater matches.
First, the crater database was created based on the DOM data. The main parameters of the crater database are as follows.
(1) Crater locations and diameters of the DOM. In this study, 78 craters with radii in the range of 5.5-175.6 m were selected from the DOM. The distributions of these craters are shown in Fig. 5(a).
(2) List of crater triangles and triangular invariants. Three craters form a crater triangle, and all possible crater triangles were catalogued in a list. Each catalogued crater triangle contains triangular invariants as follows: the smallest, middle, and largest angles; cosine of the largest and smallest angles; three crater diameter ratios; and the sense of rotation for each triangle. In particular, the crater diameter ratios were calculated as the ratio of the crater diameters to the longest leg of the corresponding crater triangle. The sense of rotation of the crater triangle indicates a clockwise or counterclockwise path, which is produced by tracing the triangle from the vertex of its smallest angle to that of the middle angle and then to the largest angle [16]. As an illustration, Fig. 5(b) depicts an image of a typical crater triangle. In Fig. 5(b), the smallest, middle, and largest angles are denoted as α s , α m , and α l , respectively. The diameters of craters i, j, and k were divided by the length of R ij to  obtain the crater diameter ratios. The sense of rotation of the crater triangle in Fig. 5(b) is counterclockwise. Second, the craters were extracted from the OOAS images, and a list of crater triangles and triangular invariants for the OOAS images was made (the OOAS list of crater triangles and triangular invariants was the same as that of the crater database).
Third, using the triangle-matching method (see Refs. [16,17] for a description of the algorithm), the crater triangle list of OOAS images was matched with that of the crater database. Thus, we obtained the matching crater pairs between the OOAS images and DOM. The matching crater pairs of the 1st and 24th OOAS images are shown in Fig. 6.
Finally, using the centers of the crater matches, we obtained a group of matching point pairs between the image coordinates and georeferenced coordinates. The geometric transformations H nD (n = 1, 2, · · · , 24) between the OOAS images and DOM were computed by matching point pairs [18]. With decreasing spacecraft altitude, crater matching becomes unreliable owing to the decreasing number of matched craters. For the 25th to the 36th images, which were captured below 3 km, the feature-matching method was used to construct the geometric transformations H nD (n = 25, 26, · · · , 36) between the OOAS images and DOM. The feature-matching process is shown in Algorithm 2.
The feature-matching process can be described as follows. First, by selecting and matching the point Algorithm 2 Process of feature matching Input: OOAS images (24-36), geometric transformation H 24D between the 24th OOAS image and DOM Output: Geometric transformations H nD (n = 25, 26, · · · , 36) between OOAS images and DOM data Algorithm: Step 1: Select and match point features between OOAS images n and n−1. The matching point pairs p n and p n−1 can be acquired.
Step 2: Compute homography matrices H n(n−1) from the matching points in Step 1.
Step 3: Align H n(n−1) with the same coordinates, which is the camera coordinate system of the 24th OOAS image taken. Subsequently, the unified homography T n can be determined.
Step 4: Convert the unified homography T n to the georeferenced coordinate by H 24D . The geometric transformations H nD (n = 25, 26, · · · , 36) between the OOAS images and DOM are obtained. features [19] from the OOAS images n and n − 1 (n = 25, 26, · · · , 36), the matching point pairs p n and p n−1 between the adjacent OOAS images were obtained. The matching points between OOAS images 24 and 25 are shown in Fig. 7. Second, homography matrices were computed from the matching point pairs. The homogeneous coordinates of p n and p n−1 are P n and P n−1 , respectively. The corresponding relationship is as Eq. (1): where H n(n−1) can be solved by the least-squares method when there are multiple sets of matching point pairs. Third, the 24th OOAS image and its camera coordinate system were used as the reference image and reference coordinate system, respectively. The OOAS images can be aligned with the reference image as Eq. (2): where T n is the unified homography between the camera coordinate system of images 25-36 and the reference coordinate system. With T n , we projected the OOAS images (25-36) onto the 24th OOAS image. The projection results are shown in Fig. 8. Finally, with the geometric transformation between the reference coordinate system and georeferenced coordinate system H 24D , the geometric transformations H nD between the OOAS images (25-36) and DOM can be acquired by By combining the geometric transformations obtained by crater matching and feature matching, we obtained all geometric transformations H nD (n = 1, 2, · · · , 36), which were then used in the following landing site positioning and descent trajectory reconstruction processes.

TW-1 landing site positioning
Because the lander adopted a vertical-descent method for approaching the landing zone during the final period of the descent, we set the center of the last OOAS image as the landing site. The OOAS images were projected onto the DOM with the geometric transformations H nD (n = 1, 2, · · · , 36), as shown in Fig. 9, where the center of the last OOAS image is marked with a red cross in DOM. According to the range of DOM, the longitude and latitude of the landing site were confirmed as 109.9259 • E and 25.0659 • N. To verify the accuracy of the landing site location, the result was compared with the location of landing site seen in the high-resolution image from the China National Space Administration [20]. Using the four craters near the landing site as references, the landing site obtained using the proposed method was seen to be consistent with that in the real image, as shown in Fig. 10. Moreover, when the high-resolution image was projected onto the DOM (Fig. 11), the calculated landing point differed by only 1.56 m from that observed in the photograph.

TW-1 descent trajectory reconstruction
To evaluate the autonomic-control strategy and obstacle-avoidance method, the landing process of TW-1 was recovered and is described in this section. With the geometric transformation H nD = [h 1 , h 2 , h 3 ] (n = 1, 2, · · · , 36) obtained in Section 3, the rotation matrix where h 1 , h 2 , h 3 and r 1 , r 2 , r 3 are the column vectors of H nD and R, respectively. K is the internal camera parameter matrix, which is obtained from Table 1. z c is a scale factor computed from Eq. (4).
Thus, the position of the camera in the georeferenced coordinate system T ms is According to the installation relation T si between the OOAS and the lander, the landing trajectory of the lander in the georeferenced coordinate system X is obtained as Eq. (6): The reconstructed trajectory is shown in Fig. 12, and Fig. 13 shows the variation in lander altitude with time.  to 6083.48 m above the Martian surface, TW-1 moved 1200.00 m to the southeast. Subsequently, TW-1 adjusted its direction to the east and continued to fly for 1216.50 m. When the lander altitude decreased to 1353.26 m, the backshell was ejected and backshell-avoidance algorithm executed. To avoid the ejected backshell, TW-1 moved 285.00 m to the northeast between the altitudes of 1353.26 and 100.00 m and entered the hover and precise obstacle-avoidance phase. At an altitude of 100.00 m, TW-1 hovered and selected a safe landing site. Then, the lander moved to the northeast and crossed a crater with a diameter of 18.00 m, during which the altitude decreased to 69.57 m. Finally, TW-1 descended slowly and landed safely. As seen in Fig. 14, the horizontal distance between the landing points of the lander and backshell was 350.00 m and the corresponding distance obtained by the navigation camera was 350.82 m [21]. This shows that the landing trajectory was obtained with an acceptable precision.
To verify the accuracy of the determined trajectory quantitatively, the RMS of the reprojection error and the horizontal error of the ground control points (GCPs) were calculated [22]. In this study, 19 ground control points (GCPs) near the landing site were manually selected on the DOM and OOAS images according to the principle that the preferred choices are smaller craters and rock centers and those distributed evenly in the landing area. The distributions of the GCPs in this study are plotted in Fig. 15. With the geometric transformations calculated in Section 3, the GCPs in the DOM were projected onto the OOAS images, and the corresponding projected positions in the OOAS images are regarded as estimation values. The GCPs manually selected in the OOAS images were considered as the true values. The difference between the true and estimated values corresponds to the reprojection error. The horizontal error of the GCPs is the result of multiplying the OOAS image resolution with the   Table 2. From Table 2, we can see that the reprojection error of the control points is 1.42 pixels, and the horizontal error is 1.79 m, which is considered to be the accuracy of the TW-1 descent trajectory determination.

Conclusions
In this study, a hybrid-matching approach was employed to precisely locate the landing site and reconstruct the trajectory of TW-1 using sequential images captured by its OOAS. First, the precise location of the TW-1 landing site was determined to be 109.9259 • E, 25.0659 • N within an accuracy of 1.56 m. This can serve as geographic data support for subsequent Mars-exploration missions. Furthermore, the landing trajectory of TW-1 was reconstructed, and the horizontal error of the trajectory was verified by the selected control points as being 1.79 m. The reconstructed trajectory can be applied to analyze and optimize obstacle-avoidance and autonomous-control strategies. In the future, with the release of higher quality images and greater scientific data, more automatic methods can be utilized to obtain more precise results with increased efficiency.