Topographic Reconstruction of the “Tianwen-1” Landing Area on the Mars Using High Resolution Imaging Camera Images

High-resolution optical cameras have always been important scientific payloads in Mars exploration missions, and the Mars topographic data produced by their detection data can provide support for scientific research on Mars geomorphology and geological structure evolution. As of December 2021, there are still relatively few high-resolution image data at the submeter level on the Martian surface, with about 2.6% global coverage and even more limited stereo coverage (just about 0.4%). At the same time, there are still some difficulties in data acquisition and terrain reconstruction processing methods for high-resolution Mars images that need to be solved. This article described how we designed the in-orbit stereo imaging strategy based on the characteristics of the high-resolution imaging camera (HiRIC) of China’s first Mars exploration mission (Tianwen-1), studied the technical solutions for HiRIC stereo image photogrammetry processing, and produced a topographic dataset for the “Tianwen-1” landing area, including a digital orthophoto map (DOM) with a ground sample distance (GSD) of 0.7 and 3.5 m and a digital elevation model (DEM) with a GSD of 3.5 m. Precision analysis results show that these topographic data have good consistency in planar position and elevation compared with the existing Mars terrain data and have advantages in spatial resolution and terrain detail expression, which will be widely used in the geological background study of the “Tianwen-1” landing area, as well as landing site positioning, Martian surface remote operation planning, and other Mars scientific research and engineering tasks.


31
M ARS topographic information is the most important 32 data to quantitatively describe the basic characteristics 33 of the Martian surface, and it is indispensable for scientific 34 research on surface topography, mineral and rock distribution, 35 geological structure and evolution, gravity field and internal 36 structure, atmospheric dynamics, and so on. At the same time, 37 the high-precision topographic data can also support engineer-38 ing tasks, such as probe orbit optimization design, landing site 39 positioning, Mars rover exploration path planning, and so on. 40 Therefore, the detailed reconstruction of Mars terrain is an 41 important and fundamental work in Mars exploration missions. 42 High-resolution satellite images and laser altimeter have 43 been two major categories of datasets for topographic mod-44 eling. To date, international Mars exploration missions have 45 carried a variety of high-resolution optical cameras and laser 46 altimeters, such as the Mars orbiter camera (MOC) [27], [29] 47 and the Mars orbiter laser altimeter (MOLA) [39] on Mars 48 global surveyor (MGS), the high-resolution imaging science 49 experiment (HiRISE) [30] and context camera (CTX) [28] 50 on Mars reconnaissance orbiter (MRO), the high-resolution 51 stereo camera (HRSC) [15], [34] on Mars express, the color 52 and stereo surface imaging system (CaSSIS) [42] on ExoMars 53 trace gas orbiter (TGO), and so on. The images and detection 54 data acquired by these payloads have been used to produce a 55 rich set of mosaics and topographic data of the Mars [3], [16], 56 [29], [31], [33], [35], [39], [43], which accurately characterize 57 the topography of the Martian surface and provide a geo-58 graphic benchmark for Mars exploration missions, exploration 59 data processing, applications, and comprehensive scientific 60 research. 61 In terms of technical characteristics, laser altimeter has a 62 remarkable higher vertical accuracy and global consistency, 63 allowing direct access to topographic information on the 64 Martian surface; however, in the horizontal direction, the 65 strip-like arrangement of points in widely spaced profiles 66 suffers from severe differences of density, which reduce the 67 effective resolution. The current "gold standard" for Martian 68 topographic data, MOLA, has collected data globally with 69 astonishingly high accuracy, but the sample spacing of this 70 dataset is only about 300-m along track, and, in many places 71 near the equator, adjacent MOLA ground tracks are separated 72 by gaps of one to several kilometers [39]. In contrast, images 73 This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ can provide much better horizontal resolution of terrain details, and so on [45]. HiRIC has the highest ground sample distance 132 (GSD) of 0.5 m at 265 km, which is an effective supplement to 133 the current submeter high-resolution Mars image data. At the 134 same time, it can provide basic data for engineering tasks, 135 such as landing site selection, landing site positioning, Mars 136 surface exploration path planning for the "Zhurong" rover, and 137 scientific tasks, such as geological background study of the 138 landing area [26], [44]. 139 In this article, we used the high-resolution stereo images 140 of the "Tianwen-1" landing area acquired by HiRIC during 141 the Mars parking orbit to produce the topographic dataset 142 with a coverage range of 150 × 60 km by photogrammetry 143 processing, including a digital orthophoto map (DOM) with a 144 GSD of 0.7 and 3.5 m and the DEM with a GSD of 3.5 m. 145 The absolute position of this dataset is in good agreement with 146 the existing Mars topographic data, such as the CTX mosaic 147 [6] and the HRSC/MOLA Blended DEM [8], and can present 148 richer topographic details to meet the needs of future scientific 149 research and engineering tasks. HiRIC adopts an off-axis three-mirror astigmatic (TMA) 153 optical system with a focal length of 4669.794 mm, 154 an F-number of 12, and a field of view (FOV) of 2 • × 155 0.635 • . It contains three linear-array time delay and inte-156 gration charge-coupled devices (TDI CCDs) and two-area-157 array complementary-metal-oxide semiconductors (CMOSs) 158 on one imaging focal plane to produce push-broom and frame 159 imaging, respectively. TDI CCDs can work in panchromatic 160 and four spectral bands, including blue (B1), green (B2), 161 red (B3), and near infrared (B4), so that HiRIC can obtain 162 both panchromatic and color images. The pixel number for 163 the panchromatic band of each TDI CCD in the cross-scan 164 direction is 6144 with a pixel size of 8  Because of the highest GSD of 0.5 m at 265 km on the Martian surface, panchromatic band CCD push-broom imaging will be the default in-orbit detection mode for HiRIC [45].   Fig. 2) reaches 2 • , periareon orbital maintenance is 204 required. The periareon orbital arc of the Mars parking orbit 205 transits over the landing area of "Zhurong" rover from north to 206 south, whose center coordinates are the nominal landing site. 207 Before the rover landed on the Mars, HiRIC performed high-208 resolution stereo imaging of the landing area in the periareon 209 orbital arc. Using these stereo images, topographic dataset 210 of the landing area was produced through photogrammetric 211 processing within one month.

