Yutu-2 radar observation of the lunar regolith heterogeneity at the Chang’E-4 landing site

Context. The lunar penetrating radar (LPR) carried by the Yutu-2 rover performed the ﬁrst in situ measurement of the subsurface structure and physical properties of the subsurface materials on the far side of the Moon. It provides an unprecedented opportunity to study the formation and evolution of the lunar surface. Aims. This paper aims to quantitatively estimate the heterogeneity of the lunar regolith using the high-frequency Yutu-2 radar observation and constrain the modeling parameters (e.g., autocorrelation length) on a radar simulation. Methods. The heterogeneity of the lunar regolith was quantiﬁed by comparing the simulation and observation acquired by the high-frequency Yutu-2 radar within the ﬁrst 17 lunar days after its landing. The radar simulation was determined by the numerical calculation of the stochastic regolith model. The change in the autocorrelation length to the modeling was derived by calculating the coarseness of the model. Results. The disturbance range of the lunar regolith with a thickness of ∼ 12 m at the Chang’E-4 landing site is constrained to be ∼ 0.20 ± 0.06m, indicating a high self-similarity. The stochastic model’s spatial disturbance is controlled by the autocorrelation length and is also scaled by the model size, and the radar scattering echo strength decreases with the increase in autocorrelation length. Conclusions. We conclude that the heterogeneity of lunar regolith is positively related to the geological age. The application of the disturbance range at the decimeter scale might provide a valuable reference to assist in interpreting the radar observation data of the Moon (e.g., Arecibo radar, Min-SAR and Mini-RF, and in situ LPR). key to the model. This also discusses the autocorrelation length effect on radar simulations. Our results suggest that the disturbance is controlled by the autocorrelation length and is also scaled by the model size, and the radar scattering echo decreases with the increase in autocorrelation length. be


