The size–frequency distribution (SFD) of rocks has been widely investigated on the previous Mars exploration missions (Golombek et al., 1997; Golombek and Rapp, 1997; Golombek et al., 2003; Ward et al., 2005; Golombek et al., 2006; Grant et al., 2006; Golombek et al., 2008; Grant et al., 2011; Craddock and Golombek, 2016; Golombek et al., 2021; Wu et al., 2021; Chen et al., 2022; Wu et al., 2022). Since rock distributions result from the geological process, including impacts, volcanic activities, aeolian processes, and even water transportation, the SFD of rocks could be an important target for the Martian geological history (Craddock and Golombek, 2016; Chen et al., 2022) and past climate analysis (Kereszturi, 2012). In addition, it plays an important role in quantifying the hazards of spacecraft landing and evaluating the risks for rover moving (Golombek et al., 2003; Golombek et al., 2021). Furthermore, it is significant for planning the potential infrastructure construction on Mars for the future.

The Tianwen-1 probe is China’s first autonomous Mars exploration mission, which was launched on 23 July 2020 and then entered the Mars orbit on 10 February 2021. After circling Mars for more than 3 months, the Zhurong rover, which was onboard the Tianwen-1 lander, finally successfully landed in the lowland area of southern Utopia Planitia on Mars on 15 May 2021. Among all the six scientific instruments mounted on the Zhurong rover, the Navigation and Terrain Cameras (NaTeCams) are used to take photographs of the Martian surface to provide support for the guidance, navigation, and control of the rover. In addition, these photographs could also be used for scientific observations. The rock SFD analysis in this paper is all based on the pictures captured using NaTeCams.

During the landing process, due to the impact of the retro-rocket engine on the Martian surface, an excavated depression was formed just below the lander. Compared with other depressions caused by the plume surface interaction (PSI), this depression is the deepest with an excavated depth of ∼ 40 cm (Xu et al., 2023), exposing the shallow depth of the soil here and the geological information beneath the surface. Since the high-resolution imaging camera (HiRIC) on the orbit and NaTeCams on the ground both provide images of the surface, the high-frequency and low-frequency radar detect the zones with a depth of several meters to tens of meters below the surface, the depth range of tens of centimeters has not been detected yet, and the investigation on this exposed depression just makes up for the research gaps in this depth. Most importantly, the structure of soil and rocks at the subsurface will definitely be more important for future drilling missions (Baker and Carter, 2019; Altieri et al., 2023), such as the future ExoMars rover with 2-m-deep subsurface access (Kereszturi et al., 2016). It is inferred that subsurface sampling will also be significant in the prospective Mars exploration, for example, according to the China National Space Administration, the collection and return of samples on Mars will be accomplished around 2030, and the preliminary research has already begun (Wei, et al., 2018; Hou et al., 2019). Thus, the geological structure at this depth can provide an important basis for the selection of future landing sites (Kereszturi et al., 2016) and drilling missions.

Apart from the excavated depression, the region around the lander was also disturbed during Zhurong touchdown, which was called the blast zone. It is known that the blast zone is brighter, having extremely different photometric properties from the surroundings, which may result from the smoothing of the surface and the redistribution of fine particles (Clegg et al., 2012). This means that small particles in the blast zone might be blown away by the plume, which affected the rock distribution here.

Before Zhurong went dormant, it had driven 1,921 m in total for the last 361 Sols and kept taking photographs along its route. These photographs enlarge our view of the southern Utopia Planitia and provide more information about the Martian surface, including the rock characteristics and distributions. Therefore, in this paper, the rock SFD is calculated and compared inside the excavated depression, within and out of the blast zone, and the influence of touchdown is also analyzed. Finally, we attempt to discuss the fractal dimension of rock sizes and its potential indication briefly and preliminarily.

Below the lander, there is a depression excavated by the retrorocket during touchdown. The NaTeCams mounted on the Zhurong rover are binocular stereo cameras, which can capture a pair of photographs of the observed object through their two camera lenses, CamB and CamA, from slightly different perspectives at the same time. Such a pair of photographs of the depression (shot on sol 12) is shown in Figures 1A, B. After the adjustment of contrast, the depression of the pictures is enlarged so that it can be seen more clearly (Figures 1C, D). It is worth noting that only one large columnar-like rock block with a height > 15 cm (Xu et al., 2023) stands separately in the deep center of the depression and that other small rocks are all embedded in the wall.