212
The landing area of the "Zhurong" rover is 150 km 213 (north-south direction) × 60 km (east-west direction), 214 as shown in Fig. 2(b). Since the imaging swath width of 215 a single-track push-broom image is only 9 km at 265 km, 216 to meet the demand of topographic data production, HiRIC 217 needs to image the same area of the landing area along the 218 north-south direction with two angles from the forward and 219 backward views in adjacent two orbits. The imaging baseline-220 height ratio should be ≥0.4, i.e., the "Tianwen-1" orbiter needs 221 to adjust ±11.31 • in pitch direction for forward and backward 222 views [ Fig. 2(a)]; in addition, the overlap coverage of adjacent 223 images along the east-west direction should be no less than 224 15%, so that the orbiter needs to adjust its attitude in roll 225 direction.

226
Based on the analysis, we divided the landing area into nine 227 imaging strips ( Fig. 2(b); see also [45]), each with an east-west 228 width of about 10 km and a north-south length of about 229 160 km, with an overlap rate of 15% between adjacent imaging 230 strips. During the in-orbit detection, HiRIC will complete the 231 forward-and backward-view imaging of each imaging strip 232 one by one. Before each imaging session, according to the 233 orbital altitude of the orbiter during its transit over the landing 234 area, the distance of the subsatellite point of the orbiter from 235 the center of the imaging strip and other information, the 236 attitude angles of the orbiter in pitch and roll direction during 237 HiRIC imaging are solved, together with its yaw angle. The 238 orbiter attitude is controlled to align the HiRIC FOV with 239 the imaging area, and the imaging parameters settings are 240 completed to achieve high-resolution imaging of the imaging 241 strips, such as integration level, image motion compensation 242 parameters, and so on.

243
On 26 and 28 February 2021, test imaging sessions were 244 conducted to verify the in-orbit operating status, image quality, 245 and stereo imaging strategy of HiRIC. After fully verifying 246 the data quality, HiRIC completed a total of 18 formal imag-247 ing sessions for the nine imaging strips from forward and 248 backward views to achieve stereo coverage during 2 March 249 to 10 April, with each imaging lasting about 42 s, imag-250 ing swath width about 11 km, and north-south length of 251 about 170 km. Since the periareon of the Mars parking orbit 252 shifts eastward over time and moves away from the nominal 253 landing site, the roll angle of the orbiter is different each 254 time HiRIC images the specified imaging strip (as shown in 255  Table II), and the farther the periareon is from the center of the 256 imaging strip, the larger the roll angle required for imaging. 257 According to the analysis results of HiRIC imaging quality 258 and the orbit determination accuracy of the "Tianwen-1" 259 orbiter, the roll angle of the orbiter should not exceed ±20 •   Table II), which can satisfy both the consistency of the  [40]).

