Properties Analysis of Lunar Regolith at Chang’E ‐ 4 Landing Site based on 3D Velocity Spectrum of Lunar Penetrating Radar

: The Chinese Chang’E ‐ 4 mission for moon exploration has been successfully completed. The Chang’E ‐ 4 probe achieved the first ‐ ever soft landing on the floor of Von Kármán crater (177.59°E, 45.46°S) of the South Pole ‐ Aitken (SPA) basin on January 3, 2019. Yutu ‐ 2 rover is mounted with several scientific instruments including a lunar penetrating radar (LPR), which is an effective instrument to detect the lunar subsurface structure. During the interpretation of LPR data, subsurface velocity of electromagnetic waves is a vital parameter necessary for stratigraphic division and computing other properties. However, the methods in previous research on Chang’E ‐ 3 cannot perform velocity analysis automatically and objectively. In this paper, the 3D velocity spectrum is applied to property analysis of LPR data from Chang’E ‐ 4. The result shows that 3D velocity spectrum can automatically search for hyperbolas; the maximum value at velocity axis with a soft threshold function can provide the horizontal position, two ‐ way reflected time and velocity of each hyperbola; the average maximum relative error of velocity is estimated to be 7.99%. Based on the estimated velocities of 30 hyperbolas, the structures of subsurface properties are obtained, including velocity, relative permittivity, density, and content of FeO and TiO 2 .


Introduction
The Chinese Chang'E-4 mission for moon exploration has been successfully completed. The Chang'E-4 probe achieved the first-ever soft landing on the floor of Von Kármán crater (177.59°E, 45.46°S) of the South Pole-Aitken (SPA) basin on January 3, 2019 [1][2][3]. SPA basin is the broadest basin on the Moon. This ancient basin was born from asteroid impacts 4 billion years ago, recording the evolutionary history of the far side of the Moon, and is of great significance for researching the internal materials and structures of the Moon [4,5,6]. Von Kármán crater is one of the primary craters in the SPA basin, with a diameter of 186 km. Recent studies have revealed that the ejecta from adjacent craters have various contributions to the subsurface material of Von Kármán crater, which results in the complex subsurface structure at the Chang'E-4 landing site [1,3,4].
In order to patrol and investigate the lunar surface, a Yutu-2 rover is carried by Chang'E-4 probe. Yutu-2 rover is mounted with several scientific instruments containing a lunar penetrating radar (LPR) which has been verified to be an effective device for detecting the lunar subsurface structure [7,8]. The LPR is equipped with two types of channels (CH-1 and CH-2), the center frequencies of which are 60 MHz and 500 MHz [9,10]. In addition, the CH-2 possesses one transmitting antenna and two receiving antennas of different offsets which are also known as CH-2A and CH-2B [10][11][12]. With different frequencies, CH-1 and CH-2 have different detecting missions. The objective of CH-1 is to detect the deep structure of Von Kármán crater [13]; CH-2 is to map the detail of near surface layers, leading to the property analysis of lunar regolith [14,15].
During the interpretation of LPR data, subsurface velocity of electromagnetic waves is a vital parameter necessary for stratigraphic division and computing other properties. However, the velocity analysis of field LPR data faces many difficulties. Complex subsurface structure and interference of noise always result in the incomplete, interlaced, and amplitude-varying hyperbolas. In previous studies of Chang ' [17]. However, these methods need humans to select the hyperbolas, which is highly subjective. Thus, in this article, we applied 3D normalized velocity spectrum to estimate the velocity [18]. This method can automatically and objectively select hyperbolas and analyze the velocities; the normalization processing can solve the error brought by different amplitudes of different positions on hyperbolas; subsequently, during the computation, we applied a variable horizontal computation window along longitudinal direction to satisfy the field situation that rock sizes increase vertically.
This paper is organized as follows. Section 2 introduces basic theory of 3D velocity spectrum and properties analysis. In section 3, firstly a model test is performed to verify the feasibility of 3D velocity spectrum; subsequently the method is applied to the CH-2B data analysis. Subsequently, based on the 3D velocity spectrum, we obtain the positions and velocities of the points where hyperbolas exist; then the property structures of lunar regolith are computed including velocity, relative permittivity, density, and content of FeO and TiO2. In section 4, we discuss the computation error of each hyperbola. Section 5 is the conclusion.