The calculation of the rock size and area includes four steps. The first step is the rock identification. Here, the regions of the excavated depression and inside rocks are outlined manually by brown and red lines (Figure 1C) for preparation. Second, based on the digital image correlation (DIC) algorithm (Blaber, 2023), the same depression region and rocks in the CamA image are matched automatically (Figure 1D). The size of the whole picture (Figures 1A, B) is 2048 × 2048 pixel. During automatic matching, a smaller area shown in Figure 1A, which is 300 × 300 pixel covering the depression, is selected. After tracking the features of this small area using the DIC algorithm, the matching area is shown in Figure 1B, and the displacement between them is calculated. According to this displacement of the small area, the displacement of each pixel within this area is also obtained, and the matching pixels in Figure 1B are found one by one. Based on the automatic matching of pixels, the third step is followed subsequently. According to the stereo vision 3D reconstruction method (Zhang et al., 2017; Zhang et al., 2020; Zhang et al., 2022), the 3D coordinates (X, Y, and Z) in the CamB coordinate system of each pixel with the corresponding root mean square error are calculated. To ensure the accuracy of the calculation, the pixel coordinates whose root mean square errors are larger than 1.0 mm were deleted and replaced based on the interpolation method. Furthermore, if the distance between the pixel and the CamB projection center along the optical axis (optical axis distance, represented by coordinate Z) is larger than 8.0 m, then this pixel was also excluded from the calculation. At last, the areas and sizes of rocks are computed. Since there is no notable evidence proving that the rocks are arranged in a specific direction, the observed apparent width is considered an average sample of the actual rock diameter (Golombek et al., 2021), which is the length in the horizontal direction of each rock. We chose this length as the rock size. According to the 3D coordinates, the areas of the three-dimensional surface of the outlined rocks and depression are also calculated. Thus, an important parameter reflecting the richness of rocks—rock abundance—is obtained, which equals the ratio of total rock areas to the depression area.

The blast zone (Figure 2A) is the region around the lander that was disturbed during descent (Clegg et al., 2012). For Tianwen-1, it was approximately 88.21 m wide and 129.60 m long, according to a rough estimate based on the threshold segmentation of the blast zone image (Figure 2B), captured by the High Resolution Imaging Science Experiment (HiRISE). After the touchdown, Zhurong started its tour mainly toward the south on Mars. This journey began from the center of the blast zone (Site A, Figures 2A, 3A), passing through its wing (sites B and C, Figures 2A, 3B, C), and reached the edge (Site D, Figures 2A, 3D) on 13–15 June (sol 30–32). After that, the rover kept walking, passing through a sand dune (Site S, Figures 2A, 4A) and also observing some normal and typical scenes (Site N, Figures 2A, 4B). Affected by the touchdown, the characteristics of rocks at sites A, B, C, and D in the blast zone varied. For comparison, the rock characteristics in unaffected regions, for example, at sites S and N out of the blast zone, were also analyzed. The number of these image pairs captured using the NaTeCams of each site is shown in Supplementary Table S1.

The calculation of rock sizes and areas out of the depression also takes the similar four steps. The difference is that when identifying the rocks at the first step, the Cascade region-based convolutional neural network (Cascade R-CNN) (Cai and Vasconcelos, 2018) is used instead of outlining manually, given the large number of rocks in the pictures (Figures 3, 4). This is a multi-stage object detection algorithm. In this study, ResNet-50 is used as the backbone of the model, which is augmented with a feature pyramid network for robust feature extraction (Lin et al., 2017). Then, a region proposal network (RPN) and a differentiable region of interest (RoI) pooling layer (Girshick, 2015) are used to generate proposal feature maps. For training and validation purposes, a dataset comprising 4,007 images with annotated bounding boxes highlighting Martian rocks is curated. We used weights pretrained on ImageNet (Deng et al., 2009) as the initial model and fine-tuned on our datasets to obtain the final model. Based on this method (Wang et al., 2022), the rocks are identified automatically. To identify as many rocks with different diameters as possible, only the rocks in the lower half of the picture are used for recognition and calculation because the upper half is far away from the camera, and small-sized rocks are easy to be ignored. The identification results of rocks are shown in Supplementary Figures S1, S2. We checked the recognition results of the rocks in the five images by a human, and the accuracy turns out to be 94.6% on average. The picture size in Figures 3, 4 is 2048 × 2048 pixel. Among these six pictures, the longest optical axis distance from the center of the picture to the projection center of the CamB is 7.97 m in Figure 3A, and the shortest distance is 3.62 m in Figure 4A.

