Rock Location and Property Analysis of Lunar Regolith at Chang’E-4 Landing Site Based on Local Correlation and Semblance Analysis

The Lunar Penetrating Radar (LPR) onboard the Yutu-2 rover from China’s Chang’E-4 (CE-4) mission is used to probe the subsurface structure and the near-surface stratigraphic structure of the lunar regolith on the farside of the Moon. Structural analysis of regolith could provide abundant information on the formation and evolution of the Moon, in which the rock location and property analysis are the key procedures during the interpretation of LPR data. The subsurface velocity of electromagnetic waves is a vital parameter for stratigraphic division, rock location estimates, and calculating the rock properties in the interpretation of LPR data. In this paper, we propose a procedure that combines the regolith rock extraction technique based on local correlation between the two sets of LPR high-frequency channel data and the common offset semblance analysis to determine the velocity from LPR diffraction hyperbola. We consider the heterogeneity of the regolith and derive the relative permittivity distribution based on the rock extraction and semblance analysis. The numerical simulation results show that the procedure is able to obtain the high-precision position and properties of the rock. Furthermore, we apply this procedure to CE-4 LPR data and obtain preferable estimations of the rock locations and the properties of the lunar subsurface regolith.


Introduction
The surface of the Moon is covered by regolith, which records at least a 4-billion-year history of meteoroid impacts and implantation of solar wind. The thickness of the regolith varies from about 2 m beneath the lunar maria to more than 10 m beneath the lunar highland, respectively, which is correlated with the age of the lunar surface [1]. The estimation of regolith thickness usually uses indirect measurements including seismic experiments [2], microwave remote sensing [3], and impact crater morphology and frequency distribution of the crater diameters [4]. Unlike the Apollo Lunar Sounder Experiment [5] and Lunar Radar Sounder [6], Chang'E-3 conducted the first in situ Lunar Penetrating Radar (LPR) exploration with a wide frequency band and high spatial resolution in the Imbrium basin, showing the detailed structure of lunar regolith [7,8]. In addition, Chang'E-4 (CE-4) is conducting the first in situ exploration on the farside of the Moon in the Von Kármán crater at the South Pole-Aitken (SPA) basin [9,10].
The LPR aboard on the Yutu-2 Rover of CE-4 is a nanosecond imaging radar which is carrier-free and operates in the time domain [11], and has two channels (CH-1 and CH-2) with center frequencies at 60 and 500 MHz, respectively [10]. CH-1 is used to map the structure of the shallow lunar crust with meter-level resolution and the CH-2 is used to detect the structure of regolith with a depth resolution of 0.3 m, which contains one transmitting antenna and two receiving antennas of different offsets-namely, CH-2A and CH-2B [12]. The LPR measurements suggest an underestimation of the global lunar regolith thickness by other methods and reveal a vast volume from the last volcanic eruption [8]. Yutu-2 has obtained a 425 m LPR profile along the rover tracks in the first 16 lunar days, displaying clear and complex subsurface structures of the landing area, which would help us to reveal the fine structure of lunar regolith.
The subsurface velocity of electromagnetic waves is a vital parameter for stratigraphic division, rock location estimates, and calculating rock properties including relative permittivity, density, and content of FeO and TiO 2 during the interpretation of LPR data. Meanwhile, the location of rocks in regolith can cause diffractions in LPR data, and the vertex position of the diffractions can help estimate the possible location of rocks. This is different from the main reflections from the layers, which are useful in stratigraphic analysis. Analysis of the rock location and the properties of the lunar regolith can help to reveal the formation and evolution history of the landing site. Previous studies on geological stratification and parameter inversion of the lunar regolith were mainly based on the CH-2B data only [3,8,12]. For example, Feng et al. proposed a hyperbola fitting method for radar velocity analysis in the CH-2B radar-gram [12]. Dong et al. calculated the parameters of the regolith by relative reflection amplitudes [13]. Lai et al. acquired the radar velocity based on the two-way delay method [14]. Hu et al. proposed an adaptive rock extraction method based on local similarity constraints to achieve the rock location and quantitative analysis for regolith [15]. Dong et al. analyzed the regolith properties from the LPR data of Chang'E-4 based on the 3D velocity spectrum [16]. However, a complex subsurface structure and interference of noise always result in incomplete, interlaced, and amplitude-varying hyperbolas which cause the inaccurate velocity analysis of field LPR data [16], and the previous studies cannot reach a satisfied accuracy. In addition, the determination of the rock location or the vertex position of the diffractions is the key procedure to obtain a high-precision diffraction velocity analysis spectrum. Although Zhang et al. used two sets of CH-2A and CH-2B data to estimate the electrical parameters and the iron-titanium content of regolith [17], the subsurface velocity estimation is not as accurate as expected due to the limitation of their method.
In this paper, two sets of CH-2A and CH-2B data are used to estimate velocity and the relative permittivity; thus, the other properties of lunar regolith including density and content of FeO and TiO 2 can be calculated based on the relative permittivity. Firstly, we consider the measurement similarity [18,19] between the two sets of CH-2A and CH-2B data based on the local correlation, and develop an extraction method to determine the vertex position of the diffractions and estimate the rock location. Secondly, we apply the normalized velocity spectrum based on the semblance analysis to estimate the subsurface velocity of the lunar regolith [20]. Finally, we utilize the interpolation method to obtain the velocity of the subsurface structure and derive the distribution of the relative permittivity.