293
The level 2B data product of HiRIC images, after relative 294 radiometric correction and geometric positioning [40], is the 295 main data source for this article. In addition, the orbital 296 ephemeris and attitude data of the "Tianwen-1" orbiter are 297 used to solve the external orientation elements of the HiRIC 298 images. The orbital ephemeris of "Tianwen-1" orbiter is 299 provided by the Beijing Aerospace Control Center (BACC), 300 using the J2000 coordinate system, with the precision orbiting 301 accuracy better than 1 km (3σ ) and the extrapolated orbit 302 accuracy better than 3 km (3σ ), in which the precision orbital 303 ephemeris is used for data processing in this article; the atti-304 tude data of the orbiter are obtained by telemetry parameters, 305 using the J2000 coordinate system too, with the measurement 306 accuracy better than 0.05 • (3σ ) and the stability better than 307 0.0005 • /s (3σ ), which will cause a positional deviation of 308 280 m from 320-km altitude.  tures [41]. Tie points between the forward-and backward-375 view images within the same imaging strip and in the overlap 376 area between adjacent imaging strips are extracted to construct 377 the affine transformation function between images, and then, 378 exact matching and coarse difference rejection on the initial 379 tie points is performed [13], [14]. These tie points are used as 380 adjustment points for photogrammetry processing. At the same 381 time, the CTX mosaic and the HRSC/MOLA Blended DEM 382 mentioned earlier are used as the planer position and elevation 383 data, respectively, to extract the adjustment control points. 384 Uncontrolled block adjustment technology of optical remote 385 sensing images based on "alternating direction method of mul-386 tipliers (ADMM)" [4] and the global least squares adjustment 387 method is adopted to optimize RFM model parameters and 388 affine transformation function coefficients and to realize seam-389 less mosaic and high precision positioning of HiRIC images. 390 Using the optimized RMF model and affine transformation 391 function, the epipolar line constraint relationship of the stereo 392 image pairs is established, and the pixel-by-pixel matching of 393 homologous points of the HiRIC forward-and backward-view 394 stereo images is performed to obtain the dense point cloud 395 data [47]. Inverse distance weighted (IDW) is used to produce 396 a regular grid DEM with a GSD of 3.5 m based on the dense 397 point cloud data. Finally, a DOM with a GSD of 3.5 and 0.7 m 398 was produced using DEM and HiRIC images.

400
The RFM is a generalized expression of various sensor 401 imaging geometric models, which is independent of specific 402 sensors, simple in form, and so on. It can effectively overcome 403 the shortcomings of traditional, strict imaging models based on 404 colinear equations with strong correlation between orientation 405 parameters, more solution parameters, and unstable numerical 406 solutions, and meet the requirements of transparency of sensor 407 parameters, generalization of imaging geometric models, and 408 high-speed intelligence of processing, which has been widely 409 used in sensor model calculation. The RFM is fixed based on 410 a priori orientation data, and the adjustment will be done in 411 terms of affine functions in image space.

412
In the RFM model of HiRIC, the relationship between 413 image pixel coordinates and their corresponding ground coor-414 dinates is expressed by rational polynomials, where all coordi-415 nates are normalized. For descriptive convenience, let x n and 416 y n be the normalized image pixel coordinates, and (ϕ n , λ n , h n ) 417 be the normalized geographic coordinates (longitude, latitude, 418 and ellipsoidal height); then, the mathematical expression of 419 the RFM model is as follows: where f t (ϕ n , λ n , h n ), t = 1, 2, 3, 4 is a cubic polynomial, 423 which contains a total of 80 rational polynomial coefficients 424 (RPCs; the constant term of the two denominators is fixed 425 at 1.0) and 10 normalization parameters, which are called 426 RPCs. The RPCs of HiRIC are computed using a topography-427 independent solution method proposed by Tao and Hu [41]. 428 The method first uses the orbital attitude model of HiRIC 429 images to establish a set of virtual 3-D grid points as control 430 points, and then solve the RPCs by the $L2$ regularization, 431 which can achieve a high-precision RFM fitting of the HiRIC 432 rigorous geometric imaging model for the subsequent pho-433 togrammetric processing.