Introduction
The Moon is the closest satellite to the Earth and the prime target for human exploration of extraterrestrial celestial bodies. It is the end member of the evolution of terrestrial planets . Studying the geological process and evolution history of the Moon helps us to understand the past and future of the Earth's evolution better. In the exploration of the Moon, the lunar regolith is one of the most direct and primary scientific subjects (McKay et al. 1991). Generally, surface topography, chemical composition, and physical properties of the lunar regolith can be observed by the optical instruments (Lucey et al. 2000;Pieters et al. 2000;Blewett et al. 2005;Jin et al. 2015). However, they are unable to measure the subsurface materials directly, whereas the lunar ground penetrating radar can provide an unparalleled advantage Xiao et al. 2015;Li et al. 2020). (Huang et al. 2018;Xiao et al. 2021).
In 1972, the Apollo lunar sounder experiment (ALSE) onboard Apollo 17 was the first space-borne radar to detect the subsurface of the Moon (Porcello et al. 1974). Subsequently, the Japanese launched the SELenological and ENgineering Explorer (SELENE) spacecraft that carried the Lunar Radar Sounder (LRS) to the Moon in 2007, which is similar to the ALSE radar in terms of its operating frequency and detection modes (Ono et al. 2009). The difference is that the range resolution of the LRS radar (∼75 m at free space) is higher than that of the ALSE radar (Ono et al. 2009;Porcello et al. 1974). However, both the LRS and ALSE radars cannot detect the small details of the interior structure of the lunar regolith due to the limitation of range resolution.
In 2013, the Chinese Chang'E-3 mission carried the first in situ lunar penetrating radar (LPR) onboard the Yutu rover to probe the Moon's surface. The LPR provides an unprecedented opportunity to study the interior structure and physical properties of lunar subsurface materials at high range resolution Fa et al. 2015;Xiao et al. 2015;Feng et al. 2017;Ding et al. 2020a). Afterward, the Chinese Chang'E-4 mission successfully landed at the bottom of the Von Kármán crater on the far side of the Moon in 2019 ( Fig. 1A; Li et al. 2020;Ding et al. 2020b). Ray patterns of the ejecta from the Finsen crater extended to the Chang'E-4 landing site and modified the local surface materials Fig. 1B). The digital elevation model (DEM) data (Fig. 1C.) of the landing area illustrates that the elevation from the northwest to southeast is higher than that from the northeast to southwest. The identified linear geological features from the northeast to southwest suggest that they were formed and/or modified by the ejecta of the Finsen crater (Di et al. 2019;Xiao et al. 2021), shown in Fig. 1C.
In general, the kilometer-scale roughness of the Moon reflects the remnants of large craters or volcanic eruptions (Kreslavsky et al. 2013). The 100 m-scale roughness reflects the transformation of the Moon's surface by small impact craters (Kreslavsky et al. 2013). The meter-scale roughness is more directly related to the weathering degree of the interior lunar regolith and the degeneration of small craters (Helfenstein & Shepard 1999;Cai & Fa 2020). However, the decimeter-scale roughness of the interior lunar regolith has yet to be studied further (Fa et al. 2011). The roughness of the lunar surface can be derived using the optical images and the DEM data (Fa et al. 2011), but a quantitative study of the interior heterogeneity of the lunar regolith is beyond the capability of optical means because there is no penetration there. Here, the high-frequency channel of the Yutu-2 radar has a range resolution better than 0.3 m Su et al. 2014), which provides an opportunity to study it at the decimeter scale.
The lunar regolith is a globally distributed loose material layer over the bedrock on the Moon's surface (McKay et al. 1991). Its thickness distribution on the entire lunar surface is highly uneven; for example the average thickness of the regolith is ∼5 m in the mare region, but ∼12 m in the highland region (Shkuratov & Bondarenko 2001). A broad definition of the lunar regolith is that it encompasses all the weathered materials covering the bedrock, even including rocks several meters in diameter (McKay et al. 1991;Wilcox et al. 2005). A narrow definition of the lunar regolith would be that it includes the materials with the particle size < 1 cm in diameter, but > 1 cm are classified as rock fragments (McKay et al. 1991;Wilcox et al. 2005). The lunar regolith is a substance that develops on top of bedrock (e.g., basalt on the mare region) or ejecta materials from craters (McKay et al. 1991;Crawford et al. 2021). In general, the lunar regolith at the old geological unit is fine and homogeneous (McKay et al. 1991), and the lunar regolith at the young geological unit is a mixture of regolith and rock fragments (e.g., local impact debris and foreign ejecta materials) (McKay et al. 1991;Basilevsky et al. 2015). However, the recent studies on the lunar regolith indicate that its interior material is not homogeneous (Xiao et al. 2015;Feng et al. 2017;Ding et al. 2021). This is different from our past perception that it is highly homogeneous and even "transparent" to the electromagnetic sounding in the Apollo era (Carrier III et al. 1991).
The heterogeneity of the lunar regolith strongly affects the radar sounding over the lunar surface, especially for the radar operation on a high-frequency channel (e.g., ground-based Arecibo radar, space-borne Min-SAR and Mini-RF, and in situ LPR) (Fa & Wieczorek 2012;Shkuratov & Bondarenko 2001;Spudis et al. 2010;Stickle et al. 2016;Ding et al. 2020c;Zhang et al. 2021). The high heterogeneity of the lunar regolith leads to radar scattering being more complicated; for example, it elevates the radar circular polarization ratio (CPR) value, which might be cause to misjudge the existence of water ice in the anomalous craters on the Moon (Fa & Eke 2018;Liu et al. 2018). Therefore, it is an obstacle to interpreting the radar data more precisely, especially for a radar simulation since we do not know the background of the heterogeneity of the lunar regolith (e.g., the autocorrelation length) (Fa et al. 2011).
In this paper, the disturbance range is constrained to be ∼0.20 ± 0.06 m, and the empirical formula of the disturbance range and the geological age is determined. In addition, we also first report the influence of the autocorrelation length on modeling the stochastic regolith and the radar simulation. Our results show that the stochastic model's spatial disturbance is controlled by the autocorrelation length and it is also scaled by the model size; also, the radar scattering echo strength decreases with the increase in the autocorrelation length. The high-frequency LPR acquired the radar observation used in this study within the first 17 lunar days after the landing of Chang'E-4. The general structure of the paper is as follows: first, the high-frequency LPR data collection and radar description are introduced in Sect. 2.1. Second, an approach to construct the stochastic medium model is applied to characterize the interior heterogeneity of the lunar regolith (see Sect. 2.2). An algorithm for calculating the coarseness of the model is introduced in Sect. 2.3. The relationship between the coarseness and the autocorrelation length is determined in Sect. 3.1, and the effect of the autocorrelation length on the radar simulation is discussed in Sect. 4.1. Third, the disturbance range of the lunar regolith heterogeneity is obtained by comparing a radar simulation and observation seen in Sect. 3.3. Finally, the disturbance range of the lunar regolith is inferred to associate it with the geological age, and then its indications for lunar regolith heterogeneity are discussed in Sect. 4.2.

