N Tanabe

Y Cho

E Tatsumi

T Ebihara

K Yumoto

T Michikami

H Miyamoto

T Morota

C Honda

P Michel

K A Otto

O S Barnouin

K Yoshioka

H Sawada

Y Yokota

N Sakatani

M Hayakawa

R Honda

S Kameda

M Matsuoka

M Yamada

T Kouyama

H Suzuki

K Ogawa

S Sugita

Development of image texture analysis technique for boulder distribution measurements: Applications to asteroids Ryugu and Itokawa

Keywords:

Recent asteroid missions have revealed that ma ny sub-kilome ter asteroids are rubble piles. Large parts of their surfaces are covered with boulders larger than tens of centime ters. An evaluation of the abun danc e and size distribution of boulders provides clues to un derstand surface processes on boulder-covered asteroids. Here we report a new me thod that automa tically me asures the abun danc e of sm all boulders (sub -pixel to a few pixels), wh ose boun daries cannot be recognized with visual inspection, by quantifying the surface radiance variation that occu rs du ring the spinning of the asteroid. After validating our approach with previous boulder coun ting data, we apply this me thod to images of the asteroids Ryugu and Itokawa, wh ich were visited by JAXA's Hayabu sa and Hayabu sa2, and obtain a global distribution of the boulders larger than 0.75 -3 m, wh ich corresponds to 1.5-6 pixels. We find that the boulder nu mbe r density of this size range is sm aller (1 ) on the western bu lge than on the eastern hemisphere and (2 ) on the equatorial ridge than on the higher latitudes, both of wh ich exceed the nu mbe r density of boulders > 5 m by an order of ma gnitude. The boulder size distribution at 1.25 -20 m shows that the boulders sm aller than 1 m are more abun dant at the equator than at mid-latitudes, wh ile those larger than 1 m in diam eter are more abun dant at midlatitudes than at the equator. This contrast suggests size-depend ent migration of boulders in the latitudinal direction. We also find that the typical boulder size (the size reaching the cu mu lative areal coverage of 50 %) is 1.9 m at the equatorial region (1 0°S-10 °N) wh ile it is 2.6 m at mid-latitudes (4 0°S-50 °S, 40 °N-50 °N). The typical boulder size is also sm aller in the western bu lge (2 .0-2.2 m) . We construct global ma ps of the power-law index of the size frequency distribution of boulders and find minor variations over the entire surface of Ryugu (-2.53 ± 0.03 ) for boulders larger than 1.25 m. This sm all variation suggests homogeneous size sorting processes on Ryugu. Surface roughness does not show a significant correlation with the v-band albe do bu t shows a high anti-correlation (R = -0 .73) with the cu rrent geological slope on the eastern hemisphere. Our me thod is useful enhanc eme nt of sm ooth area detection and boulder distribution ch aracterization that will be applicable to other planetary explorations in the future, includ ing those of Pho-

bos and other asteroids.

Introduction

Recent near-Earth asteroid (NEA) missions have revealed that many sub-kilometer asteroids, such as Itokaw a, Ryugu, and Bennu, are rubble piles with high abundance of boulders (Fujiwa ra et al., 2006;Watanabe et al., 2019;Lauretta et al., 2019). These observations are consistent with theoretical predictions that many small asteroids would be formed from catastrophic disruption of larger parent bodies and subsequent re-accumulations of fragments (Walsh et al., 2018;[START_REF] Michel | Collision and gravitational reaccumulation: possible formation mechanism of the asteroid Itokawa[END_REF]Durda et al., 1998). These spacecraf t observations also revealed great details of these small asteroids (e.g., Miyamoto et al., 2007;Sugita et al., 2019). For example, average particle size on Itokaw a, Ryugu and Bennu, are several tens of meters to several meters in diameter (Saito et al., 2007., Sugita et al., 2019;Michikami et al., 2019;Walsh et al., 2019). Thus, analysis of such sizes of boulders is useful for understanding the surfa ce processes and properties, such as mass movement, resurfa cing, spectral change, and thermal inertia.