434
It is important to note that the "Tianwen-1" orbiter is  2) Holistic matching based on image salient edge features. 496 The basic idea of this method is Hough transform image 497 matching [14], that is, using the spatial geometric constraint 498 relationship between the features extracted from HiRIC images 499 and orthophoto images to construct affine transform functions, 500 voting in the affine transform parameter directly for fast 501 matching between images, which can quickly and effectively 502 eliminate the systematic errors existing in the original RFM 503 orientation model. The systematic errors of about several 504 hundred or even thousands of pixels between the HiRIC 505 images and the planar position data can be optimized to obtain 506 a matching result of about 20-30 pixels, and finally, the tie 507 points between the HiRIC images are obtained. The holistic 508 matching takes only a few seconds to complete, and the 509 matching results enable the initial positioning of HiRIC images 510 and provide initial values for the exact matching positioning 511 of tie and control points. The affine transformation function 512 for an image point is given by where (x, y) are the image point coordinates, and (x, y) 516 is their corresponding corrections, (ϕ, λ, h) is the 3-D coordi-517 nates of the image point on the ground, and (a 1 , a 2 , a 3 ) and 518 (b 1 , b 2 , b 3 ) are the six affine transformation parameters for 519 automatic extraction of tie and control points between adjacent 520 images. For each HiRIC image, first, the Harris operator [13] 521 is used to extract evenly distributed feature points; second, 522 to effectively eliminate the influence of image geometric 523 distortion caused by factors, such as terrain fluctuation, orbital 524 position, and attitude errors on the matching results, in the 525 automatic matching process of tie and control points, extract 526 a feature point P on the reference image or existing geographic 527 information data, define a matching image window W with P 528 as the center, project the window onto a small surface element 529 of the corresponding matching window W on the image to be 530 matched based on the elevation information interpolated by the 531 DEM and the imaging model, and optimize the affine trans-532 formation parameters based on holistic matching of window 533 W and W , and finally, the image near the predicted point in 534 window W is resampled and matched with the image window 535 W in the reference image or existing geographic information 536 data; then, the image coordinates of the points to be matched 537 is inversed by matching results. We call this process Image-538 Reshaping (seen in Fig. 3). After the image matching is 539 completed, the random sample consensus (RANSAC) algo-540 rithm [9] is used to remove the wrong matching points in the 541 automatic matching to further improve the accuracy of image 542 registration. [as shown in formula (2)], the block adju as follows: In the formula, B TP and B GCP represent the coefficient

563
Based on the least squares principle, the corresponding 564 normal equation matrix is For the characteristics of HiRIC imaging geometry, 568 we choose an easy-to-parallelize and efficient uncontrolled 569 combined block adjustment method by combining the well-  2) Assuming that the ground coordinates (ϕ, λ, h) of all tie 579 points are known, the affine transformation parameters (a 0 , a 1 , 580 a 2 ) and (b 0 , b 1 , b 2 ) of each image can be calculated, which 581 is called "backward intersection," and the normal equation of 582 this step is only a 6 × 6 matrix for each image. 583 3) Repeat 1) and 2) until the iterative change of all 584 unknowns is less than a predefined threshold to obtain the 585 optimal initial values of all unknowns. 586 4) Using the image affine transformation parameters 587 and the tie point ground coordinates, the image residuals 588 (v x , v y ) i,i=1,2,...,n of all tie points can be calculated, where n 589 means that the tie point appears in the n-view images. Let i=1,2,...,n ), and set a 591 priori weight of p = 1.0/v 2 max for each tie point. 592 5) Based on the ground coordinates of the tie points and 593 the corresponding a priori weights, the global least-squares 594 adjustment for RFM parameters is executed based on (1)-(4), 595 and the small-scale gross errors that still exist are removed 596 by the iterative weight of tie points calculated based on their 597 image residuals after adjustment.