Data collection and Yutu-2 radar description
The Chang'E-4 spacecraft is a backup plan for the Chang'E-3, so that the scientific instruments of both missions are similar, including the LPR (Jia et al. 2018). The LPR onboard the Yutu-2 rover consists of two channels of different operating frequencies . The low-frequency channel operates at a central frequency of 60 MHz with a range resolution of ∼1 m (Xiao et al. 2015), and its penetration depth is > 100 m in the ground glacier experiment (Zhang et al. 2014). The scientific goal of the low-frequency LPR is to explore the subsurface structure of the lower lunar crust ). On the other hand, the highfrequency LPR operates at a central frequency of 500 MHz with a range resolution of <0.3 m in the regolith-like material (Zhang et al. 2014). It is mainly used to detect the interior structure and thickness of the lunar regolith ; for example the high-frequency LPR revealed three stratigraphic structures within the upper ∼10 m at the Chang'E-3 landing site (Fa et al. 2015;Xiao et al. 2015;Zhang et al. 2019;Xing et al. 2017), and the average permittivity within the upper ∼4 m was estimated to be ∼3.2 (Ding et al. 2020c), which is similar to that of typical lunar regolith . For the Chang'E-4 scenario, layering structures within a thickness of ∼45 m below the surface of the Moon have been observed by the high-frequency LPR system Xiao et al. 2021;Zhou et al. 2022).
The Yutu-2 radar had operated for 32 lunar days until July 2021 since its landing on the surface of the Moon. It acquired a large volume of scientific data archived in the Ground Application System, National Astronomical Observatories of the Chinese Academy of Science. Here, we only adopted the Level 2B (L2B) data of the high-frequency LPR to achieve our objective, which can be freely accessed online 1 . In addition, the Radargram obtained by the high-frequency LPR onboard the Yutu-2 rover within its first 17 lunar days. The x-axis is the radar survey distance (m). The right y-axis is the radar echoes' corresponding two-way traveling time (ns). The left y-axis is the derived depth (m), where the velocity of the radar pulse propagation into the lunar subsurface material was estimated to be 0.16 m/ns using the hyperbolic fitting method ). (A) The high-frequency LPR L2B radargram. (B) Processed radargram after background removal, band-pass filtering (250MHz-750MHz), and applied SEC gain for amplitude compensation. (C) The black dashed line plots the boundary between the upper fine-grained lunar regolith and the ejecta materials from the nearby craters. Our manual calculation of the dashed line is based on the amplitude of the interface reflection which has become sharply stronger, indicating the boundary of difference of materials. The average thickness of the upper fine regolith corresponds to ∼12 m. low-frequency LPR observation data are not considered in this study due to their radar wavelength being too long to detect the lunar regolith, especially for describing their interior materials in detail Su et al. 2014). During the LPR operation in the first 17 lunar days, the radar surveyed a total of 443.35 m along the Yutu-2 traveling route (e.g., Fig. 1D). We carried out 11,698 valid traces after repeating data removal, and they are arrayed, forming a radar image which is plotted in Fig. 2A. There are strong surface reflections and multireflections caused by multiple bounces of the radar pulse between the lunar surface and the Yutu-2 rover body (seen in Fig. 2A). These noise reflections submerged the subsurface reflections and resulted in the bad quality of the radargram. Therefore, we conducted data preprocessing to improve the data quality and enhance the radar echoes' signal-to-noise ratio (S/N). The preprocessing methods applied to the L2B radar data include background removal, bandpass filtering (with a band-pass of 250-750 MHz), and amplitude compensation with the spherical and exponential compensation (SEC) gain, where the background removal was obtained by subtracting the average value of each A-scan from the unprocessed data horizontally (Ding et al. 2020c). The preprocessing methods are similar to those we performed on the Chang'E-3 and the Chang'E-4 radar dataset in our previous studies (Ding et al. 2020c(Ding et al. ,b, 2021. Fig. 2B shows the processed radargram.

Stochastic medium model construction
The stochastic medium model describes the physical properties of natural lunar regolith better than the traditional homogeneous model ). Its numerical simulation has been successfully applied to the interpretation of seismic wave and ground penetrating radar data (Askar & Cakmak 1988;Irving & Holliger 2010;Jiang et al. 2013). Recently, the Chinese Chang'E-3 and Chang'E-4 in situ radars have been successfully deployed on the Moon's surface to probe the lunar regolith (Xiao et al. 2015;Li et al. 2020). One approach of the stochastic medium modeling for the lunar regolith has been introduced to assist in the interpretation of the LPR data Li et al. 2017;Wang et al. 2021). The stochastic medium model is capable of describing large-scale and small-scale heterogeneities (Ikelle et al. 1993;Ergintav & Canitez 1997;Ding et al. 2017). The large-scale heterogeneity indexes the average properties of the medium, and the small-scale one is determined by a random disturbance, indexing the model parameters of the medium, such as the autocorrelation length and the standard deviation disturbance Ergintav & Canitez 1997). Here, the stochastic regolith model is built to constrain the heterogeneity of the lunar regolith.
On the hypothesis of the stationary stochastic process, the stochastic permittivity model ε (r) of the lunar regolith can be expressed as wherer = xe x + ze z is the spatial position vector in the model, ε m is the average relative permittivity of the background material, δ (r) is the standard deviation disturbance, ε f (r) is the random autocorrelation function.
Steps to construct the stochastic regolith model ε (r) are described as follows, and the corresponding block scheme is presented in Fig  (1) First of all, we selected an Gaussian autocorrelation function r (r) to obtain the random autocorrelation function ε f (r) of the model, which is exhibited as where a, b are horizontal and vertical autocorrelation lengths, respectively. The Gaussian autocorrelation r (r) is plotted in Fig. 4A.
(2) Secondly, the power spectral density R k can be determined by the Fourier transformation of the r (r): wherek = k x e x + k z e z represents the vector of the frequency component, and M and N are the size of the model after discretization in a row and column, respectively.
(3) Thirdly, to obtain the random power spectral density W k , a random phase ϕ k is added to the power spectral density (Eq. (3)). The random phase follows a Gaussian distribution ( Fig. 4B) with the range of [0, 2π): (4) According to the Wiener-Khinchin theorem, the power spectral density is the Fourier transform of its autocorrelation function (Cohen 1998;Shinozuka & Deodatis 1991). Therefore, the random autocorrelation ε f (r) can be obtained by the inverse Fourier transform of the random power spectral density W k , expressed as (5) Finally, the random autocorrelation ε f (r) should be normalized (Fig. 4C), and then substituted into the Eq. (1). The stochastic medium model is finally generated and shown in Fig. 4D.

Calculation of the coarseness of the stochastic medium model
In constructing the stochastic medium model, the autocorrelation length is an important parameter and it affects the range of horizontal or vertical disturbance ). However, the influence of the autocorrelation length on modeling has rarely been studied (Li & Liu 2011). Here, we introduce texture coarseness to quantify the influence on modeling. For simplicity, the vertical and horizontal autocorrelation lengths are assumed to be the same in this paper. The relationship between the coarseness and the autocorrelation length can be obtained by running the models with different parameters.
In the stochastic medium model, the distribution of relative permittivity fluctuates with the spatial disturbance (e.g., Fig. 4D). Various forms of the stochastic medium model can be generated by setting up different autocorrelation lengths (Ergintav & Canitez 1997;Ding et al. 2017). The coarseness is proposed by Tamura et al. (1978), which is considered a basic and important texture feature according to the study of visual perception (Tamura et al. 1978). The texture feature refers to the variation in the image grayscale level (e.g., Collewet et al. 2004). The stochastic regolith model is composed of a relative permittivity variation. Similarly, the relative permittivity variation can be assimilated with the image grayscale level. Therefore, we consider the relative permittivity variation in the model to be a texture feature. The following steps are introduced to calculate its coarseness (Tamura et al. 1978).
(1) The average permittivity of the model ε k (r) is calculated with a sliding window (2 k × 2 k ), and its equation is expressed as wherer = xe x + ze z ,l = ie x + je z , k = 0, 1, 2, ..., k max , k max is the maximum scale of the sliding window, and g l is the value of relative permittivity at each point (x, z) in the model.
(2) For each point of relative permittivity in the model, we calculated the average value with content of 2 k × 2 k between the nonoverlapping range in horizontal and vertical directions. Eq. (7) ε kh represents the nonoverlapping range in the horizontal direction, and the equation (8) ε kv is the nonoverlapping range in the vertical direction, which is defined as (3) Eq. (9) is used to calculate the maximum average permittivity in step 2, which can be expressed as After we obtained the maximum value of relative permittivity using Eq. (9), its corresponding parameter k was substituted into the formula S best (r) = 2 k to derive the optimum value.
(4) Finally, the model coarseness F coar can be calculated by averaging all the optimum values, which can be expressed as where M and N are the size of the model in a row and column, respectively.