LPR Data Processing
The CE-4 LPR is a dual-frequency ground penetrating radar (GPR) system, operating at 60 MHz (low-frequency) with a frequency band of 40~80 MHz and 500 MHz (highfrequency) with frequency band of 250~750 MHz [11,21]. The CH-2 radar data were collected during the first 16 lunar days along the Yutu-2 s traverse of about 425 m, as shown in Figure 1, which is consistent with the results acquired by Lin et al. [22]. To reveal the near-surface structure of the regolith, we processed and interpreted the high-frequency LPR data (CH-2A and CH-2B) by excluding the effects of the electromagnetic coupling with the rover's metallic body. The CH-2 antenna is mounted at the bottom of the lunar rover about 0.3 m away from the ground, and the space between two adjacent antenna elements (one transmitting antenna and two receiving antennae) is about 0.16 m [23]. The receiving antenna for CH-2A data is closer to the transmitting antenna than that of CH-2B, and thus has a lower signal-to-noise ratio. Therefore, we developed a series of procedures to process the CH-2A and CH-2B data as follows: (1) data deleting; (2) time delay removal; (3) band-pass filter application; (4) removal of DC bias; (5) background removal; (6) mean spatial filter application. Then, we could obtain the high-resolution CH-2A and CH-2B data with 1958 samples with 0.3125 ns time intervals and 11,661 traces of 0.0365 m spatial intervals, as shown in Figure 2.

Local Correlation.
The processed LPR CH-2A and CH-2B data are always affected by noise significantly, and the strong coherence of the noise may cause local maximums [15]. Here, we To reveal the near-surface structure of the regolith, we processed and interpreted the high-frequency LPR data (CH-2A and CH-2B) by excluding the effects of the electromagnetic coupling with the rover's metallic body. The CH-2 antenna is mounted at the bottom of the lunar rover about 0.3 m away from the ground, and the space between two adjacent antenna elements (one transmitting antenna and two receiving antennae) is about 0.16 m [23]. The receiving antenna for CH-2A data is closer to the transmitting antenna than that of CH-2B, and thus has a lower signal-to-noise ratio. Therefore, we developed a series of procedures to process the CH-2A and CH-2B data as follows: (1) data deleting; (2) time delay removal; (3) band-pass filter application; (4) removal of DC bias; (5) background removal; (6) mean spatial filter application. Then, we could obtain the high-resolution CH-2A and CH-2B data with 1958 samples with 0.3125 ns time intervals and 11,661 traces of 0.0365 m spatial intervals, as shown in Figure 2. To reveal the near-surface structure of the regolith, we processed and interpreted the high-frequency LPR data (CH-2A and CH-2B) by excluding the effects of the electromagnetic coupling with the rover's metallic body. The CH-2 antenna is mounted at the bottom of the lunar rover about 0.3 m away from the ground, and the space between two adjacent antenna elements (one transmitting antenna and two receiving antennae) is about 0.16 m [23]. The receiving antenna for CH-2A data is closer to the transmitting antenna than that of CH-2B, and thus has a lower signal-to-noise ratio. Therefore, we developed a series of procedures to process the CH-2A and CH-2B data as follows: (1) data deleting; (2) time delay removal; (3) band-pass filter application; (4) removal of DC bias; (5) background removal; (6) mean spatial filter application. Then, we could obtain the high-resolution CH-2A and CH-2B data with 1958 samples with 0.3125 ns time intervals and 11,661 traces of 0.0365 m spatial intervals, as shown in Figure 2.