598
In the uncontrolled adjustment method proposed in this 599 article, an "average" virtual control point is constructed by 600 the "ADMM" to solve the ill-conditioned (rank deficient) 601 problem of the normal equation matrix encountered in the 602 uncontrolled adjustment. On the one hand, this method can use 603 parallel computation to quickly calculate the initial values of 604 the unknowns to be adjusted, and on the other hand, it also uses 605 the "good" initial values to finally perform the least squares 606 overall adjustment. A dense, accurate, and reliable high-precision dense point 610 cloud matching algorithm for linear-array images is the key 611 to generating DEM data. In this article, a multi-image cor-612 relation matching algorithm (geometrically constrained cross 613 correlation (GC3) [47] is proposed based on 3-D geometric 614 constraints for HiRIC linear-array images on the basis of 615 the traditional image correlation matching algorithm. The 616 algorithm can organically combine the matching results of 617 various matching elements (feature points, feature lines, and 618 grid points) and make full use of the local and global texture 619 information of the image, so that the automatically matched 620 feature points and lines can fully express the important features 621 of the topography, which can be applied to the dense grid point 622 matching in the area without texture or with little texture to 623 meet the requirements of Mars image data processing.

624
In the dense point cloud matching process, one of the several 625 images is selected as the reference image, and the others 626 are all images to be matched. Using the Image-Reshaping 627 introduced earlier, we can obtain the matching results of the 628 reference image feature points on the search image pixel 629 by pixel using the precise orientation parameters obtained 630 after block adjustment. Similarly, distortions caused by the 631 terrain undulation, imaging geometry, and imaging scale can 632 be automatically compensated during the image matching 633 process. as the elevation data. The DOMs include two types of products 643 with a GSD of 0.7 and 3.5 m, respectively, and the GSD of the 644 DEM product is 3.5 m. The DOM product with 0.7-m GSD 645 has a total of 15 subframe products, and each subframe covers 646 a range of 30 × 20 km, with an overlap of ten pixels between 647 adjacent subframes (as shown in Fig. 4).       For the planar position precision, 36 evenly distributed 706 checkpoints were selected from the "Tianwen-1" landing area 707 DOM and compared with the planar position at the corre-708 sponding checkpoints in CTX mosaic. The results show that 709 the average planar position deviation between the "Tianwen-1" 710 landing area DOM and CTX mosaic is 24 m, with an rms value 711 of 15 m. It can be noted that a deviation of 15 m is much larger 712 than the control points standard deviations. The reason is that, 713 on the one hand, there would be a certain matching error, 714 because check points were manually extracted; on the other 715 hand, only the feature points with small horizontal deviation 716 are retained as control in the adjustment process. The position 717 precision of the control points reflects the internal consistency 718 of the topographic dataset, while the position precision of the 719 check points reflects the actual position deviation between 720 the topographic dataset and the existing terrain data. The 721 precision analysis shows that although the CTX dataset is an 722 uncontrolled mosaic, "Tianwen-1" landing area topographic 723 dataset registers well to the control source, which indicates that 724 our implementation works at the indicated level of precision 725 and accuracies on this order can reasonably be expected if a 726 more accurate control source is used.

727
For the elevation precision, the DEM of the landing area 728 was resampled to a spatial resolution of 200 m and subtracted 729 from the HRSC/MOLA Blended DEM to analyze the elevation 730 deviation between them. The results show that the average 731 elevation deviation between "Tianwen-1" landing area DEM 732 and the HRSC/MOLA Blended DEM is −7 m with an rms 733 value of 10 m, in which the deviation in 95.70% of the area 734 is less than 10 m (seen in Fig. 5). It should be noted that 735 the actual resolution of the HRSC/MOLA Blended DEM is 736 poorer than 200 m because of the large gaps between altimetry 737 profiles. Comparing the DOM of the landing area [seen in 738 Fig. 4(a)], it can be seen that the areas with larger deviation 739 are mainly concentrated in the bottom of craters, the boundary 740 of the suspected volcanoes, and so on, which is caused by the 741 large spatial resolution difference between the two types of 742 topographic data. The positional deviation statistics are shown 743 in Table V.

785
At present, the "Tianwen-1" orbiter is conducting the global 786 remote sensing exploration, and the main task is to achieve 787 global coverage of Martian surface using the moderate-788 resolution imaging camera (MoRIC), during which the HiRIC 789 mainly performs high-resolution imaging of typical terrain and 790 landforms on the Martian surface, such as valleys, volcanoes, 791 craters, and so on. Because HiRIC stereo coverage requires 792 frequent adjustment of the orbiter attitude, to ensure the global 793 remote sensing exploration mission, currently, no HiRIC stereo 794 images are acquired again for the time being. After the end of 795 the global remote sensing exploration mission, stereo imaging 796 will be arranged according to the exploration requirement 797 during the extended mission to supplement the stereo coverage 798 area of meter-scale high-resolution images on the Martian 799 surface. In addition, due to the acquisition of high-resolution 800 images of different seasons and different geographical loca-801 tions, HiRIC is also expected to provide basic data for the 802 solution of atmospheric correction, dust elimination, photo-803 metric correction, color correction, and other problems [5], 804 [18], [19], [32] in Mars high-resolution images.

805
In summary, the topographic dataset produced in this article 806 has achieved seamless mosaic and high-precision positioning 807 within the "Tianwen-1" landing area. It has been applied in 808 scientific and engineering missions, such as the geological 809 background study of the landing area, the positioning of the 810 "Zhurong" rover landing site, and the Martian surface remote 811 operation planning of the rover.