3D Velocity Spectrum
The velocity spectrum analysis method is originally based on common middle point (CMP) ground penetrating radar (GPR) data [19,20], but it can also be applied to common offset data processing such as LPR data processing. In order to estimate velocity, a stacked amplitude is computed as follows: where f is the LPR data in t-x domain; Ni denotes the selected horizontal computation region size at ith time; to satisfy the field situation that rock sizes increase vertically, we make Ni variations along longitudinal direction; nt, nx, and nv are the number of sampling points of each trace, number of traces, and number of velocities used in computation; xj represents the horizontal distance between the jth point and the extreme point of the hyperbola; ti,j,k is the two-way time of the jth points of the hyperbola [21,22], which can be obtained using the formula below:  (2) where ti is the two-way time of extreme point of the hyperbola; vk represents the velocity used in computation.
However, the method of (1) is not always effective in the field as the hyperbolas are incomplete, interlaced, with varying amplitude. To solve the error brought about by these situations, the normalization form of stacked amplitude is applied: where L is the time gate; the average in the time gate is adopted because the hyperbolas on radargram are not curves but regions with hyperbolic shapes. The normalized form can guarantee that the result Ci,j,k is not affected by the varying amplitude. Importantly, 1/Ni is added in this form to compensate the energy differences caused by using different windows Ni along longitudinal direction. Ci,j,k is 3D data, whose local maximum values indicate the time, horizontal position, and velocity of the points with hyperbolas.

Properties Computation
After we obtain the velocities from the 3D spectrum, the properties can be computed. Without considering magnetic permeability, the relative permittivity can be easily derived using the following formula: where c equals 0.3 m/ns.
According to the studies of lunar regolith samples of Apollo, there is a relation between the relative permittivity and density of lunar regolith [9,17]: The tanδ is the loss tangent which represents the ratio of the imaginary part of the dielectric constant to its real part. The density of lunar regolith is shown in Figure 14. Subsequently, previous studies show that the loss tangent is also related to the TiO2 and FeO weight percentage values [9,17]: Combining (6) and (7), we can obtain the TiO2 and FeO weight percentage of lunar regolith.

Model Test
Before the processing of LPR data, a model test was performed to verify the feasibility of 3D velocity spectrum. The relative permittivity model used in test and its simulated GPR profile are shown in Figure 1a and 1b, respectively. In the simulation, the frequency of electromagnetic wave is 100 MHz; the time window and the sampling interval are 100 ns and 0.32 ns; the trace interval is 0.1 m; the transmitting and receiving antennas are placed on the ground surface. The radargram shows there are three pairs of hyperbolas (H1, H2, and H3); the double hyperbolas generate from the upper and bottom surface of the scatters. Subsequently, the 3D velocity spectrum is computed, whose result is shown in Figure 1c.
In fact, for a point (ti,xj), there is only one true velocity; this true velocity corresponds to the maximum value along velocity axis; the maximum value can help to find the position of hyperbolas on (ti,xj) profile. Thus, the maximum values Ci,j along velocity axis are figured out using (8). Subsequently, in order to suppress the noise interference and obtain the accurate positions of hyperbolas, we apply a soft threshold function [23] as (9) shows. The value of ε is set to 0.3 and the result is shown in Figure 1d. The energy stacks in three rectangles indicate the positions of hyperbolas. The irregular energy stack in the ellipse is generated by the intersection of two hyperbolas. After the x position is located, the velocity spectrum slices can be selected out to estimate the velocities. The spectrum contour slices of the three hyperbolas are shown in Figure 2. From Figure 2, the velocities can be read out. The estimated results and errors are shown in Table

Lunar Penetrating Radar (LPR) Data Processing
The walking route of Yutu-2 rover for the first day and the second day on the Moon is shown in Figure 3. The length of the path is about 105 m. Subsequently, the multi-segment data are spliced to achieve the CH-2B profile. After a series of processing ( Table 2) the LPR profile is obtained, which is shown in Figure 4. Table 2. Details of data processing.

Processing Details
Traces amending Adjusting the longitudinal displacement of traces, based on the phase of a strong reflection event