The size and spatial distributions of boulders were investigated by a number of previous studies by visual counting (e.g., Dombard et al., 2010;Pajola et al., 2015;Michikami et al., 2008). These analyses focused on the spatial density of boulders (Dombard et al., 2010), sizefrequency distributions (Michikami et al., 2008), and the shape and orientation of boulders (Nakamura et al., 2008;[START_REF] Jawin | Global pa tterns of recent ma ss mo vement on as tero id (101955) Bennu[END_REF]. When multiband images were available, the color of boulders wa s obtained to discuss their compositions, degree of space weathering and solar heating (Sugita et al., 2019;Walsh et al., 2019). On Itokaw a, spatially heterogeneous boulder number density is thought to have been caused by granular motion, which resulted from repeated seismic wa ves generated by impact events (Miyamoto et al., 2007). On Ryugu, local highresolution images show an evidence of regolith migration from equatorial ridge towa rds the current low potential region, such as imbricated boulders and run-up regolith on larger boulders (Sugita et al., 2019).

The size-frequency distribution (SFD) of boulders on planetary bodies is used to investigate the mass flow of surfa ce materials, size sorting effect on the surfa ce, and formation process of regolith and boulders. Assuming that the boulders are fragments of larger disintegrated boulders, the SFD of boulders follows a power law. The fragmentation of boulders can be described as a fracture process by collisional cascade, in which self-similar fragments are produced [START_REF] Turcotte | Fr actals an d Chao s in Geol og y an d Geophysi cs (Cup)[END_REF]. The power-law index varies depending on many reasons, such as different mechanical strengths of boulders (Ta naka et al., 1996) or the fragmentation process of boulders (Hartmann, 1969). Surfa ce processes that sort the size of boulders (e.g., mass flow, convection, displacement by impacts (Mottola et al., 2015;Miyamoto et al., 2007;Da ly and Schultz, 2016)) subsequently af fect the index. The actual power-law index of the SFD of boulders ranges from -2 to -4 on planets, moons, asteroids and comets (e.g., [START_REF] Turcotte | Fr actals an d Chao s in Geol og y an d Geophysi cs (Cup)[END_REF]Thomas et al., 2001). On Itokaw a, for example, Michikami et al. (2008) reported regional variations in the power law index: 3.1 ± 0.1 (global average), -2.8 ± 0.1 (head region), -3.2 ± 0.1 (body region), suggesting a material difference between the two bodies that constitute the head and body of Itokaw a, which may have been originally separated and then eventually merged together. This two-body merge scenario wa s also supported by the density contrast between each body (Lowry et al., 2014;Kanamaru and Sasaki, 2019). Ryugu exhibits an abundance of boulders larger than 20 m and regional variation of boulder number density between the eastern and western hemispheres; furthermore, lower latitudes host fewer boulders (Michikami et al., 2019). The power-law index of boulder distribution is smaller (gentler) than that of Itokaw a. This difference wa s interpreted as the result of the burial of small boulders by even smaller boulders, because large boulders (~10 m) were buried by small boulders (~1 m) (Michikami et al., 2019). However, except for some high-resolution images acquired locally during Hayabusa2's touchdown operations, the boulder distribution wa s measured for the size range larger than 5 m, while the global coverage of such boulders wa s only 20-35%. The boulders larger than 1 m cover 80-90%. The distribution of boulders including such smaller boulders (sub-meter to a few meters) on Ryugu wa s not investigated globally because of the insufficient image resolutions.

Because boulders smaller than ~5 pixels cannot be recognized in a reliable wa y (Murdoch et al., 2015), "smoother" or "rougher" textures found on images wa s used to infer the surfa ce roughness or the abundance of small boulders (sub-pixel to a few pixels) (Fig. 1). However, there is no method to convert the qualitative textures of images to actual boulder number densities. Moreover, counting hundreds of thousands of boulders on the entire asteroid surfa ce requires an enormous amount of time and effort (e.g., Michikami et al., 2019for Ryugu, Walsh et al., 2019 for Bennu). This poses a problem when, for example, a quick evaluation/selection of the safest landing site from the entire asteroid surfa ce is needed (Kikuchi et al., 2020). For such lander missions, including recent and future sample-return missions, analyzing the distribution of sub-meter boulders is important for ensuring spacecraf t safety and efficient sample collection (see Fig. 2).

Thus, in this study, we develop a method that estimates the abundance of boulders smaller than the resolution limit (~5 pixels) from the texture of the surfa ce. Our method calculates the standard deviation of radiance of a target area using images from different viewing angles. We apply this novel technique to the images of Ryugu and Itokaw a to evaluate the abundance of boulders, which were indistinguishable with the visual counting approach. We obtain maps of boulders coverage of 0.75-40 m, which cover 80-90% surfa ce of Ryugu, and a map of power-law index of SFDs to investigate the degree of local size-sorting on this asteroid.

Image analysis method

In this section, we outline a method to estimate the abundance of small boulders based on the qualitative texture of the images.

Measurement principle

We quantify the texture, or the "smoothness of an area", using the standard deviation of radiance in multiple images of a given target area. Because Hayabusa2 observed Ryugu while hovering above the surfa ce, the solar incidence angle and emission angle to/from a certain area changed as the asteroid spins. Within this region, the inclination and shape of individual particles or boulders generate an observable radiance variation due to the varying reflection from individual particles/boulders. For example, on a surfa ce that hosts boulders significantly larger than pixel scales, the asteroid's spin changes its radiance significantly because illumination and observation conditions, such as incidence, emission, and phase angles, change according to the direction and angle of individual boulder fa cets. In contrast, if a number of small particles cover the area of interest (i.e., the area appears smooth in an image), the inter-image deviation of radiance within the area is much smaller than that of a rougher area with larger particles. This is because the more particles (regolith, pebbles and boulders) of different orientations and shapes cover, the more homogeneous is the photometry of the pixel as a whole. We did not take the morphology or surfa ce texture of individual boulders into account in this automated process. Although rougher boulders make σ sp higher, the effect of such surfa ce features on σ sp is bracketed in the error of the boulder coverage.

Below, we show that the standard deviation of radiance of the same region observed with at different geometries indicates the particle size on the target region. To formulate this relation, we define σ sp as the standard deviation of radiance of each pixel due to the spin of the asteroid. The value of σ sp qualitatively shows the size of the boulders and will be correlated with actual boulder size and its coverage in section 2. 4.

( : number of images with different observation angles, : radiance of a pixel for each viewing angle, : average of radiance of the pixel of interest, observed in N images, σ sp : standard deviation of radiance of the pixel of interest).

First, we measured the radiance of individual pixels in an image using asteroid Ryugu L2b image data. The standard deviation of the radiance over different images with different viewing angles was obtained by the tracking of each pixel. For example, if a pixel on a specific boulder is in the area of interest, the specific pixel appears at different pixels in different images. The location of the pixel on the boulder is tracked over different images by assigning each pixel in the image with a polygon ID of the shape model (SHAPE_SFM_3M_v20180804; Watanabe et al., 2019). The radiance of the pixel was measured in each image to derive the standard deviation of the radiance over multiple images. For this corregistration process, we used the information on the position and attitude of the spacecraft at the time of image acquisition. Then the image of the shape model was reproduced from the viewing angle exactly the same as the real ONC-T images. The difference in the asteroid rim position was less than 1 pixel between the real image and simulated image at low-to mid-latitudes (45°S -45°N). The shape model we used had relatively large uncertainties in high-latitudes because of the lack of the images used to create the model (Watanabe et al., 2019). Thus, we calibrated the coverage of boulders between 45°S and 45°N. Our results on the even higher latitudes (45-60°) were extrapolated based on the calibration established at < 45°, although the position of asteroid's rim in the simulated images matched those of the actual images within 1 pixel on average. The boulder number density measurements at the high latitudes require further study with an updated shape model. The size of the facet of the shape model was comparable or smaller than each pixel in the image (~0.6 m/pix for Ryugu, ~0.4 m/pix for Itokawa) because the shape model consists of 3 million triangular facets with a side of ~0.4 m. If one pixel in the image covered multiple facetss on the shape model, the image pixel was assigned to the facet with the center closest to the pixel center. Thus, every pixel in an asteroid image was assigned to a specific facet on the shape model. The image pixels assigned to the same polygon (facet) ID were regarded as observing the same location of the asteroid.

We calculated the incidence/emission angles of each facet from the shape model. The number of images N to obtain the standard deviation of the radiance of each pixel was 5. Note that is irrelevant to absolute boulder size at this stage because our approach relies only on pixel scales. The size of boulders of interest can be changed by changing the resolution of image sets by binning. This feature is used in Section 3.2 to derive the power-law index of the boulder sizefrequency distribution. From the power-law index and size-frequency distribution, we calculate the typical boulder size (the size reaching the cumulative areal coverage of 50%) for each region.

Dataset

We used the images acquired during multiple observation campaigns at different altitudes: so-called home position (HP, 20 km above the surfa ce, 2.0 m/pix, Ryugu spins 3°per image; 127 images obtained on July 10th, 2018) and mid-altitude covering 50°S-50°N (5 km, 0.5 m/pix; 91 images obtained on August 1st, 2018). We conducted the same analysis using 60 images of Itokaw a (ST_2420699984-ST_2426434412), which were obtained while the asteroid rotated at an interval of 3°between subsequent images. The resolution of the Itokaw a images is 0.4 m/pix. We acquired these images from JAXA's Da ta Archives and Transmission System (DARTS). The shape models of Itokaw a (Stereo Photoclinometry (SPC) model version 5.04) and Ryugu (SHAPE_SFM_3M_v20180804) were used to calculate incidence/emission angles conditions for the images and to derive geologic slopes of the asteroid surfa ce.

Model caveats

We note two caveats for potential bias in due to the photometric effects caused by different incidence/emission angles and shadows of boulders. First, different degrees of variations in the incidence/ emission angles result in wider radiance variations of pixels in low latitudes than in high latitudes. Because our method does not require photometric correction, the different degrees of the photometric effects lead to an overestimation in boulder number density at lower latitudes: the equatorial regions tend to show larger variation in pixel radiance because the variation in the photometric effects tends to be larger around the equator (Fig. 3a). Because of the effects, widely different incidence/emission angles yield a large variation in apparent pixel radiance among the images. Without any corrections, the degree of overestimation would be 0.5% in regardless of the abundance of boulders, when the average of was 5-10%. To reduce the radiance deviation due to the photometric effect, we made a criterion on incidence and emission angles toward and from individual facets: and ( : incidence angle, : emission angle). Based on this criterion we removed the deviation by extreme photometric conditions: 7% of facets at 45°S-45°N were excluded. The facets that did not meet the criteria were not used for our σ sp calculations. With this correction, the apparent radiance variation due to the photometric effect is within the level of the calibration accuracy (section 2.5). The thresholds of the angles were set to make the number of radiance data larger than 1000. If the conditions of the incidence/ emission angles are too strict, the number of data points (facets) becomes very small. For example, if then the number of pixels at 0-10°N falls to only ~100, hindering the assessment of global boulder number distribution.

The second caveat with our method is involved in the processing of shadows cast by boulders. Without any corrections, the boulder number density would be estimated to be lower than the actual value in areas where numerous boulders cast many shadows. This is because our algorithm removes the shadow associated with boulders from calculations. We excluded facets showing a radiance of less than 1% of the global average. Because numerous boulders cast shadows as the asteroid spins on a boulder-rich surface, one boulder on the surface could be in the shadow of nearby boulders in another image. As a result, relatively boulder-poor areas in the boulder-rich region were selectively used for our calculation. This procedure led to a systematic underestimation of boulder coverage by 1.0-1.5% in the northern hemisphere, where boulder coverage is 20-25% (boulders of 5-20 m). Nevertheless, such underestimation would result in 1.0-1.5% in coverage, which is within the level of our calibration ac-curacy (section 2.5). With this procedure, ~6% of the area in 25-45°N and ~2% of 25-45°S were excluded from the calculation. The north-south contrast occurred because the northern hemisphere was in winter in our dataset and tended to show shadows that are more abundant.

Calibration of

In this section, we calibrate , the standard deviation of the radiance among five images, with the boulder number density obtained by visual counting using higher-resolution images. First, we calculated globally with the low-resolution image sets obtained from the spacecraf t home position (pixel resolution ~2 m/pix). We then correlated the coverage of the pixels showing σ sp higher than a threshold with the coverage of boulders (%) obtained from visual boulder-counting data of selected regions (Fig. 4).

For this calibration, we counted boulders in higher-resolution images acquired from the mid-altitude operation (~0.5 m/pix). We used the images of the 7 candidate sites (100 m × 100 m) that were selected for Hayabusa2's first touchdown (Kikuchi et al., 2020). The smallest boulders counted were 3 m (= 6 pix) in diameter. Because few or no boulders larger than 12 m were found in the regions we used for calibration, we calculated the coverage of boulders of 3-12 m, which would occupy 1.5-6 pixels in the home-position images used for making the σ sp maps. The coverage of the 3-12 m boulders in the seven 100 × 100 m areas were between 12 and 27, providing the "ground truth" for the calibration. Then, we constructed σ sp maps in these regions using the lower-resolution images taken from the HP. Each region wa s expressed by 50 × 50 pixel because the pixel resolution wa s 2 m/ pix. Note that boulders smaller than 12 m were no longer resolvable with these HP images.

To measure the cumulative boulder number density in a certain area (km -2 ), which is usually measured using the visual counting approach, we calculated the ratio of pixels that showed a σ sp larger than a threshold (σ 0 ), to the total pixel number in the areas. The area coverage of σ sp > σ 0 in the 50 × 50 pixel area wa s calculated in each area and correlated with the boulder coverage of 1.5-6 pixel (0.75-3 m) measured with our boulder-counting analyses (Fig. 4). We found that the σ 0 = 0.1 yields the highest correlation between the coverage of pixels showing σ sp > σ 0 and the actual coverage of boulders (Fig. 5). The threshold could differ for different cameras or optical setups, therefore this calibration should be performed when a different instrument is used. The fitting line between "the coverage of the region exhibiting σ sp > 10%" and "actual coverage of boulders (3-12 m)" wa s used as the calibration line to convert any σ sp distribution to the coverage of boul- 5. Concept of our method to convert σ sp of individual pixels to the coverage of boulders. We calculated σ sp from images for each region. We binarized the σ sp map with different thresholds σ 0 and calculated the coverage of pixels that show σ sp higher than σ 0 . The pixel coverage is compared with actual boulder coverage of the region. The same images as Fig. 4 were used. The figures on the right show the binary maps and their context. ders (Fig. 6). The size of the area in which we calculate coverage can vary depending on the objective of the analysis. The root mean square of the coverage prediction wa s 1.85% when boulder coverages were 10-30% (Fig. 6a). This value wa s used in this study as an intrinsic uncertainty in our measurements.

Validation

To validate our method, we derived the boulder coverage with our calibration line and compared the results with the visual boulder counting data reported in a previous study (Michikami et al., 2019). We calculated σ sp that corresponds to the boulders of 5-20 m to match the size range reported by Michikami et al. (2019). To conduct this analysis, we changed the effective resolution of the images of Ryugu from 2.0 m/pix to 3.3 m/pix by binning the images with 0.5 m/pix resolution. We found that the coverage values obtained with our calibration line were consistent with the visual counting results (Fig. 7) within the range of 2σ. The RMS error in coverage wa s 1.5% (7.9% relative error when normalized by the median coverage of 19%) when compared to the visual counting data. The amount of error in our analyses wa s comparable to the intrinsic error in visual counting. This amount of error includes those from systematic overestimation by photometric conditions and shadows discussed in Section 2.1. This result demonstrates that our method can derive the coverage of boulders (>1.5 pix) from the texture of the surfa ce, including the boulders that are too small to be visually distinguished. The high correlation between the coverage of boulders analyzed by our method and those analyzed by visual counting also demonstrates the capability of changing the size range of boulders by artificially changing the resolution of images (i.e., binning). Such a spatial resolution control will be used to derive the power-law exponent of the boulder SFD in section 3.2.

Boulder distribution on ryugu

In this section, we show the global distribution of boulders (0.75-3 m in diameter) on Ryugu, which were too small to investigate globally in previous studies. We produced global maps for boulder number densities, derived the distribution in the longitudinal and latitudinal directions of boulders larger than 0.75 m, and compared the result with previous studies. We then construct a global map of the power-law index of boulder SFDs by producing lower-resolution images (~1.7, 3.3, 6.7 m/pix) and repeating our method on these images. We compared these parameters with the v-band albedo and surfa ce slope obtained in previous studies.

Distribution of boulders in diff erent longitudes and latitudes

Fig. 8a shows the global distribution of small boulders (0.75-3 m), which wa s not mapped in previous studies because of insufficient image resolutions. We investigated the features of boulder distribution in the longitudinal and latitudinal directions at a boulder scale smaller than that visually counted in previous studies (Figs. 6 and7). First, the longitudinal distribution of the boulder number density indicates that the coverage of boulders on the eastern hemisphere is larger than that on the western bulge (Fig. 9a). Our study reveals that the dichotomy in boulder number density observed at > 5 m (Michikami et al., 2019) extends to sub-meter boulders (Fig. 9b). The area of smallest boulder coverage (0.75-3 m) is found on the ~100 m × 100 m area on the trough (Fig. 8c).

Second, the latitudinal distribution of boulders of 0.75-3 m shows that fewer boulders are present at lower latitudes (20°S-20°N) (Fig. 10a), extending the trend found for boulders (>5 m) (Fig. 10b). This result is consistent with the geomorphological evidence that boulders at low latitudes are covered with finer boulders (Michikami et al., 2019). One interpretation to account for this distribution is the different mobility due to the particle size (e.g., Phillips et al., 2006).

Power-law index of boulder size-frequency distribution

In this section, we obtain the power-law index α for boulder SFDs globally in the size range 10 times smaller than those obtained studied in previous studies. We then produce α maps of Ryugu and Itokaw a, to investigate the size distribution of boulders on these asteroids.

Distribution on ryugu

We investigated the spatial distribution of the power-law index and the size distribution of boulders on Ryugu. By 2 × 2 or 4 × 4 binning of the images obtained at the altitude of 5 km (0.5 m/pix), we reduced the effective pixel resolution and obtained the coverage of boulders corresponding to 1.5-6 pix in each image. The power-law index was calculated for each facet with the four cumulative coverage data at 1.25, 2.5, 5.0 and 10 m (Fig. 11).

The global map of the power-law index shows a very small variation in α over all of Ryugu's surface (Fig. 12). Note that the map tends to show the smaller features of α in the high latitudes because the area specified by 5 °× 5 °is almost half of that in the equatorial region. However, this latitudinal dependence of area size does not affect the derivation of σ sp or α because σ sp is calculated at a much finer scale (individual facets ~0.4 m). The average and standard deviation of the index were -2.53 and 0.03, respectively. This small variation indicates that the lateral boulder size-sorting process is either weak or very homogeneous on Ryugu, producing no particularly smooth regions and other rough regions as observed on Itokawa. We also found that the power-law index decreases (gentler) with decreasing boulder size globally. This global trend is consistent with the result of visual counting in a few local areas, where boulders of several to tens of centimeters were measured with high-resolution images (Michikami et al., 2019). The power-law index for the boulders of 1.25-10 m is consistent with that reported for the boulders larger than 5 m (Michikami et al., 2019) equal to -2.53 ± 0.05, but shallower than the result inferred by the data for those larger than 20 m (Sugita et al., 2019) because of the different size ranges (Table 1).

The latitude-dependent SFDs using the boulders of 1.25-20 m show that boulders smaller than 1 m are more abundant at the equator than in the mid-latitudes, while those larger than 1 m are more abundant in the mid-latitudes than on the equator (Fig. 11). This different crossover in boulder SFD depending on latitude could be because of the size-

Table 1

Power-law indices of boulder size-frequenc y distributions in low latitudes, mid latitudes, eastern hemisphere (0 -13 5°and 33 0-36 0°), and western bu lge (1 80 -30 0°) (Sugita et al.,20 19 ). The longitudes 13 5-18 0°and 30 0-33 0°are wh ere the fossae are located and not categorized either east or west. dependent mobility of boulders which may fa vor the travel of larger distances of smaller boulders, although further study is required to assess how such a size-dependent mobility works on small asteroids. One possible scenario is that the boulders smaller than 1 m flowed towa rds the equator in the paleo high-spin-rate era of Ryugu, while larger (5-10 m) ones were left at mid-latitudes; the <1 m boulders are still on the equatorial ridge because the current geopotential gradient is not steep enough for the small boulders to flow down to the mid-latitudes.

Another hypothesis to explain this latitudinal distribution would be that the small boulders are buried by the larger ones in the mid latitudes due to regolith migration processes, such as the Brazil nut effect (Miyamoto et al., 2007). These smaller particles would have been deposited before/during the paleo high-spin-rate or broken down in place.

Distribution on itok awa

We applied our method to the images of Itokaw a (0.4 m/pix) and obtained the roughness (σ sp ) map and power-law index map (Fig. 13). The lowest σ sp occurs in the smooth regions named Muses Sea and Sagamihara. The averages of the power-law index for global, head, body, west and east regions were -3.6, -2.8, -2.9 and -3.2, respec-tively. These values are consistent with those obtained by previous studies (Michikami et al., 2008 (global: 3.1 ± 0.11, head: 2.8 ± 0.1, body: 3.2 ± 0.1, east: 2.9 ± 0.1, west: 3.2 ± 0.1), Ma zrouei et al., 2014 (head: 3.1 ± 0.4, body: 3.6 ± 0.2), Ta ncredi et al., 2015 (head: 2.96 ± 0.36, -3.6 ± 0.2)), further validating our method based on radiance variation. The standard deviation of calculated for the entire surfa ce of Itokaw a is 0.83, which is much larger than that of Ryugu. This large variation reflects the boulder distribution contrast in smooth regions with few boulders (Muses Sea, Sagamihara) and other rougher surfa ces on Itokaw a. This difference could be due to granular motion, which is thought to have transported small particles to the lowpotential regions (Miyamoto et al., 2007). In contrast, our results suggest that such an effect is not prominent on Ryugu possibly because of different particle sizes and densities, cohesive forces, the symmetric shape, the large size or the high porosity of Ryugu.

Typical boulder size on ryugu

The typical boulder size (i.e., the size at which the cumulative coverage reaches 50%) wa s calculated from the SFDs of different longitudinal/latitudinal bins. We found that the typical boulder size is 1.9 m at the equatorial region (10°S-10°N) and 2.6 m at mid-latitudes (40°S-50°S, 40°N-50°N) (Fig. 14a). We also found an east-west contrast in typical boulder size (Fig. 14b).

Comparison with v-band albedo and topographic slope

Because our analysis yields the global maps of small-scale roughness and power-law index of boulder size distribution, they can be compared with different observables that characterize the asteroidal surfa ces (Fig. 15). For example, the v-band albedo shows little correlation with σ sp (R = 0.34), while the albedo shows a higher correlation with the power-law index (R = 0.58) (Fig. 16). The reason for the higher correlation between the power-law index and the v-band albedo is unknown, but one hypothesis would be that size-sorting processes on Ryugu, which would change the values of α, made the variation of vband albedo. Second, σ sp and the current geological slope exhibit a high anti-correlation (R = -0.73) on the eastern hemisphere (latitudes 0-45°) (Fig. 13), while the western bulge shows a weaker correlation (R = -0.49). The substantial negative correlation in the eastern hemisphere indicates that fewer boulders occur on slopes and more boulders sit on the plain. This distribution is consistent with the interpretation that the boulders migrated to the current gravity low, suggesting that the relocation of boulders took place af ter the spin state slowed down to the current level (Watanabe et al., 2019). The reason for the low correlation in the western bulge is unclear, but it might be because the western bulge is younger than the eastern hemisphere (Cho et al. in revision) and boulders are currently gravitationally unstable. This trend suggests that the distribution of boulders on the equatorial region slid down along with the current geopotential, not with that of the paleo high-spin-rate era, and it created the mass flow features reported by Sugita et al. (2019).

Conclusions

We developed a novel method to estimate the abundance of boulders of 1.5-6 pixels by measuring the standard deviation of radiance (σ sp ) that individual pixels exhibit as an asteroid spins. Then we calibrated the coverage of regions having σ sp > 0.1 with the coverage of boulders. This calibration yielded a boulder coverage consistent with those obtained in a previous study of Ryugu, which wa s conducted by visual counting (Michikami et al., 2019). From a comparison with this published boulder coverage, the root mean square error of our prediction wa s determined to be 1.85%, when the median boulder coverage wa s 30%. This level of uncertainty is comparable to the intrinsic uncertainty of visual boulder counting.

Using this method, we obtained the global distribution of boulders of 0.75-3 m (1.5-6 pixels), which is otherwise not possible to extract visually due to the resolution limit. We found that the boulder number density is smaller (1) on the western bulge than on the eastern hemisphere and (2) on the equatorial ridge than on the higher latitudes, both of which extend the global trend reported for boulders > 5 m (Sugita et al., 2019;Michikami et al., 2019) to a sub-meter scale. The latitudedependent SFDs using the boulders of 1.25-20 m in size show that boulders smaller than 1 m are more abundant at the equator than at midlatitudes, while those larger than 1 m are more abundant at midlatitudes than at the equator (Fig. 11). This different crossover in boulder SFD depending on latitudes could be because of the size-dependent mobility of boulders.

The SFD of boulders of 1.25-10 m wa s derived globally by changing the pixel resolution through binning the original images. The global map of the power-law index shows a very small variation in α over the entire Ryugu (Fig. 12). The small variation in α (-2.53 ± 0.03) suggests that the lateral boulder size-sorting process is either weak or very homogeneous on Ryugu, producing no particularly smooth regions and other rough regions as observed on Itokaw a. The power-law index for the boulders of 1.25-10 m is consistent with that reported for the boulders larger than 5 m on Ryugu (Michikami et al., 2019).

Using the SFDs of different longitudinal/latitudinal bins, we found that the typical boulder size (i.e., the cumulative coverage reaches 50%) is 1.9 m at the equatorial region (10°S-10°N) and 2.6 m at midlatitudes (40°S-50°S, 40°N-50°N). We also found an east-west contrast in typical boulder sizes.

The application of our technique to Itokaw a images demonstrated that smooth areas, such as Muses Sea and Sagamihara, can be detected by our approach. The power-law index of the SFD showed a remarkable difference between the smooth areas and other rougher areas.

The roughness index σ sp and power-law index α were compared with the v-band albedo and geologic slopes on Ryugu. We found that the vband albedo shows little correlation with σ sp (R ~0.34), but shows a higher correlation with the power-law index (R ~0.58). Although a different power-law index suggests that this is the result of size sorting (e.g., very steep α in the smooth region on Itokaw a), the reason for the higher correlation between the power-law index and v-band albedo is unknown: it is possible that size-sorting processes af fect the v-band albedo. The surfa ce roughness σ sp and current geological slope show a high anti-correlation (R ~-0.73) on the eastern hemisphere. This anti- correlation is consistent with the interpretation that the relocation of boulders occurred af ter the spin rate of Ryugu slowed down to the current speed.

The results of this study demonstrated that the newly developed method of boulder analysis without extensive visual counting can be used for the smooth area detection and boulder distribution characterization. This approach would be useful for other planetary explorations in the future, including other asteroids as well as Ma rs' moon Phobos and other asteroids.

Fig. 1 .

 1 Fig. 1. Image of Itokawa taken by Hayabu sa. Note the sm ooth surface texture at Mu ses Sea and the surroun ding rough surface texture. (Image ID: ST_249 30 31 59 4) .

Fig. 2 .

 2 Fig.2. Concept of our me thod. We use images of the same region taken from different angles wh ile the asteroid spins. The position of a region of interest occu pies different pixels in each image, bu t the pixels corresponding to the region are tracked and the standard deviation of their radiance is me asured.

Fig. 3 .

 3 Fig. 3. Number of pixels showing individual incidence/emission angle conditions for the images of (a) low latitudes (5°S-5°N) and (b) higher latitudes (15°S-25°N); (c) reflectance of Ryugu's surface of calculated by a photometric model (Tatsumi et al., 2020). The region we use is surrounded by black facets. Note the distribution contrast in the photometric conditions between the different latitudes. The combination of the angles (i.e., photometric conditions) occupy different area between low-and mid-latitudes. Because the combination changes the apparent reflectance of each pixel, such a different distribution would lead to a systematic bias of if no corrections were made. The shape model of Ryugu (SHAPE_SFM_3M_v20180804) was used to calculate the incidence/emission angles.

Fig. 4 .

 4 Fig. 4. Schematic diagram of our calibration method. Based on low-resolution images (a, hyb2_onc_20180710_081932_tvf_l2b is shown as an example), a map of the standard deviation of pixel radiance (σ sp ) was produced (b). The σ sp map was then calibrated with the results of visual counting with high-resolution images (c, hyb2_onc_20180801_141202_tbf_l2b).

Fig.

  Fig.5. Concept of our method to convert σ sp of individual pixels to the coverage of boulders. We calculated σ sp from images for each region. We binarized the

Fig. 6 .

 6 Fig. 6. (a) Relationship of σ sp obtained in our method and coverage of boulders based on visual counting. The calibration line is used for converting σ sp to the coverage of boulders. (b) Relationship of the threshold of σ sp (%) and coefficient of determination. R 2 peaks at the threshold of 10%, indicating that the number of pixels showing σ sp > 10% is most correlated with the coverage of boulders larger than 1.5 pix.

Fig. 7 .

 7 Fig. 7. Comparison of coverage of 5-20 m boulders in different regions obtained with our method and visual counting (Michikami et al., 2019). (a) Comparison in different latitudinal bins and in (b) longitudinal bins. The two results are consistent within the uncertainties.

Fig. 8 .

 8 Fig. 8. (a) Global map of the standard deviation of radiance, which indicates the coverage of the boulders of 0.75-3 m. The map was constructed with the images taken at mid-altitude (August 1, 2018) and shape model of Ryugu. One pixel in the map is 0.5 m × 0.5 m. (b) Context image of Ryugu's surface. (c) Smoothest area found on the trough (200°E, 30°N) (hyb2_onc_20180801_180109_tvf_l2b). (d) Roughest area found at mid-latitude on the eastern hemisphere (110°, 50°N) (hyb2_onc_20180801_230638_tvf_l2b).

Fig. 9 .

 9 Fig. 9. (a) Longitudinal boulder distribution and profile of coverage of boulders (0.75-3 m). (b) Profile compared with Michikami et al. (2019) for boulder size in the range 5-20 m. Darker and lighter shades show 1σ and 2σ errors, respectively. The global map of σ sp shown in Fig. 8 was used to calculate the coverage of boulders (0.75-3 m and 5-20 m).

Fig. 10 .

 10 Fig. 10. (a) Latitudinal profile of boulders (0.75-3 m) obtained with our method. Lower latitudes exhibit lower boulder coverage. Darker and lighter shades show 1σ and 2σ errors, respectively. (b) Coverage of boulders in latitudinal bins. The results obtained using our method are compared with those obtained by Michikami et al. (2019) in the size range of 5-20 m. The global map of σ sp shown in Fig. 8 was used to calculate the coverage of boulders (0.75-3 m and 5-20 m).

Fig. 11 .

 11 Fig. 11. SFD of boulders in low-(20°S-20°N) and mid-latitudes (40°N-20°N, 40°S-20°S). Cumulative boulder number density data obtained in this study are plotted with those of previous studies. The data by Michikami et al. (2019) are plotted at the diameters 5 and 10 m. The power-law fitting lines (dashed lines) are constructed using the data points obtained in this study. Note the crossover at < 1 m. The global map of σ sp shown in Fig. 8 was used to calculate the coverage of boulders (0.75-3 m and 5-20 m).

Fig. 12 .

 12 Fig. 12. Distribution of power-law indices on Ryugu. One pixel in the map is 5 °× 5 °.

Fig. 13 .

 13 Fig. 13. (a) σ sp map for Itokawa. The regions with the lowest σ sp (~5%) coincide with smooth regions called Muses Sea (30-120°E, 0-60°S) and Sagamihara (160-240°E, 15°S-90°N). (b) Power-law index map for Itokawa. The smooth areas show high (steep) absolute value of power-law index (~-4.5). Power-law index of Itokawa is much more heterogeneous than that of Ryugu. One pixel in the maps is 3 °× 3 °. The images of Itokawa (ST_2420699983 -ST_2426434412) were used to produce these global maps.

Fig. 14 .

 14 Fig. 14. Typical boulder size (the size reaching the cumulative areal coverage of 50%) in different (a) latitudes and (b) longitudes.

Fig. 15 .

 15 Fig. 15. (a) Geological slope and (b) v-band albedo of Ryugu (Sugita et al., 2019). These parameters are compared with σ sp and α. The geological slope was calculated from the shape model (SHAPE_SFM_3M_v20180804). The v-band albedo map was produced by the global images of Ryugu obtained at HP (Sugita et al., 2019).

Fig. 16 .

 16 Fig. 16. Correlation between (a) power-law index and v-band albedo (45 °S to 45 °N), (b) σ sp and geological slopes in the two hemispheres. The global maps of power-law index α (Fig.8) and σ sp (Fig.12) were correlated with the geological slope derived from the shape model of Ryugu (SHAPE_SFM_3M_v20180804).

Ac knowledgements

This study wa s supported by KAKENHI from the Japanese Society for Promotion of Science (JSPS) (Grant Numbers 19K14778, 20H04607) and the JSPS Core-to-Core program "International Network of Planetary Sciences." P.M. acknowledges support from the French space agency, CNES and from the European Union's Horizon

Unci ted References

Michel and Richardson, 2014,

Decl aration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. 2020 research and innovation program under grant agreement No 870377 (project NEO-MAPP), and from Academies of Excellence: Complex systems and Space, environment, risk, and resilience, part of the IDEX JEDI of the Université Côte d'Azur.