Relationship between autocorrelation length and coarseness
According to the modeling method described in Sect. 2.2, we established four model scenarios with different sizes. The sizes of the models are 2 m × 2 m, 4 m × 4 m, 6 m × 6 m, and 8 m × 8 m. The interval is set to 0.01 m. The value of autocorrelation length starts from 0.01 m to the length of the model size. Figure 5 shows the relationship between the autocorrelation length and the coarseness with the different sizes of the models since we can see that the variation of coarseness shows a similar trend in the four scenarios of the models (Fig. 5). Here, Fig. 5A is taken as an example to describe the coarseness variation. As the autocorrelation length increases, the value of coarseness increases rapidly until the peak position (e.g., ∼0.34 m in Fig. 5A). Then, the coarseness beyond the peak position exponentially decreases until the autocorrelation length at the position of ∼0.69 m. Finally, the coarseness reaches a steady trend with a slight decline.
The coarseness variation of the four scenarios of the model suggests that the coarseness is not only controlled by the autocorrelation length, but it is also scaled by the model size.

Multiple reflections between the lunar surface and the bottom of the Yutu-2 rover
The high-frequency LPR is mounted at the bottom of the Yutu-2 rover and is 0.3 m above the lunar surface Li et al. 2020). The transmitter antenna, which is not encased in a shell, sends a Gauss-like pulse down to the lunar surface . When the radar pulse touches the Moon's surface, one part of the energy is reflected by the surface, and the other part is transmitted. When the reflected energy bounces to the Yutu-2 rover, it is reflected again, caused by the bottom of the rover, which is covered by a metal-like material. The radar antennas can easily acquire these multiple reflections. Surface roughness at a centimeter scale has a slight effect on the radar echoes due to the short wavelength of the high-frequency LPR . The geometry between the lunar surface and the bottom of the Yutu-2 rover is relatively fixed. As a result, the jamming echoes show, as horizontal features, spacing at a fixed interval in the radargram (Fig. 6C). Here, we use the radar data without background removal to compare the simulation to observation because the amplitude contrast of the radargram is degraded after background removal. Therefore, to avoid interference due to multiple reflections, we chose to use only the radar observation with the two-way traveling time >50 ns for the comparative analysis (e.g., the green arrows point at a time delay of ∼50 ns; Fig. 6B). The red dashed line in Fig. 6B identifies the interface between the lunar regolith and its underlying coarse-grained ejecta materials at a depth of ∼12 m Xiao et al. 2021). Here, we extracted the envelope of the radar echoes using the Hilbert transform, amplified the weak signal, and binarized the radargram to make the interface reflections more distinctive. The processed radargram is shown in Fig. 6A.