Local Correlation.
The processed LPR CH-2A and CH-2B data are always affected by noise significantly, and the strong coherence of the noise may cause local maximums [15]. Here, we

Local Correlation
The processed LPR CH-2A and CH-2B data are always affected by noise significantly, and the strong coherence of the noise may cause local maximums [15]. Here, we introduce local correlation [18] to determine the rock position in order to take advantage of both LPR CH-2A and CH-2B data simultaneously.
The global correlation coefficient between two discrete signals, a i and b i , can be defined as [18]: which can be represented as two least-squares inverse based on a linear algebra notation: where N is the length of a signal, a and b are vector notions for a i and b i . Localizing Equations (3) and (4) were used to add regularization to inversion. Using shaping regularization, the scalars can turn into vectors [18]-i.e., scalars γ 1 and γ 2 turn into vectors c 1 and c 2 , defined as: where λ 1 = a T a and λ 2 = b T b are scaling controls, A and B are two diagonal operators composed of the elements of a and b, and S is a shaping operator, such as Gaussian smoothing, with an adjustable radius. The componentwise of vectors c 1 and c 2 defines the local correlation as c, whose elements are given by The local correlation is a measure of the similarity between two signals [15,19].

Semblance Analysis
The high-frequency electromagnetic wave will be reflected, diffracted, and refracted simultaneously when propagating the interface with discontinuous dielectric properties, such as boulders, voids, and soil inhomogeneities [24]. The LPR CH-2A and CH-2B data represent the common offset profiles, and the responses analyzed in semblance analysis are diffraction hyperbola instead of reflection hyperbola in the common midpoint data. The normal moveout equation was reconfigured for diffraction trajectories in the common offset profile, thus the traveltime was approximated as: where x is the position along the common offset profile, x 0 is surface position vertically above the diffracting target, t(x 0 ) is the two-way traveltime to the target when the antenna is positioned at x 0 [25], and υ st is the stacking velocity that is similar to the root-mean-square velocity. Velocity analysis is a reliable way to estimate stacking velocity and the accuracy depends on the signal-to-noise ratio of data, which influences the quality of picking. The semblance analysis was performed along the rover's traverse, which corresponds to the diffraction hyperbola in the LPR profile. The calculated stacking velocity by semblance analysis could be transferred to interval velocity through Dix inversion [26] as: where υ interval,n is the interval velocity of n-th layer and υ rms is the root-mean-square velocity.
Semblance analysis provides a measure of the coherency of energy along trial trajectories defined by the substitution of trial pairs of υ st and t(x 0 ), which was conducted over an aperture of L traces and 2M + 1 samples as: where L is 11,661 in LPR data, i and j are time sample indices, k is a trace number and Q j,k is the amplitude of the j-th sample of the k-th trace [27]. The selection of M would not influence the peak location of contours for semblance analysis. A proper selection of M would gain a better tradeoff between resolution (traveltime and velocity) and noise compression. Our sensitivity analysis indicates that a value of M = 3 would give a reasonable resolution and less noise level, so here we set M = 3.