Traces selecting
The rover might stop at some points on the way to collect other scientific data but LPR never stops measurement, resulting in repeated acquisition of multiple traces at the same location. We average the repeated traces.

Time lag adjustment
There is a 28 ns delay for the start time of the transmitting antenna compared to the receiving antenna. Useless data deleting Signals after 500 ns are removed since these signals are not reliable.

Attenuation compensation
Conducting automatic gain control (AGC) to make deep information more visible Background removal Removing the average data along the rover path.

Band-pass filtering
Adopting band-pass filtering to suppress the low-frequency and highfrequency noise.

3D Velocity Spectrum
Based on (1)-(3), the 3D velocity spectrum is computed; the range of values are limited to 0 to 1; the 3D result without compensation of 1/Ni is shown in Figure 5. In order to locate the positions of hyperbolas on GPR profile, the maximum values Ci,j along velocity axis are also figured out, which is shown in Figure 6a. Subsequently, 1/Ni is adopted to compensate the energy differences caused by using different windows Ni along longitudinal direction, the result is shown in Figure 6b. The result shows the energies in shallow regions are compensated.
Subsequently, the soft threshold function is also applied. The value of ε is set to 0.6. Subsequently, based on the result of soft threshold processing, the locations of local maximums are selected which are shown in Figure 7.  Based on Figure 7, we can search for the velocities in 3D velocity spectrum in Figure 5. However, not all the points in Figure 7 are useful, there are three reasons: (1) For a highly flat interface, the estimated velocity will reach to 0.3 m/ns, which is obviously wrong; the lunar regolith velocity is less than 0.2 m/ns [24], so we delete the points with velocities close to 0.3 m/ns.
(2) There may be a large stacked energy even if there are no hyperbolas but interlaced regions of several hyperbolas; for this situation, we should delete the non-hyperbolic points.
(3) The hyperbolas on radargram are not curves but regions with hyperbolic shapes, so there will be several excess points at one hyperbola, which should be deleted.
After filtering using the above method, the 30 preserved hyperbolas are noted in Figure 8a; the velocity spectrum slices of these hyperbolas are shown in Figure 8b.

Analysis of Radargram
For Figure 8a, the hyperbolas located inside of the red rectangle, the region of 320-500 ns cannot be analyzed, so we focus on the region of 0-320 ns. Based on the reflection energy, several distinct layers and special regions are noted on the radargram (Figure 9). Layer 1 is a surface layer; a distinct reflection energy difference can be seen at 25 ns. Layer 3 is the bottom layer; although its adjacent upper and lower areas show few reflections, the interface shows strong reflection energy. Layer 2 is a complex layer. The zone of 25-150 ns is relatively more uniform than the deeper zone; it contains few reflections indicating its high weathering degree. The zone of 150-300 ns is more complex possessing strong reflections; the complex signals are formed by the interlaced hyperbolas generated by the scattering of waves on rocks. Four obvious regions are selected out. The strong reflections within Regions 1 and 3 are distinct from the adjacent areas. The reflection energy in Region 2 is relatively lower than that of Regions 1 and 3, which means the material and horizontal structure changes between Regions 1 and 3. Region 4 also contains strong reflections, but its vertical size is shorter because the material and horizontal structure changes again at the position of 80 m; the material at about 250 ns in Region 3 does not extend to Region 4. Overall, the subsurface material and structure at the Chang'E-4 landing site vary both vertically and horizontally, so the simple horizontal stratigraphic division is not appropriate within the zone of 25-300 ns.

