Self-Organization Characteristics of Lunar Regolith Inferred by Yutu-2 Lunar Penetrating Radar

Most previous studies tend to simplify the lunar regolith as a homogeneous medium. However, the lunar regolith is not completely homogeneous, because there are weak reflections from the lunar regolith layer. In this study, we examined the weak heterogeneity of the lunar regolith layer using a self-organization model by matching the reflection pattern of both the lunar regolith layer and the top of the ejecta layer. After a series of numerical experiments, synthetic results show great consistency with the observed Chang’E-4 lunar penetrating radar data and provide some constraints on the range of controlling parameters of the exponential self-organization model. The root mean square permittivity perturbation is estimated to be about 3% and the correlation distance is about 5–10 cm. Additionally, the upper layer of ejecta has about 1–2 rocks per square meter, and the rock diameter is about 20–30 cm. These parameters are helpful for further study of structural characteristics and the evolution process of the lunar regolith. The relatively small correlation distance and root mean square perturbation in the regolith indicate that the regolith is mature. The weak reflections within the regolith are more likely to be due to structural changes rather than material composition changes.


Introduction
The structure of the lunar shallow surface is complex due to long-term impacting, sputtering and stacking processes; in addition, the Moon has no strong cementation due to the absence of interstitial water. Thus, the heterogeneity on the Moon can be well retained and is much stronger than that on the Earth. The remote-sensing spectrum and radar detection results show that there is a thick regolith layer on the top of the Von Kármán crater [1][2][3]. This provides an excellent opportunity for detection of the structural characteristics and the stacking mode of the lunar regolith, which are critical for understanding the formation process of the lunar regolith. However, in the past, the lunar regolith has usually been assumed to be a uniform layered medium [4][5][6][7][8]. A series of recent studies have shown that the lunar regolith is essentially inhomogeneous [9][10][11][12][13][14][15][16] but the specific structural characteristics of the lunar regolith are still unclear.
Radar is widely used in the exploration of the Moon and Mars, and has become a necessary tool to detect the underground structure of extraterrestrial bodies [17][18][19][20][21][22][23]. The exploration of the Moon, Mars, and other planets with penetrating radars may benefit from the long experience in using the same approaches (the same type of radars and models to interpret the measurements) in studies carried out on Earth. The Yutu-2 rover traveled a long distance (>500 m) with a fine spatial sampling interval (~3-5 cm) on the far side of the Moon. Such an observation provides a significant opportunity to explore the local fine structure under the rover's path. The lunar regolith at the landing site has a thickness of up to~12 m [24][25][26], which is much thicker than that detected on the near side of the Moon by the Chang'E-3 lunar penetrating radar (LPR) [27,28].
Surprisingly, even though the thickness of the lunar regolith is high (>10 m), the Chang'E-4 LPR profile still shows some weak reflections in the lunar regolith, and their amplitudes are far smaller than those of the reflections from the underlying ejecta. This indicates that the lunar regolith is not homogeneous but weakly heterogeneous. The weak heterogeneity of the lunar regolith does not seriously influence the imaging and inversion [24,29,30] because it is too weak compared with the reflections from the ejecta. However, further study of these weak heterogeneous characteristics would be helpful in revealing the process of lunar regolith formation, particularly to the most recent stage of local small impacts [31]. Of course, long-distance heavy bombardments may also convey an amount of small-particle ejecta and contribute to such an ultra-thick layer of fine regolith on the far side of the Moon [32].
Most traditional numerical simulations of LPR data assume layered models [4][5][6][7][8] and cannot well represent the heterogeneity of the lunar regolith. Recently, several works have considered the heterogeneity [9][10][11][12][13][14][15] but do not take typical models into account (e.g., cracks and ejecta). Lv et al. [16] constructed four typical models and performed numerical simulations using the time-domain finite-difference method, which can better present the wave propagation in the lunar regolith with strong heterogeneity. Some works noted the existence of weak heterogeneity of the lunar regolith [11,[33][34][35]. However, to date, no numerical simulation has been conducted specifically on the weak heterogeneity of the lunar regolith, probably due to three aspects: 1. lack of methods describing heterogeneous media; 2. technical difficulty in constructing random models; 3. excessive computational cost of numerical simulations [16] compared with traditional methods. In this study, we assumed that the weak heterogeneity of the lunar regolith obeys the self-organization model [36], which is basically a stochastic approach for dealing with scattering waves through randomly inhomogeneous structures and is well known in the seismology field; then, we used the finite-difference time-domain method [16] to carry out a series of numerical simulations using various models; next, we analyzed the influence of different parameters on the weak reflection characteristics; finally, we constrained the possible range of each controlling parameter through grid searching techniques by matching the numerical simulation results with the Yutu-2 LPR observations, in terms of waveforms, envelopes, and reflection patterns.

Modeling and Simulation
The weak heterogeneity of the lunar regolith can be described by a stochastic model [36]. We keep a small part of the underlying ejecta by adding some sputtering blocks to the stochastic model [16]. This can provide a reliable reference for checking the amplitude variations from weak reflections of the lunar regolith to strong reflections of ejecta.
First, we generate a 2D dielectric permittivity background according to the typical 1D lunar regolith model [37][38][39] (Figure 1) ε r = 1.919 ρ and ρ = 1.92(z + 12.2)/(z + 18), where ε r is the relative dielectric permittivity, ρ (g/cm 3 ) is the density, and z (cm) is the depth. The velocity of radar wave propagation within the lunar regolith can be converted from the dielectric permittivity as v = c/ √ ε r , where v is the velocity of the radar wave and c is the speed of light in vacuum [40]. Then, we add some random perturbations and some ejecta to the background model to generate the model.
The 2D model has a width of 15 m and a depth of 16 m ( Figure 2). From top to bottom, it includes a 0.3 m thick vacuum medium (with a relative dielectric permittivity of 1) between the rover bottom and the lunar surface, a 12 m thick random medium representing the lunar regolith, and a~4 m thick random medium with some rocky blocks (i.e., ejecta layer). The numerical simulation method and parameter settings are consistent with those of Lv et al. [16]. The finite-difference method [41,42] is used to discretization the spatial partial derivatives of Maxwell's curl equation [43,44], and only the 2D transverse electric pattern is considered. The radar frequency is 400 MHz to avoid visible numerical  with those of Lv et al. [16]. The finite-difference method [41,42] is used to discretization the spatial partial derivatives of Maxwell's curl equation [43,44], and only the 2D transverse electric pattern is considered. The radar frequency is 400 MHz to avoid visible numerical dispersion. We can analyze the reflection characteristics of the lunar regolith by either a snapshot of the wavefield at a given time or a waveform at a given position.

Self-Organization Model
The self-organization model is an important approach for describing stochastic media. The controlling parameters generally include the mean value, variance (or standard deviation), and autocorrelation function (ACF). The ACF and its corresponding power spectral density function (PSDF) are the two most important parameters to reflect the distribution characteristics of random media. The most widely used autocorrelation functions include Gaussian, exponential, and Von Kármán types [45][46][47][48]. In the stochastic process theory, the PSDF is the Fourier transformation of ACF [49]. To characterize the random media, we define the ACF of the fractional permittivity fluctuation as: where non-dimensional quantity ξ(x) ≡ δε r (x)/ε 0 r is the fractional fluctuation of the permittivity, ε 0 r is the mean permittivity, and δε r is the perturbation of permittivity. This provides a statistical measure of the spatial scale and the magnitude of medium inhomogeneity.
The most tractable ACF is Gaussian ACF (Figure 2a,d), as given by: where a is the correlation distance. For the 3D case, the PSDF is: where m is the wavenumber vector and m = |m|. The Gaussian ACF is used to describe the medium with poor short-wavelength components because the PSDF goes rapidly to zero for a large m. Many models are based on the Gaussian ACF due to its solid mathematical background [50,51]. Another commonly used model is the exponential ACF (Figure 2b,e): For the 3D case, the PSDF is There is an extension of the exponential ACF, called Von Kármán ACF (Figure 2c,f): where Γ(κ) is the gamma function and K κ is the κ-th order modified Bessel function of the second kind. For the 3D case, the PSDF is: Compared with the Gaussian ACF, the Von Kármán ACF (including exponential ACF) is more suitable for describing the random permittivity inhomogeneity of the lunar regolith due to the power law property of the Von Kármán ACF [36].
We assume that the randomness is spatially stationary, and we individually add Gaussian, exponential, and Von Kármán (zero-order) random fractional fluctuations to the background model. The ejecta are then added at the bottom of the lunar regolith ( Figure 2). Figure 3 shows the numerical simulation results of these three models. Because of the poor short-wavelength component of Gaussian ACF random media, their reflection characteristics show more long-wavelength structures. Compared with the Gaussian media, the Von Kármán media (including the exponential type) show more short-wavelength structures. Due to the intense sputtering and stacking process, the heterogeneity of the lunar regolith should be more consistent with the Von Kármán (including the exponential type) media [35,36,52]. However, we cannot well distinguish its difference from the exponential ACF due to the limited radar wave resolution, even though the zero-order Von Kármán model has greater heterogeneity at small scales. Therefore, we choose the exponential ACF to construct random media in the following experiments. This can simplify the model and facilitate the discussion of the results.
We assume that the randomness is spatially stationary, and we individually add Gaussian, exponential, and Von Kármán (zero-order) random fractional fluctuations to the background model. The ejecta are then added at the bottom of the lunar regolith ( Figure 2). Figure 3 shows the numerical simulation results of these three models. Because of the poor short-wavelength component of Gaussian ACF random media, their reflection characteristics show more long-wavelength structures. Compared with the Gaussian media, the Von Kármán media (including the exponential type) show more short-wavelength structures. Due to the intense sputtering and stacking process, the heterogeneity of the lunar regolith should be more consistent with the Von Kármán (including the exponential type) media [35,36,52]. However, we cannot well distinguish its difference from the exponential ACF due to the limited radar wave resolution, even though the zero-order Von Kármán model has greater heterogeneity at small scales. Therefore, we choose the exponential ACF to construct random media in the following experiments. This can simplify the model and facilitate the discussion of the results.

The Model of Lunar Regolith
The correlation distance and root mean square (RMS) perturbation are two key parameters to describe the self-organization random media; thus, we concentrate on these two parameters in the numerical simulations. These two parameters have a relatively clear physical meaning: the correlation distance can well constrain the variation range of regolith permittivity and the RMS perturbation can well present the extent of regolith permittivity variation. Different groups of correlation distance and RMS perturbation limit the physical properties and structural characteristics of the lunar regolith.
We established different models of 5-20 cm correlation distance and 1-7% RMS perturbation. Figure 4 shows three exponential ACF random media with different correlation distances, namely, 5, 10, and 20 cm, respectively. Figure 5 shows three exponential ACF random media with different RMS perturbations of 3, 5, and 7%, respectively. Sputtering stones were also added below the 12 m depth of the lunar regolith.
For the sake of simplicity, we only select the critical groups of controlling parameters by plotting them together, as shown in Figure 6, although we actually performed a

The Model of Lunar Regolith
The correlation distance and root mean square (RMS) perturbation are two key parameters to describe the self-organization random media; thus, we concentrate on these two parameters in the numerical simulations. These two parameters have a relatively clear physical meaning: the correlation distance can well constrain the variation range of regolith permittivity and the RMS perturbation can well present the extent of regolith permittivity variation. Different groups of correlation distance and RMS perturbation limit the physical properties and structural characteristics of the lunar regolith.
We established different models of 5-20 cm correlation distance and 1-7% RMS perturbation. Figure 4 shows three exponential ACF random media with different correlation distances, namely, 5, 10, and 20 cm, respectively. Figure 5 shows three exponential ACF random media with different RMS perturbations of 3, 5, and 7%, respectively. Sputtering stones were also added below the 12 m depth of the lunar regolith.
For the sake of simplicity, we only select the critical groups of controlling parameters by plotting them together, as shown in Figure 6, although we actually performed a vast number of numerical simulations. Figure 6 well presents the sensitivity of each controlling parameter. Obviously, a large correlation distance corresponds to weak reflections; in addition, a large RMS perturbation corresponds to strong reflections. Therefore, we can fit the Yutu-2 LPR data by searching for different correlation distances and RMS perturbations to obtain the subsurface features of the lunar regolith at the Chang'E-4 landing site.
vast number of numerical simulations. Figure 6 well presents the sensitivity of each controlling parameter. Obviously, a large correlation distance corresponds to weak reflections; in addition, a large RMS perturbation corresponds to strong reflections. Therefore, we can fit the Yutu-2 LPR data by searching for different correlation distances and RMS perturbations to obtain the subsurface features of the lunar regolith at the Chang'E-4 landing site.

The Diameter and Stacking Density of Rocky Blocks
A 4 m thick ejecta layer at the bottom was used to qualitatively identify the rock size and spatial density in the shallow surface of ejecta. The physical properties of the lunar regolith were further determined by the relative amplitude of weak reflections within the lunar regolith compared with the strong reflections from the ejecta layer. We set up three kinds of models: without rocky blocks (Figure 7a), with several~20 cm rocks (Figure 7b), and with a~1 m rock (Figure 7c). The simulated LPR profiles of these three models are shown in Figure 8. regolith were further determined by the relative amplitude of weak reflections within the lunar regolith compared with the strong reflections from the ejecta layer. We set up three kinds of models: without rocky blocks (Figure 7a), with several ~20 cm rocks (Figure 7b), and with a ~1 m rock (Figure 7c). The simulated LPR profiles of these three models are shown in Figure 8.
The numerical simulation results show that the purely random model does not indicate an obvious interface between the lunar regolith and the ejecta layer, and the reflection patterns in different time windows are similar (Figure 8a). After adding the rocky blocks at the bottom, an obvious interface appears (Figure 8b). In contrast, a large rock corresponds to a wide range of reflections in the shape of a hyperbola with a slight curvature (Figure 8c). Therefore, a single big rock (e.g., ≥1 m) or without rocky blocks does not match the LPR observations. There must be a stack of rocks at the bottom of the model. We further used different stacking densities (i.e., low, moderate, and high) of small rocks and constrained their diameters (Figure 9). Our results show that the rock diameters should range from ~20 to ~30 cm, and the stacking density of rocks should be about 1-2 within each square meter, as shown in Figures 9 and 10.   The numerical simulation results show that the purely random model does not indicate an obvious interface between the lunar regolith and the ejecta layer, and the reflection patterns in different time windows are similar (Figure 8a). After adding the rocky blocks at the bottom, an obvious interface appears (Figure 8b). In contrast, a large rock corresponds to a wide range of reflections in the shape of a hyperbola with a slight curvature (Figure 8c). Therefore, a single big rock (e.g., ≥1 m) or without rocky blocks does not match the LPR observations. There must be a stack of rocks at the bottom of the model. We further used different stacking densities (i.e., low, moderate, and high) of small rocks and constrained their diameters (Figure 9). Our results show that the rock diameters should range from~20 to~30 cm, and the stacking density of rocks should be about 1-2 within each square meter, as shown in Figures 9 and 10.

Results
To evaluate the structural characteristics of the regolith near the Chang'E-4 landing site, we quantitatively compared the Yutu-2 radar data with the simulated data. We adopted the envelope-stacking method [53], which discards phase information and only considers energy variations. This method has been successfully applied to seismograms for characterizing small-scale spatial heterogeneity of the mantle [54,55]. We took a small Figure 10. The simulated LPR profiles using the models shown in Figure 9. (a) With small spatial density of rocks; (b) with medium spatial density of rocks; (c) with large spatial density of rocks.

Results
To evaluate the structural characteristics of the regolith near the Chang'E-4 landing site, we quantitatively compared the Yutu-2 radar data with the simulated data. We adopted the envelope-stacking method [53], which discards phase information and only considers energy variations. This method has been successfully applied to seismograms for characterizing small-scale spatial heterogeneity of the mantle [54,55]. We took a small segment of the Yutu-2 LPR profile processed by Zhang et al. [24] as a reference. For both observed and synthetic data, we first calculated the envelope of each trace and then stacked them together. We also normalized the stacked envelopes and calculated the correlation coefficients as well as square errors between the observed and synthetic data (Tables 1 and 2). We only selected the data from 0.06 to 0.2 µs (~5-16 m) because the reflection patterns and energy of the first few meters could not be well considered by our numerical simulation method. More explanations on this can be found in the Discussion section. As shown in Figures 11 and 12, the stacked envelopes can well reflect the relative magnitude variation of reflections from the regolith and ejecta layers. As shown in Figures  13 and 14, when the correlation coefficient is large or the square error is small, the stacked envelope also fits well with the observed data ( Figure 11). Obviously, both the weak reflection of the lunar regolith and the strong reflection from the ejecta have good consistency in terms of radar profiles ( Figure 15) and stacked envelopes ( Figure 11). Therefore, we can determine the best group of the two controlling parameters: the correlation distance of the self-organization model is about 5-10 cm and the RMS perturbation is about 3%. In contrast, the other RMS perturbations (e.g., 1%, 5%, and 7%) show significant deviations from the stacked envelopes of the observed data. This indicates that the proposed method is sensitive to constraining the range of RMS perturbation. In addition, with the increasing value of RMS perturbation from 1% to 7%, the correlation coefficient changes accordingly and has a consistently varying trend. This indicates that the proposed method is feasible and robust for the evaluation of the RMS perturbation.
The correlation distance (about 5-10 cm) is reasonable because it is smaller than the spatial resolution (about 30 cm) of the high-frequency (500 MHz) channel of the LPR [56]; otherwise, we would observe a clear diffraction or reflection if the correlation distance was comparable to or bigger than 30 cm. The RMS perturbation of 3% is also reasonable for weak reflections of the lunar regolith, considering that a much higher perturbation or contrast would cause more evident reflections. Such a small RMS perturbation means that there is no significant variation of regolith materials; instead, local weak variations of structures due to an uneven deposit or gardening would be expected. The diffractions would be evident if the scattering body (with a much higher permittivity compared with that of the regolith) is close to the ground surface [16]; however, we did not observe clear diffractions within the whole regolith layer from 0 to 12 m, even around the upper surface. This indicates that the permittivity variation within the regolith should be relatively gentle. Therefore, we conclude that the RMS perturbation of 3% and the correlation distance (about 5-10 cm) well reflect the physical property of the regolith variations.      Figure 12 shows that a correlation distance of 5-10 cm is the best estimation, because that of 20 cm has a lower correlation coefficient and larger square errors compared with the observed LPR data. These results also indicate that the proposed method is not as sensitive to constraining the correlation distance within a range of 5 to 10 cm, at least compared with its high sensitivity to constraining the RMS perturbation. Considering the relatively small correlation distance and RMS perturbation in the regolith, we argue that this is more likely to be due to structural changes rather than material composition changes.
The size (or diameter) of rocks is estimated to be about 20-30 cm, which is reasonable because it is comparable with the spatial resolution (about 30 cm) of the LPR payload [56]. The stacking density of rocks is only about 1-2 per square meter, which is consistent with the expectation of the shallow structure of the ejecta [57]. Note that the total depth of this layer is up to 4 m (Figure 15). Such a thick layer with low spatial density of rocks may indicate that it experienced a completely different gardening history, because it is significantly different from the upper fine regolith layer and the typical rocky ejecta layer (at~20 m depth). It may be the transitional zone from ejecta to the regolith.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 19 Figure 13. Comparison of correlation coefficients between stacked envelopes of the observed and synthetic LPR data. These correlation coefficients are listed in Table 1. The correlation distance (about 5-10 cm) is reasonable because it is smaller than the spatial resolution (about 30 cm) of the high-frequency (500 MHz) channel of the LPR [56]; otherwise, we would observe a clear diffraction or reflection if the correlation distance was comparable to or bigger than 30 cm. The RMS perturbation of 3% is also reasonable for weak reflections of the lunar regolith, considering that a much higher perturbation or contrast would cause more evident reflections. Such a small RMS perturbation means that there is no significant variation of regolith materials; instead, local weak variations of structures due to an uneven deposit or gardening would be expected. The Figure 13. Comparison of correlation coefficients between stacked envelopes of the observed and synthetic LPR data. These correlation coefficients are listed in Table 1.  Table 2.  Table 2. Figure 14. Comparison of square errors between stacked envelopes of the Yutu-2 LPR data and the synthetic data. These square errors are listed in Table 2.  Figure 12 shows that a correlation distance of 5-10 cm is the best estimation, because that of 20 cm has a lower correlation coefficient and larger square errors compared with the observed LPR data. These results also indicate that the proposed method is not as sensitive to constraining the correlation distance within a range of 5 to 10 cm, at least compared with its high sensitivity to constraining the RMS perturbation. Considering the relatively small correlation distance and RMS perturbation in the regolith, we argue Our results show that the weak heterogeneity of lunar regolith is only about 3%, and its correlation distance is no more than 10 cm. Such a low degree of heterogeneity indicates that the lunar soil is mature in the Von Kármán crater. Consequently, it allows radar wave energy to penetrate more deeply, up to 45 m. In contrast, the lunar regolith is less mature in the Imbrium basin [27] compared with that in the Von Kármán crater. This kind of regolith would have a strong scattering effect, particularly close to the upper surface [16]. As a result, it would prevent radar waves from penetrating more deeply; thus, the LPR signals detected by the Chang'E-3 mission in the Imbrium basin are mainly concentrated within 5 m, with almost no clear deep reflections. This indicates that the regolith maturity may be reflected by the radar wave energy and the maximum penetrating depth. Therefore, in addition to the difference in chemical components, the maturity or heterogeneity of the regolith layer also contribute to the LPR reflections.

Discussion
The observed data have relatively strong reflections within the depth of~1-2 m on the top of the lunar regolith, but we could not explain them at the current stage because their simulation must consider the antenna effect and 3D models. Our code is currently able to process 2D random media but cannot well consider the antenna effect; in contrast, gprMax [58] can consider the antenna effect but it has significant difficulty in constructing random models [16]. Nevertheless, both the stacked envelope and reflection pattern are well matched between the synthetic and observed LPR data from~5 to 16 m, as shown in Figures 11 and 15. This verifies the feasibility of the proposed method.
As shown in Figures 11 and 12, the stacked envelopes of the synthetic data can well match those of the Yutu-2 radar data before 0.16 µs (corresponding to the bottom of the regolith layer), compared with those after 0.17 µs (corresponding to the top of ejecta layer). The amplitudes of reflections from rocks after 0.17 µs are also influenced by the rocks' position and depth; however, it is not necessary to better match the stacked envelopes of this part because we concentrate on the weak heterogeneity of the regolith layer, and the ejecta layer is only introduced to qualitatively evaluate the relative amplitude variations from the regolith layer to the ejecta layer. The local peak amplitudes of the stacked envelopes after 0.17 µs are generally comparable to each other; therefore, they can well present the general energy of reflections within the ejecta layer and may be a good reference to evaluate the general energy of reflections within the regolith layer.
In this paper, we only present 12 models with three correlation distances (i.e., 5, 10, and 20 cm) and four RMS perturbations (i.e., 1, 3, 5, and 7%). However, we actually tested a large number of various models to estimate the possible ranges of all parameters. Similarly, significant efforts were also made for estimating the rock diameter and spatial density. We used the 2D finite-difference time-domain method to simulate the LPR data [16,26]; unfortunately, even for a 2D model, and using an appropriate setting for the finite-difference method to avoid evident artifacts [16], processing would take about one day (~24 h) when running on a supercomputing cluster with 40 cores; that is, 3D numerical simulations are almost unaffordable at present, even though the radar waves are actually propagating through 3D random media.
The self-organization characteristics of the lunar regolith are only inferred by the Yutu-2 LPR data using numerical simulations, which may be corroborated using a drilling core. However, all existing drilling cores on the Moon are located on the near side, and the maximum length is no more than 3 m [18,59,60]. This indicates that our model cannot be verified from the perspective of geology until deeper drilling can be carried out on the far side of the Moon (to at least 5 m).

Conclusions
We used three self-organization random models with different autocorrelation functions to simulate the lunar regolith on the far side of the Moon. We also added randomly distributed rocks under the 12 m thick regolith to simulate the ejecta. Through a series of comparative analyses of the simulated and observed LPR data, we found that the structure of the lunar regolith in the landing area of Chang'E-4 can be expressed by an exponential autocorrelation random model. Numerical simulations were proven to be capable of quantitatively estimating statistical parameters of the regolith. The correlation distance was estimated to be about 5-10 cm and the RMS perturbation was about 3%. The ejecta under the lunar regolith is relatively sparse, only containing about 1-2 rocks per square meter; in addition, the diameter of the rocks is about 20-30 cm, which is comparable to but smaller than the spatial resolution of the high-frequency LPR. These results well constrain the ranges of the controlling parameters of self-organization models and provide a qualitative evaluation of the spatial density of the rocks on the top of ejecta layer. The relatively small correlation distance and RMS perturbation in the regolith indicate that the regolith is mature. The weak reflections within the regolith at the Chang'E-4 landing site are more likely to be due to local weak changes of structures rather than material composition changes.

Data Availability Statement:
The data used in this work is available on the Science and Application Center for Moon and Deep Space Exploration, Chinese Academy of Sciences (http://moon.bao.ac.cn accessed on 31 July 2021).