Constraining the heterogeneity of the interior lunar regolith
Through the construction and numerical calculation of many stochastic regolith models, the simulations were carried out and then compared to the radar observations. The optimum autocorrelation length was estimated by constraining the best fit model, which can describe the heterogeneity of the lunar regolith. It is important to note that we only compared the radar observation within the fine-grained lunar regolith to the simulations since the autocorrelation length is the semi-axis of the autocorrelation function r (r). Here we specify that the disturbance range is double the autocorrelation length.
In the modeling, the parameters of the regolith model were set as follows. The lunar regolith model size was set to 15 m in vertical depth and 7 m in horizontal distance since the thickness of the fine-grained lunar regolith at the Chang'E-4 landing site was measured to be ∼12 m (e.g., Fig. 2C; Li et al. 2020). The relative permittivity of the background medium was assumed to be ∼3.5 with a standard deviation of 10%, and the loss tangent of the model was set to 0.004, which are consistent with the results estimated at the Chang'E-4 landing site ). The autocorrelation length was set from 0.01 m to 1.5 m with an interval of 0.01 m so that a total of 150 sets of models were constructed. The autocorrelation length of 1.5 m corresponds to 21.4% of the model size in horizontal direction, which is beyond the peak value of coarseness (normally ∼15.3% of the model size). It builds all possibilities of the regolith model in this case. Figure 7 shows three examples of the lunar regolith model.
In the numerical calculation, we adopted a frequency domain method to calculate the regolith models one by one (Bitri & Grandjean 1998). The working frequency of the radar pulse was assumed to be 500 MHz, which is consistent with that of the high-frequency channel of the Yutu-2 radar. The radar transmitting wave was set to a Gaussian waveform . The time window for the simulation was set to 150 ns. The spacing interval of two adjacent traces was set to 0.035 m to get a simulation with 200 A-scan traces.
After calculating the entire regolith models (150 sets of models), we obtained the corresponding 30 000 A-scan traces. We assumed that the simulation data cover all possible responses from an observation on the lunar regolith at the Chang'E-4 landing site because all possibilities of stochastic regolith models were considered.
The radar observations used in this paper have 11 698 valid A-scan traces. The steps for calculating the optimal autocorrelation length for the best fit model are described as follows.
Step 1. Both the radar observation and the simulation are normalized with their peak value of the surface reflection echo.
Step 2. An A-scan of a radar observation must be selected, for example beginning from the first trace of the observation (11 698 traces in this case). The root mean square error (RMSE) is a loop calculating between the A-scan of the radar observation and simulation (e.g., starting from the simulation with the minimum autocorrelation length). Then the minimum RMSE is recorded as the best result. The smallest RMSE indicates the simulation's best fit for the observation.
Step 3. One should loop to the next simulation and repeat step 2. The RMSE values can be obtained after completing all 150 simulations.
Step 4. One should loop to the next A-scan of a radar observation, and redo steps 2 and 3. One should end the algorithm until all traces of the radar observation have been processed. Fig. 8A-D shows a case of an A-scan radar observation labeled 8910 to constrain the optimal disturbance range. Figure 8B is plotted concerning steps 2 and steps 3. It shows that the RMSE is a function of the disturbance range (Fig. 8A). It can be seen that the lowest RMSE occurs when the disturbance range is at ∼0.21 m in Fig. 8A. This means that the heterogeneity of the lunar regolith indicates a random disturbance range of ∼0.21 m. On the other hand, the mismatch cases between the radar simulation and observation are shown in Fig. 8C with a disturbance range of 0.1 m, and Fig. 8D with a disturbance range of 0.5 m.
After looping the calculation of all radar observations, we carried out 11 698 optimal disturbance ranges, and the corresponding histogram is plotted in Fig. 8E. The statistical result shows that the disturbance range is ∼0.20 m on median and with a 1δ standard deviation of 0.06 m, and the 14-84% percentile range is from ∼0.05 m to ∼0.30 m (Fig. 8E). In other words, the disturbance range of the lunar regolith is constrained to be ∼0.20 ± 0.06 m at the Chang'E-4 landing site. In addition, we also adopted the other metric, the relative error (RE; where s (i) is the simulated data, and o (i) is the observed data.

