**2. Methods**

#### *2.1. 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:

$$\begin{aligned} S\_{i,j,k} &= \sum\_{j=1}^{N\_i} f(t\_{i,j,k}, \mathbf{x}\_j) \\ \mathbf{x} = 1, \dots \mathbf{n} \mathbf{t}; \quad j = 1, \dots n \mathbf{x}; \quad k = 1, \dots n \mathbf{v} \end{aligned} \tag{1}$$

where *f* is the LPR data in *t*-*x* domain; *Ni* denotes the selected horizontal computation region size at *i*th 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 *j*th point and the extreme point of the hyperbola; *ti*,*j*,*<sup>k</sup>* is the two-way time of the *j*th points of the hyperbola [21,22], which can be obtained using the formula below:

$$t\_{i,j,k} = \left(t\_i^2 + \frac{4x\_j^2}{v\_k^2}\right)^{1/2} \tag{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 e ffective 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:

$$C\_{i,j,k} = \frac{1}{N\_i L} \sum\_{l=i}^{L+i-1} \left| \frac{S\_{i,j,k}^2}{\sum\_{j=1}^N f^2(t\_{i,j,k}, x\_j)} \right| \tag{3}$$

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 <sup>C</sup>*i*,*j*,*<sup>k</sup>* is not a ffected by the varying amplitude. Importantly, 1/*Ni* is added in this form to compensate the energy di fferences caused by using di fferent windows *Ni* along longitudinal direction. <sup>C</sup>*i*,*j*,*<sup>k</sup>* is 3D data, whose local maximum values indicate the time, horizontal position, and velocity of the points with hyperbolas.

## *2.2. 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:

$$\varepsilon = (\mathfrak{c}/\mathfrak{v})^2 \tag{4}$$

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]:

$$
\rho = \log\_{1.919}(\varepsilon) \tag{5}
$$

$$
\tan \delta = 10^{0.440\rho - 2.943} \tag{6}
$$

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]:

$$\tan \delta = 10^{\left[0.038 \times (\omega \text{(TiO}\_2) + \omega (\text{FeO})) + 0.312 \rho - 3.26\right]} \tag{7}$$

Combining (6) and (7), we can obtain the TiO2 and FeO weight percentage of lunar regolith.