As the distance from the camera increases, the actual size represented by each pixel (S _p ) changes regularly, as well as the size of the smallest rock (S _r ) identified at each distance. For example, at the optical axis distance (D) of 3.0 m, S _p and S _r are 1.58 mm and 9.98 mm on average, respectively; while at the D value of 8.0 m, S _p and S _r both increase to 10.81 mm and 107.89 mm on average, respectively. This means that at the optical axis distance D, the rocks smaller than S _r will not be identified in this research. In the log–log plot, S _p and S _r both show linear fits with D (Equations 1, 2 and Figure 5): log S p = 1.963 l o g D   −   0.739 , (1) log S r = 2.427 l o g D   −   0.159 . (2)

When calculating the areas of rocks out of the depression, we follow the hypothesis: the undersurfaces of the rocks are considered to be circles, and their diameters are the observed apparent width, as described previously. In addition, the area of the whole region is also calculated as a flat plane, and the rock abundance equals the ratio of the area of all rocks here to the area of this region.

There are some inevitable uncertainties in the calculation of rock sizes and areas. To start with, because of the limit of the image resolution, the small and more distant rocks are ignored. Second, some small rocks may be sheltered by the bigger ones. In addition, some rocks are partly buried or hidden in the depression wall or under the surface, and only the sides of exposed rocks facing the camera can be seen, which certainly affected the size and area calculation. Furthermore, the hypothesis that the undersurfaces of rocks are round may lead to inaccuracies in the calculation of the rock areas. The rock sizes and areas in different sites are shown in Table 1.

Fractal geometry can well-describe the regularity of irregular things in nature quantitatively, for example, the fractal characteristic of rock sizes reflects its fragmentation mode, in other words, the physical mechanism of rock evolution (Mandelbrot, 1998; Badge et al., 2002; Tu et al., 2005). This may illustrate that different fractal dimensions represent different geological processes during the formation of the rock blocks.

The fractal dimension of rock sizes could be calculated as follows: if rock size l > L, where L is a specific value, then the cumulative number of rocks N(l > L) with the characteristic size L varies as (Turcotte, 1986; Mandelbrot, 1998; Badge et al., 2002; Tu et al., 2005) N   l   >   L   ∼   L − D , (5) where D is the fractal dimension of the size distribution of the rocks. Performing a logarithmic operation to both sides of Equation 5, it can be obtained: log N l   >   L   ∼   − D l o g L . (6)

According to Equation 6, the least square method is used for linear fitting based on the rock sizes computed previously; then, the fractal dimension of rocks (D) and the corresponding adjusted determination coefficient (R ² ) are obtained (Table 2). Figure 7 shows the cumulative rock number in the depression versus the rock size in a log–log plot. It is notable that we cannot use a single power-law exponent to characterize the particle distribution, for in that case, both the relatively smaller and larger rocks are excluded (Figure 7A). Therefore, we roughly divide the rocks at each site into small rocks and large rocks to calculate the fractal dimensions (Figure 7B). The dimension of small rocks is called D ₁ with an adjusted determination coefficient R ² ₁ and that of large rocks is called D ₂ and R ² ₂ . The size that distinguishes these two kinds of rocks with different sizes is called the demarcation size (S _d ). We develop a scheme to divide the rocks—that R ² ₁ and R ² ₂ for each place are set to be as high as possible (>0.92)—with the highest average of all possible attempts. The fractal dimensions are shown in Table 2. For the rocks inside the depression, D ₁ and D ₂ are 0.43 and 2.06, respectively. The difference in fractal dimensions may reflect different fragmentation processes (Charalambous, 2015), in other words, different geological processes during the formation of rocks. This means that in the depth of tens of centimeters of the Zhurong landing site, the small rocks (12.83–21.72 mm) may experience different geological processes from the large rocks (>21.72 mm).