Implications of autocorrelation length on the radar simulation
It is generally believed that the spatial disturbance of the stochastic model is intensified by the increasing autocorrelation length (Xi & Yao 2004;Li & Liu 2011;Zhang et al. 2018). However, our results in Sect. 3.1 suggest that the spatial disturbance is not only controlled by the autocorrelation length, but it is also scaled by the model size.
The autocorrelation length is the parameter that directly relates to the description of the inhomogeneity of the stochastic medium model. To better understand the effect of the autocorrelation length on radar simulations, we proposed three problems listed as follows and discuss them each individually.
1. The effect of autocorrelation length on the stochastic model. For this purpose, we established three sets of stochastic medium models with the autocorrelation lengths of 0.31 m, 0.61 m, and 0.92 m, plotted in Figs. 9A-C. Relative permittivity along the profiles in the corresponding models are plotted in Figs. 9D-F.
The permittivity variation of the one-dimensional regolith models shows that the number of peaks in Fig. 9D, Fig. 9E, and Fig. 9F are 4, 2, and 1, respectively. The spatial disturbance region of permittivity in Fig. 9F is more significant than that in Fig. 9D and Fig. 9E. Therefore, we can infer that the permittivity variation of Fig. 9F is more aggregated than those of the others, indicating a higher self-similarity. The peak width of the permittivity variation in Fig. 9 is limited by the length of the autocorrelation, for example there are two peaks of permittivity variation in Fig. 9E, and the width of which is about 1 m. The peak width is used to describe the degree of spatial disturbance of the regolith model, which is positively correlated with A43, page 9 of 17 The third row is the exponential model whose autocorrelation length is the same as in the first row. The fourth row is for one-dimensional models extracted from the panel in the third row.
double the autocorrelation length (e.g., 0.61 m in Fig. 9E). The peak width can be described by the disturbance range as well. However, the peak width, otherwise referred to as the disturbance range, seems twice as long as the autocorrelation length in Fig. 9F, and the one-dimensional permittivity variation appears to be unsmooth compared to Figs. 9D and E. The reason for this unsmoothness at the micro-scale is possibly associated with the value of the autocorrelation length beyond the maximum coarseness of its model. The model size scales its autocorrelation length at 23%, but the maximum coarseness at ∼15.3% is obtained in Sect. 3.1. To avoid this situation, we suggest that the autocorrelation length is preferably within its model size of first ∼15.3% in the stochastic medium model construction.
To answer the question, with the increase in autocorrelation length, the self-similarity of materials in the model increases, which shows that "agglomeration" is strengthened. The microscale disturbance appears in the model, which may be related to the fact that the autocorrelation length exceeds the first ∼15.3% model size.
2. The generality of the function of coarseness and the autocorrelation length of the Gaussian stochastic model. Generally, the modeling methods of stochastic media include not only the Gaussian type, but also the exponential type (Ikelle et al. 1993;Ergintav & Canitez 1997;Ding et al. 2017). We only consider the Gaussian stochastic model in Sect. 2.2, but not the exponential method. The exponential method is suitable for the situation where the dielectric properties are discontinuous at the microscale, but there is self-similarity at the macro-scale (Ikelle et al. 1993).
Here, we use the same algorithm described in Sect. 2.3 to establish the exponential stochastic models. The coarseness concerning the autocorrelation length is shown in Fig. A.1 (see  supplementary materials). Comparing Fig. A.1 and Fig. 5, we found that the rising areas of the two are almost consistent. However, the coarseness did not decrease after reaching the maximum value, it but remained stable (Fig. A.1). To explain the possible reason for this phenomenon, we constructed three exponential stochastic models with the same autocorrelation length as in Figs and 9I. Therefore, we reasonably infer that with the increase in the autocorrelation length, the internal heterogeneity changes only within the first ∼15% of the model size of the autocorrelation length in the exponential model case. If the autocorrelation length exceeds this ratio, the heterogeneity of the model presents a steady state. Therefore, the generality of the Gaussian stochastic model is only adequate when the autocorrelation length is within the first 15% of the model size.
3. The effect of autocorrelation length on the radar scattering. The above two problems discussed the modeling part of radar simulations. Here we propose to discuss the influence of the autocorrelation length on radar scattering. To this end, we established three lunar regolith models with a size of 6 m × 6 m. The model consists of three components, vacuum, lunar regolith, and bedrock, and the autocorrelation length is assumed to be 0.03 m, 0.09 m, and 0.15 m, respectively (for the rest of the model parameters, see the caption of Fig. 10). Figures 10D to 10F show the A-scan radar simulations of models in Figs. 10A-C. We found that the overall amplitude of radar scattering echo in the regolith part (black rectangles in Fig. 10) decreases from Figs. 10D to F. The reason caused the amplitude enhancement of radar scattering is more peaks of permittivity protrusion in the stochastic regolith model, indicating the dielectric discontinuity intensified. According to the principle of Yutu-2 radar detection, once the radar pulse meets the interface of dielectric discontinuity, it produces reflection, transmission, and scattering ). If the radar pulse propagates into a regolith with a short autocorrelation length, it may encounter a high probability of dielectric discontinuity, resulting in strengthened radar scattering echoes. This also explains why the amplitude of the radar scattering within the regolith in Fig. 10D is stronger than those of others. This allows us to infer that the reflection intensity of radar scattering decreases with the increase in autocorrelation length.
To sum up, we can draw the following conclusions: a. With the increase in the autocorrelation length, the selfsimilarity of the materials in the stochastic model became stronger, showing that the agglomeration was strengthened.
b. In constructing the stochastic model, the value of autocorrelation length should not exceed ∼15% of its model scale.
c. The radar scattering intensity becomes weaker with the increase in autocorrelation length.
In a one-sentence summary, the self-similarity of materials in the model is strengthened and the radar scattering echoes decrease with the increase in autocorrelation length. It is important to note that we both use the autocorrelation length and disturbance range to describe the unhomogeneous model in this paper, but the physical meaning is the same. We chose to use the disturbance range instead of the autocorrelation length to describe the heterogeneity of the lunar regolith because it intuitively allows the potential readers to understand the spatial disturbance of the lunar regolith rather than the autocorrelation length.