Properties Analysis
Subsequently, properties analysis is applied to LPR data. Firstly, we combine the positions and velocities of selected hyperbolas in Figure 8 and obtain a velocity scatter figure (Figure 10). The region size is the same with Figure 9, the color indicates the velocity of each point.
. Subsequently, the image interpolation is applied to obtain the subsurface velocity of electromagnetic wave structure, then the smoothing is performed to remove the effects of outliers. Finally, as the velocities calculated by hyperbolic fitting method are the root-mean-square velocity (RMS), the interval velocities should be derived using the Dix formula:  (10) where vint,n and vrms,n represent the interval velocity and root-mean-square velocity of nth layer respectively, tn is the travel time to the nth layer. By using (10), the subsurface interval velocity structure is computed, and the profile is shown in Figure 11. Subsequently, based on the velocities, the relative permittivity, density, and content of FeO and TiO2 can be derived using (4)-(7), the results are shown in Figures 12-14, respectively.    In Figure 10, the points with different velocities distribute in different layers: 0.15 m/ns seem to be a watershed; the points whose velocities are bigger than 0.15 m/ns are mainly distributed in 0-150 ns zone; the points whose velocities are smaller than 0.15 m/ns are mainly distributed in 150-320 ns zone; this indicates the clear property difference between these two parts which can also be seen from Figures 11-14. These phenomena also correspond to the characteristic of radargram that 0-150 ns zone is of high weathering degree with few rock reflections and 150-320 ns zone is of low weathering degree with strong reflections. Figures 11-14 show that the velocity descends as the depth increases, and the opposite varieties for velocity and relative permittivity are due to (4). In these figures, we can also find that there is a fluctuation at about 250 ns after position of 80 m, which also corresponds to the analysis in which the material and horizontal structure changes at a position of 80 m; the material about 250 ns in Region 3 does not extend to Region 4. Previous research shows that the relative permittivity of lunar basalt is about 8 [24], so based on the velocity and relative permittivity structure we obtained, we think the region we researched (0-320 ns) is mainly the lunar regolith. The depth coordinates in Figure 9-14 are determined using average root-mean-square velocity of all traces; so the thickness of the lunar regolith is more than 20 m. In addition, the TiO2 and FeO content in the bottom of Figure 14 is about 17.5 wt%, which is close to the shallow lava flow observed in [3,4], so we speculate that the surface of the basalt is within 25 m in depth.

Discussions
In this section, the estimation errors of two-way reflected time and velocity are discussed. Take the first slice in Figure 8b for instance, a 3D surface plot can be computed, which is shown in Figure  15a. The maximum value in this figure indicates the two-way reflected time and velocity of this hyperbola. However, some error will generate from the noise interference and incomplete hyperbola; the correct point is probably in a region; a flat at an amplitude of 0.8 is selected to split out this region. Two profiles, which are the red rectangle and black rectangle in Figure 15, are extracted out to compute the error intervals of velocity and reflected time, respectively. The two intersections of the blue curve and gray line around the maximum value are selected to compose the error limit. In Figure 15, the error ranges of velocity and two-way reflected time are 0.161-0.179 m/ns and 55.8-56.84 ns, respectively. So, the real velocity and time may be within these two ranges. A maximum relative error (MRE) can be adopted to quantitatively describe the error: (11) where v and t are the estimated values using 3D velocity spectrum, vran and tran represent the error ranges. So, the MREv and MREt in Figure 13 are 5.29% and 1.05%, respectively.
Subsequently, the errors of all hyperbolas are computed. The average MREs of time and velocity are 0.68% and 7.99%, respectively. It is clear that the error of two-way reflected time can be omitted, which is far less than the error of velocity. Furthermore, we computed the distributions of velocities and error ranges in horizontal and vertical directions which are shown in Figure 16.

Conclusions
This paper applies the 3D velocity spectrum to property analysis of lunar penetrating radar data from Chang'E-4. The result shows that 3D velocity spectrum can automatically search for hyperbolas; the maximum value along velocity axis with a soft threshold function can provide the horizontal position, two-way reflected time and velocity of each hyperbola. Based on the estimated velocities of 30 hyperbolas, the subsurface properties' structures are obtained, including velocity, relative permittivity, density, and content of FeO and TiO2. Combining these properties and the radargram, we can conclude that 0-12 m is the fine-grained regolith; 12-20 m is the coarse regolith of different sources; the thickness of the lunar regolith is more than 20 m; the surface of the shallow basalt is within 25 m in the depth. Importantly, the subsurface material and structure at the Chang'E-4 landing site vary both vertically and horizontally, so the simple horizontal stratigraphic division is not appropriate within the zone of 25-300 ns. Finally, the estimation error is discussed, the errors of horizontal position and two-way reflected time are small which can be omitted; the average maximum relative error of velocity is 7.99%. the distributions of velocities and error ranges in horizontal and vertical directions are obtained.