Property Calculation
The relative permittivity without considering the electric conductivity and magnetic permeability can be approximated by the propagating velocity as [8]: where c is the speed of light in vacuum, which is 0.3 m/ns. The relation between the relative permittivity and density of lunar regolith [28] is: where ρ is the bulk density with the unit of (g/cm 3 ). Laboratory measurements show that the loss tangent represents the ratio of the imaginary part of the relative permittivity to its real part, which depends on bulk density. The (FeO+TiO 2 ) abundance of lunar regolith could be estimated based on the loss tangent and bulk density, which is beyond the scope of this study.

Simulation Data Results
We verify the effectiveness of the proposed method in a synthetic model using the finite difference time domain (FDTD) method [29], as shown in Figure 3a. The synthetic model consists of a heterogeneous background and several anomalous rocks with random relative permittivities less than 5.0 [24]. The dominant frequency of a Ricker wavelet is 500 MHz, the time windows and the sampling interval are 120 and 0.02 ns, respectively, and the trace interval is 0.005 m. Based on the picking positions in Figure 3d, we could scan the stacking velocity within the local scope of the simulation data. We display the contours from the semblance analysis for the velocity ranges from 0.1 to 0.3 m/ns within a 4 ns time window in   Figure 3b,c, respectively. The double diffraction hyperbolas were generated due to the upper and bottom surfaces of the rocks [16]. We focused on Remote Sens. 2021, 13, 48 7 of 14 the upper diffraction hyperbola to reduce the error of velocity estimation. Furthermore, we demonstrate the results by applying the local correlation based on the frequency wavenumber (FK) fan filtered data in Figure 3d, where the red-cross indicates the vertex position of the diffractions. Dong et al. determined the vertex positions of the diffractions by directly applying a 3D velocity spectrum for the whole data [16], which would cause a relatively large error because the maximum value position of semblance analysis would be significantly influenced by the noise. The vertex position of the diffraction determinations based on local correlation would have higher accuracy than the direct 3D velocity spectrum. Accurate determination of vertex positions of the diffractions is the basis of high-precision semblance analysis for the velocity.
Based on the picking positions in Figure 3d, we could scan the stacking velocity within the local scope of the simulation data. We display the contours from the semblance analysis for the velocity ranges from 0.1 to 0.3 m/ns within a 4 ns time window in Figure 4. The vertex positions of the diffractions are shown in the upper right corner of the graph. The stacking velocity and the calculated depth of each rock are listed in Table 1   The scanning positions are well consistent with the actual positions of rocks except the two rocks in the shallow layer ( Figure 5). The inconsistency may be due to the high velocity in the shallow layers which results in flat diffraction hyperbola, thus the accuracy of semblance analysis for velocity is insufficient. Thus, the simulation shows that our proposed procedures based on local correlation and semblance analysis is effective to obtain the subsurface parameters.  The scanning positions are well consistent with the actual positions of rocks except the two rocks in the shallow layer ( Figure 5). The inconsistency may be due to the high velocity in the shallow layers which results in flat diffraction hyperbola, thus the accuracy of semblance analysis for velocity is insufficient. Thus, the simulation shows that our proposed procedures based on local correlation and semblance analysis is effective to obtain the subsurface parameters. The scanning positions are well consistent with the actual positions of rocks except the two rocks in the shallow layer ( Figure 5). The inconsistency may be due to the high velocity in the shallow layers which results in flat diffraction hyperbola, thus the accuracy of semblance analysis for velocity is insufficient. Thus, the simulation shows that our proposed procedures based on local correlation and semblance analysis is effective to obtain the subsurface parameters.

LPR Data Results
The amplitude of echoes in the first~150 ns (length N = 501 with 0.3125 ns time interval), with a depth of about 10 m for the LPR CH-2A and CH-2B data, is displayed in Figure 6a,b, respectively. The local correlation for the filtered CH-2A and CH-2B data is shown in Figure 6c, where the red-cross indicates the vertex position of the diffractions after selection.