Indications of the disturbance range to regolith heterogeneity
The disturbance range in the lunar regolith at the Chang'E-4 landing site is constrained to be ∼0.20 m in Sect. 3.3. Some geological activities are related to this disturbance range to heterogeneity in the lunar regolith; for example, the geological activities that modified the surface of the Moon include weathering, volcanic activity, and tectonic movement (Head III 1976;Jaumann et al. 2012). The lava flow produced by volcanic activity mainly covers the lunar regolith (Head III 1976;Heiken et al. 1991), and it slightly affects its interior, so the situation was excluded. The tectonic movement of the Moon mainly changes the activity of the lunar crust at a macro-scale (Byrne 2020), so it was not considered either. Weathering has been a long-term and continuously geological surface process on the Moon, which plays a decisive role in the modification of the lunar regolith (McKay et al. 1991;Xiao et al. 2013). The meter-scale roughness of the lunar surface is directly related to the weathering degree of the interior lunar regolith and the degradation of small craters by the optical observation (Helfenstein & Shepard 1999;Cai & Fa 2020). This implies that the disturbance range of the lunar regolith may be related to its weathering degree. The lunar regolith is the product of long-term weathering of surface materials of the Moon, and its degree of weathering is related to the exposure time (McKay et al. 1991). Therefore, a function of the exposure time and the disturbance range of the regolith may be obtained by a radar measurement on the different geological units.
The optical techniques cannot directly quantify the disturbance range of the interior lunar regolith except for the Chang'E-3 radar. In 2013, the Chang'E-3 mission landed northeast of Mare Imbrium, and the geological age of the landing area was constrained to ∼2.5 Ga (Xiao et al. 2015;Ding et al. 2020c). This provides a possibility to discuss the relationship between the autocorrelation length of the lunar regolith and the geological age. Figure 11A shows the radar observation of the Chang'E-3 high-frequency channel within the first two lunar days. Using the same algorithm, we carried out the distribution of the disturbance range of the lunar regolith at the Chang'E-3 landing site (see in Fig. 11D). The statistical results show that the disturbance range derived by the Chang'E-3 radar observation is ∼0.13 m on the median, and 1σ standard deviation is 0.06 m. In addition, we also adopted the relative error instead of RMSE to calculate the deviation of the radar observation and simulation (see in Fig. A.3). The estimated result shows almost no difference from those of Fig. 11D.
The Yutu radar surveyed on the continuous ejecta blanket of the ZiWei crater (∼450 m in diameter; Xiao et al. 2015). Although the Ziwei crater was formed ∼26.7 Ma ago, its ejecta is dominated by the lunar regolith, which develops on top of the Eratosonian-aged basalt (Sharpton 2014;Ding et al. 2020c). The geological age of Chang'E-4 is dated to ∼3.7 Ga (Huang et al. 2018;Chang et al. 2021), which is older than that of Chang'E-3 A43, page 12 of 17 C.    (∼2.5 Ga). Therefore, we infer that the disturbance range of the lunar regolith might be positively correlated with age. It is worth noting that neither the Chang'E-4 radar observation (e.g., Fig. 8B) nor the Chang'E-3 radar observations (e.g., Fig. 11C) are completely matched in comparison to the simulations. It is consistent with the internal scenario of the lunar regolith. Since the lunar regolith is a product of impacts (McKay et al. 1991), it is formed and developed on top of the basalt and/or ejecta materials. Therefore, there is the inevitable existence of rock fragments in it (McKay et al. 1991). However, the stochastic model hypothesizes that the lunar regolith does not contain rock fragments. Let us consider adding the rock fragment factor into the mathematical modeling algorithm. This makes it impossible to predict the model because it is scarcely possible to determine the buried rock fragment's spatial position, geometrical shape, size, dielectric properties (e.g., relative permittivity and loss tangent), etc.
The lunar crater chronology suggests that the lunar impact flux was nearly constant after ∼3.7 Ga (Robbins 2014;Xie et al. 2021). Therefore, we use a linear fitting for the relationship between the geological age and the disturbance range. The linear fitting result is y = (0.0536 ± 0.0139) x, R 2 = 0.99 where x represents the geological age and y is the disturbance range of the lunar regolith. This empirical formula can predict the disturbance range of the lunar regolith associated with the geological units of different ages. However, the applicability and generality of the empirical formula still need more observation data to improve. Figure 12 illustrates a possible scenario of the empirical formula application at the Chang'E-4 landing site. The top of basalt (∼3.7 Ga) on the bottom of the Von Kármán crater is covered by the ejecta materials from the Orientale basin and/or Finsen during the period of ∼3.68 Ga to ∼3.0 Ga (Figs. 12A and B). During this period, the lunar regolith had been weathered on top of the ejecta for ∼0.68 billion years. The disturbance range of the heterogeneity of the lunar regolith can be estimated to be ∼0.04 m using the empirical formula.
After that, no large amount of foreign material was transported to the landing area by the large impacts within ∼2 Ga until the early Copernican (∼1.0 Ga). However, most of the small impacts only excavated local materials and did not even penetrate the bedrock (McKay et al. 1991). The materials of the local lunar regolith are self-weathered, and its disturbance range is speculated to be ∼0.14 m (Fig. 12C), similar to that of the Chang'E-3 scenario.
Since entering the Copernican period, the shallow regolith on the Moon's surface has been continuously grinded and modified by impacts (mostly small impacts), and the weathering inside the lunar regolith has continued. Under the crushing and grinding by impacts, the local regolith material gradually becomes more fine-grained. It leads to an increase in material internal similarity in space. Finally, the disturbance range is speculated to be ∼0.20 m (Fig. 12D), suggesting a high self-similarity.

