**3. Results**

#### *3.1. 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,b, 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*,*<sup>j</sup>* 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.

$$\mathbf{C}\_{i,j} = \max\_{k} \mathbf{C}\_{i,j,k} \tag{8}$$

$$\mathbf{C}\_{i,j} = \begin{cases} \mathbf{C}\_{i,j} - \varepsilon & \mathbf{C}\_{i,j} > \varepsilon \\ 0 & \mathbf{C}\_{i,j} < \varepsilon \end{cases} \tag{9}$$

**Figure 1.** Model test results. (**a**) The model used in simulation; (**b**) the simulated GPR profile; (**c**) 3D velocity spectrum of (**b**); (**d**) maximum value along velocity axis with the soft threshold function processing.

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 1. The results show that the estimation errors of upper hyperbolas are about 1%, and the estimation errors of bottom hyperbolas are about 10%. The errors are in acceptable range; the relatively higher errors of bottom hyperbolas result from the effects of scatters; the low velocities of electromagnetic waves inside scatters reduce the estimated velocities.

**Figure 2.** Velocity spectrum contour slices of three hyperbolas. The black crossed lines indicate the positions of maximum value.


**Table 1.** The estimated results and errors of model test.

#### *3.2. 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.



**Figure 3.** The walking route of Yutu-2 rover for the first day and the second day on the Moon.

**Figure 4.** The CH-2B profile of lunar penetrating radar (LPR) data.

#### *3.3. 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.

**Figure 5.** 3D velocity spectrum.

In order to locate the positions of hyperbolas on GPR profile, the maximum values *Ci*,*<sup>j</sup>* 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.

**Figure 6.** Maximum value along velocity axis. (**a**) Before compensation of 1/*Ni*; (**b**) after compensation of 1/*Ni*.

**Figure 7.** Initial positioning of hyperbolas based on soft threshold.

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:


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.

**Figure 8.** Result of hyperbola selection. (**a**) Hyperbolas on radargram; (**b**) velocity spectrum slices of 30 hyperbolas—the black crossed lines indicate the positions of maximum values.

#### *3.4. 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.

**Figure 9.** Analysis of radargram.

## *3.5. 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.

**Figure 10.** Velocity scatter figure.

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:

$$w\_{int,n} = \sqrt{\frac{\upsilon\_{rms,n}^2 t\_n - \upsilon\_{rms,n-1}^2 t\_{n-1}}{t\_n - t\_{n-1}}} \tag{10}$$

where *vint*,*<sup>n</sup>* and *vrms*,*<sup>n</sup>* represent the interval velocity and root-mean-square velocity of *n*th layer respectively, *tn* is the travel time to the *n*th 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.

**Figure 12.** The subsurface relative permittivity structure.

**Figure 13.** The density structure of lunar regolith.

**Figure 14.** The TiO2 and FeO content of lunar regolith.

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 Figures 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.