LPR Data Results
The amplitude of echoes in the first ~150 ns (length N=501 with 0.3125 ns time interval), with a depth of about 10 m for the LPR CH-2A and CH-2B data, is displayed in Figure 6a,b, respectively. The local correlation for the filtered CH-2A and CH-2B data is shown in Figure 6c, where the red-cross indicates the vertex position of the diffractions after selection. The semblance analysis contour of diffraction hyperbolas for the LPR CH-2B data was based on the rock location extraction by applying the local correlation, and the results are shown in Figure 7. The red-cross in the contour maps represents the maximum value of the semblance analysis. The stacking velocity for the rocks was obtained by picking the maximum value of the semblance analysis, and the calculated depths are The semblance analysis contour of diffraction hyperbolas for the LPR CH-2B data was based on the rock location extraction by applying the local correlation, and the results are shown in Figure 7. The red-cross in the contour maps represents the maximum value of the semblance analysis. The stacking velocity for the rocks was obtained by picking the maximum value of the semblance analysis, and the calculated depths are listed in Table 2.    Meanwhile, we compared the diffraction hyperbola for No. 30 and No.34 rocks with the CH-2B data. As shown in Figure 8 (Figure 8b), respectively. The trend of the hyperbola shows great consistency with the CH-2B data, which indicates the semblance analysis procedure can obtain the stacking velocity with high precision.
We applied the spline image interpolation method based on the 40 selected rocks to obtain the stacking velocity of the substructures. The stacking velocity and the rootmean-square velocity approach to the same when the source-receiver offset approaches zero for the isotropic layered model. The interval velocity can be calculated using the Dix formula as shown in Figure 9a. Subsequently, the estimated relative permittivity of the subsurface structure along the Yutu-2 rover's traverse, as shown in Figure 9b, could be obtained based on the interval velocity distribution. The subsurface velocity and the relative permittivity are important parameters for calculating other physical properties.
Based on the relative permittivity of the subsurface structure without considering the echo power attenuation, we could acquire the subsurface properties of lunar regolith, which are important for geological interpretation. We applied the spline image interpolation method based on the 40 selected rocks to obtain the stacking velocity of the substructures. The stacking velocity and the root-mean-square velocity approach to the same when the source-receiver offset approaches zero for the isotropic layered model. The interval velocity can be calculated using the Dix formula as shown in Figure 9a. Subsequently, the estimated relative permittivity of the subsurface structure along the Yutu-2 rover's traverse, as shown in Figure 9b, could be obtained based on the interval velocity distribution. The subsurface velocity and the relative permittivity are important parameters for calculating other physical properties. Based on the relative permittivity of the subsurface structure without considering the echo power attenuation, we could acquire the subsurface properties of lunar regolith, which are important for geological interpretation.   We applied the spline image interpolation method based on the 40 selected rocks to obtain the stacking velocity of the substructures. The stacking velocity and the root-mean-square velocity approach to the same when the source-receiver offset approaches zero for the isotropic layered model. The interval velocity can be calculated using the Dix formula as shown in Figure 9a. Subsequently, the estimated relative permittivity of the subsurface structure along the Yutu-2 rover's traverse, as shown in Figure 9b, could be obtained based on the interval velocity distribution. The subsurface velocity and the relative permittivity are important parameters for calculating other physical properties. Based on the relative permittivity of the subsurface structure without considering the echo power attenuation, we could acquire the subsurface properties of lunar regolith, which are important for geological interpretation.

Discussion
The noise interference and incomplete hyperbola will produce errors in semblance analysis. Additionally, the region selection of the data will influence the determination of the vertex position of the diffractions, which will affect the accuracy of semblance analysis. The determination of rock locations (i.e., the vertex position of the diffractions) is the most important procedure to obtain reliable semblance analysis results. Our analysis indicates that the semblance analysis in the shallow layers of the simulation model always has a larger bias than that in the deep layers. One possible reason is that the flat diffraction hyperbola would influence the accuracy of semblance analysis. Thus, a high-precision semblance analysis method such as the weighted AB semblance [30] or new migration methods [31] should be developed to improve the accuracy in the future studies. Additionally, LPR data processing should consider the local attributes to separate the data with higher accuracies.
We display the forward result with small Gaussian noise (mean is 0 and variance is 0.005) added and the corresponding local correlation results in Figure 10 to evaluate the validity of the combined method with noise. The results show that small noise has little influence on the local correlation results, and the combined method is valid to derive accurate velocity, permittivity, and bulk density distribution.