Conclusions
The high-frequency lunar penetrating radar onboard the Chang'E-4 mission provides an opportunity to study the heterogeneity of the lunar regolith at the decimeter scale. The radar A43, page 13 of 17 A&A 664, A43 (2022) simulation is obtained by numerically calculating the stochastic medium model, which is an effective method and suitable for describing complicated scenarios such as the lunar regolith compared to the traditional homogeneous model. The dielectric variations in the stochastic medium model are controlled by the modeling parameters, such as the background relative permittivity, autocorrelation lengths, and standard deviation. Choosing an appropriate value for the autocorrelation length is key to generating the final model. This paper also discusses the autocorrelation length effect on radar simulations. Our results suggest that the spatial disturbance is controlled by the autocorrelation length and is also scaled by the model size, and the radar scattering echo strength decreases with the increase in autocorrelation length.
The heterogeneity of fine-grained lunar regolith within the upper ∼12 m was estimated by comparing the radar observation and simulation. The disturbance range of the Chang'E-4 regolith was constrained to be ∼0.20 ± 0.06 m. The disturbance range is positively associated with the geological age. Its relationship with the geological age is discussed after the autocorrelation length of Chang'E-3 was constrained. Meanwhile, an empirical formula was derived to predict the disturbance range associated with the geological units of different ages, but its generality needs to be improved using more observations. The application of the disturbance range within the lunar regolith at the decimeter scale might provide a valuable reference to assist in interpreting the ground-based and space-borne radar data of the Moon (e.g., Mini-SAR and Mini-RF). Our results can also be broadly interesting to other airless bodies.