Through the fractal dimensions of rock sizes, we can infer that different geological processes have been involved; however, it is difficult to confine what process it exactly is. On the Earth, based on the particle size distribution of soil, the depositional environments could be inferred, according to the judgment formula established by Sahu (1964). However, this work is built on a large number of statistical results of particle sizes and the corresponding known environments. Unfortunately, the depositional environments on Mars are still unclear for now, and the particle size distribution is confined to only several landing regions. Thus, it is extremely hard to determine the geological process just through Mars’s rock distribution. This requires much other geological research studies and evidence, which is very important work in the future.

The fitting curves of the other five regions out of the depression are shown in Supplementary Figures S3–S7, and the fractal dimensions are also presented in Table 2. A similar phenomenon is observed. Since at each distance from the camera, the rocks smaller than a certain size cannot be identified (Figure 5), such small rocks are excluded in the calculation of fractal dimension D ₁ , which makes this value not precise. Thus, only the dimension D ₂ is used to analyze. For the touchdown-disturbed Blast Zone Wing, D ₂ is 2.54, and in the center of it, where the rocks had a greater influence by the touchdown, D ₂ is as high as 2.92. For the undisturbed regions like sites S and N, D ₂ is 2.00 on average, and it is similar for rocks in the depression. However, although the Blast Zone Edge is also affected by the descent, considering that it has largely maintained the original distribution of rocks, D ₂ is closer to that of sites S and N.

Similarly, S _d shows the same regulation as D ₂ . For the Blast Zone Center and Wing, where the rock distribution is greatly influenced by the touchdown, S _d of small and large rocks is 44.64–52.63 mm; however, for sites S and N where the rock distribution is undisturbed, this size is 13.17–14.82 mm, which is a lot smaller. With regard to the Blast Zone Edge, S _d is also much closer to that of sites S and N, which is 16.83 mm. This means that the rocks in the Blast Zone Center smaller than 52.63 mm may be disturbed by the touchdown, as well as the rocks smaller than 44.64 mm in the Blast Zone Wing. In addition, the fractal dimension for large rocks at the disturbed regions is usually higher than that of the undisturbed regions.

In this paper, the rock size–frequency distributions are calculated in six regions of the Zhurong landing site, Mars, including regions inside the depression, within and out of the blast zone. Among them, the distribution of rocks in the depression is the first time to be obtained, before which the shallow subsurface on Mars with a depth of tens of centimeters has not been investigated. Apart from that, the fractal dimension is calculated to quantitatively describe the size distribution and reflect the geological process experienced by them in an attempt. The following conclusions are drawn:

1) The rocks in the depression and out of the blast zone are original, without disturbance by touchdown. The rock abundance in the depression is 4.15%, smaller than that of the surface out of the blast zone (the rock abundance at sites S and N is 7.66% and 6.95%, respectively).

For future Mars exploration missions, the analysis of the rock SFD in this research will provide important reference, especially for landing, drilling, and sampling. According to Golombek et al. (2021), rocks whose sizes are larger than 700–900 mm are considered to pose a potential threat to landing. Based on this research, on the surface we investigate, only one rock reaches a size of ∼400 mm; others are all < 250 mm. To be more conservative, even if we consider any rock larger than 400 mm to be a landing hazard, this place is still relatively safe. This could also be proved by the low density of larger rocks. The density of rocks with a size of > 100 mm is 0.5–1.3/m², and that for rocks > 200 mm is 0–0.3/m², which means the probability that the spacecraft will encounter a rock larger than 100 mm is relatively low. However, the excavated depression under the Tianwen-1 lander exposes a large rock, which might be in situ (Xu et al., 2023). It shows that the rock size of the shallow subsurface is not uniform, and it is likely that there are large blocks buried in it, which may cause deflection and obstruction of the drill bit during drilling and sampling.