Discussion
The noise interference and incomplete hyperbola will produce errors in semblance analysis. Additionally, the region selection of the data will influence the determination of the vertex position of the diffractions, which will affect the accuracy of semblance analysis. The determination of rock locations (i.e., the vertex position of the diffractions) is the most important procedure to obtain reliable semblance analysis results. Our analysis indicates that the semblance analysis in the shallow layers of the simulation model always has a larger bias than that in the deep layers. One possible reason is that the flat diffraction hyperbola would influence the accuracy of semblance analysis. Thus, a high-precision semblance analysis method such as the weighted AB semblance [30] or new migration methods [31] should be developed to improve the accuracy in the future studies. Additionally, LPR data processing should consider the local attributes to separate the data with higher accuracies.
We display the forward result with small Gaussian noise (mean is 0 and variance is 0.005) added and the corresponding local correlation results in Figure 10 to evaluate the validity of the combined method with noise. The results show that small noise has little influence on the local correlation results, and the combined method is valid to derive accurate velocity, permittivity, and bulk density distribution. The local correlation and semblance analysis were combined to estimate the velocity structure, which can provide a good estimation of abnormal rock location and relative permittivity. The future multiple-input multiple-output radar onboard the lander from China's Chang'E-5 (CE-5) mission can help the detection of the regolith structure and the position of the underlying abnormal rocks in the landing site. The combined method by using local correlation and semblance analysis could be helpful for further processing on the CE-5 radar data. Furthermore, the method can be used in the site selection of the International Lunar Research Station (ILRS) [32] and similar missions in the future.

Conclusions
We derived the similarity between filtered LPR CH-2A and CH-2B data based on the local correlation, which would help to reveal the rock location distribution in the lunar regolith. Subsequently, we applied the semblance analysis method for the diffraction hyperbola to obtain the subsurface velocity structure and the relative permittivity dis- The local correlation and semblance analysis were combined to estimate the velocity structure, which can provide a good estimation of abnormal rock location and relative permittivity. The future multiple-input multiple-output radar onboard the lander from China's Chang'E-5 (CE-5) mission can help the detection of the regolith structure and the position of the underlying abnormal rocks in the landing site. The combined method by using local correlation and semblance analysis could be helpful for further processing on the CE-5 radar data. Furthermore, the method can be used in the site selection of the International Lunar Research Station (ILRS) [32] and similar missions in the future.

Conclusions
We derived the similarity between filtered LPR CH-2A and CH-2B data based on the local correlation, which would help to reveal the rock location distribution in the lunar regolith. Subsequently, we applied the semblance analysis method for the diffraction hyperbola to obtain the subsurface velocity structure and the relative permittivity distribution. The other properties of lunar regolith including bulk density and FeO and TiO 2 abundance can be calculated based on the relative permittivity. A procedure that combined the local correlation and semblance analysis has an advantage over traditional methods [12][13][14][15][16] in terms of accuracy and robustness in the interpretation of the high-frequency LPR data.
The velocity structure and the properties of the subsurface were derived from diffraction hyperbolas in the high-frequency LPR data, which could provide a good structural constraint for understanding the formation and evolution of the lunar regolith at the landing site. The combined procedure is an excellent choice for estimating the subsurface velocity and the other properties of the lunar regolith. Compared with the methods by analyzing Apollo's samples [1] and by the microwave remote sensing techniques [3] to derive the physical parameters of lunar regolith, the application of LPR has great advantages both in terms of the range and the accuracy. Therefore, it is very important to develop parameter estimation methods based on the LPR data, especially its high-frequency data. The procedure combined with the local correlation and semblance analysis, proposed in this paper, is of great significance in the interpretation of LPR high-frequency data.