Abstract
The Lunar Penetrating Radar (LPR) is one of the important scientific payloads in China’s Chang’E-3 (CE-3) to image within 100 m below the lunar surface. The acquired LPR data is significant for the research of lunar geological structure. Based on the sedimentary mechanism of lunar regolith, the regolith contains many rocks with different sizes. These local anomalies appear as diffraction in LPR data, which reduces the data quality and limits the structural analysis of lunar regolith. According to the kinematics characteristics of rock caused diffraction, we transform these problems to a problem of steep dip decreasing. To reach this goal, we adopt a data preprocessing workflow to improve the quality of the radar image, firstly. Then, a dip filter based on adaptive f-x empirical mode decomposition (EMD) is proposed to extract the rocks in the regolith and the corresponding removed IMF map indicates the degree of rock enrichment and highlights regolith-basement interface. Both simulation and LPR CH-2 data present a great performance. Finally, according to the processed result, we locate the position of each rock and highlight the contact interface of regolith and the basement rock.
1. Introduction
Chang’E-3 landed at 340.4875°E, 44.1189°N on the Moon on 14 December 2013 on a new region that has not been explored before in the largest basin, the Mare Imbrium [1]. The dual-frequency Lunar Penetrating Radar aboard the Yutu Rover provides a unique opportunity to map the subsurface structure to a depth of several hundreds of metres from the low-frequency channel (CH-1, 60 MHz) and the near-surface stratigraphic structure of the regolith from the high-frequency channel (CH-2, 500 MHz). The LPR also provides an accurate detection result with high resolution from high-frequency observations [2].
LPR data processing and initial results are first presented by NAOC [3]. Initial analysis of the LPR observations, especially that from the CH-1, indicates that there are more than nine subsurface layers from the surface to a depth of ~360m [1]. The onboard Lunar Penetrating Radar conducted a 114-m-long profile, which measured a thickness of ~5 m of the lunar regolith layer and detected three underlying basalt units at depths of 195, 215, and 345 m. The radar measurements suggest an underestimation of the global lunar regolith thickness by other methods and reveal a vast volume from the last volcanic eruption [4]. Fa et al., Lai et al., and Zhang et al. speculated the near surface structure by processing the raw CH2 data [5–7]. Dong et al. and Zhang et al. calculated the parameters of regolith [8, 9].
The previous papers mainly studied on the geological stratification and parameter inversion. The quantity and location of rocks in regolith have not been researched. The quantity and location of rocks in regolith not only help to understand the evolution of regolith on the landing site, but also provide a priori information for the further CE-5 plan of regolith collection. The rocks in lunar regolith are the break point, causing diffractions in LPR data. The kinematics characteristics of these diffractions are quite different from those of main reflections, which is expressed as hyperbola [10, 11]. The vertex position of the diffractions can indicate the location of rocks. Therefore, we can suppress high wavenumber components of diffractions to locate the rock.
Huang et al. [12] use empirical mode decomposition (EMD) to prepare stable input for the Hilbert transform. The aim of EMD is stabilizing a nonstationary signal and decomposing the nonstationary signal into fast and slow oscillation components, called intrinsic mode functions (IMFs). 1D EMD can be an adaptive band-pass filter, dividing the dataset into the IMFs with different frequency range. Due to the property of EMD, Bekara and van der Baan [13] propose f-x EMD to attenuate the random and coherent noise. Cai et al. [14] propose the guideline of t-f-x EMD denoising. Chen et al. [15] add autoregressive (AR) model to f-x EMD to improve the applicable conditions and propose an EMD based dip filter. Chen et al. [16] apply randomized-order EMD to attenuate multiple reflections using randomized-order EMD. Chen et al. [17–19] introduce EMD into the framework of Seislet transform and low-rank approximation, which expanded the applications of EMD, such as signal enhancement, rank reduction, etc.
In this paper, a processing workflow is designed to image the LPR CH-2 data. Then, an adaptive f-x EMD, which can extract the rock in the regolith, is adopted. We apply the removed IMF map as a reference in highlighted regolith-basement interface. An integrated regolith model helps to verify that this method is suitable for LPR data. Finally, according to the processed result, we locate the position of each rock and highlight the contact interface of regolith and the basement rock.
2. Data Preprocessing
Yutu Rover was released by CE-3 on 14 December 2013. The Yutu Rover explored the surface and subsurface of the landing site in the northern part of the Mare Imbrium using its four main instruments: the Panoramic Camera, Lunar Penetrating Radar (LPR), Visible–Near Infrared Spectrometer (VNIS), and Active Particle-Induced X-ray Spectrometer (APXS). Its track extends to 114.8 metres near a young crater. Figure 1 demonstrates Yutu’s path on the Moon (reproduced from Bin Hu et al [20]).
The dual-frequency Lunar Penetrating Radar aboard the Yutu Rover provides a unique opportunity to map the subsurface structure to a depth of several hundreds of metres from the low-frequency channel (CH-1, 60 MHz) and the near-surface stratigraphic structure of the regolith from the high-frequency channel (CH-2, 500 MHz). CH-2 data is selected to determine the near-surface stratigraphic structure of the regolith. According to the acquisition parameters, the actual situation, and the data quality, the LPR CH-2 data preprocessing workflow is designed (Figure 2). A radar image with high resolution is accessible after data preprocessing.
Analysis of the stratified structure of regolith layer and the location of each anomalous rock will help us understand the geological information of the CE-3 landing site and carry on the future lunar exploration program. However, there are some problems of the LPR data (Figure 2).
(1) The structure of the regolith layer is not clear. The current consensus is that generally the thickness of the regolith is ~5 m in the mare areas but the average thickness is approximately 10-15 m in older highland regions [21]. Theoretically, a uniform fine-grained regolith over an intact basement rock would produce an intense reflection, but there is nothing in LPR data.
(2) The location of each rock is not easy to determine. And the number of the rocks for each area will not be easily quantified.
Based on the above problems, it is considered that further processing is required. Adaptive f-x EMD is a reasonable approach that will help us deal with the above problems.
3. Methods
3.1. 1D t-x Domain EMD
EMD can provide an empirical decomposition of a nonstationary signal. These decomposed subsignals are separated based on oscillation frequency and called IMFs. A stable IMF has a constant instantaneous frequency and narrow-band waveform, satisfying two conditions [12]: (1) the numbers of extremes and zero crossings in the data series either are equal or differ by one, and (2) at any point, the mean value of the envelope defined by the local maxima and the local minima is zero.
In 1D case, the signal is decomposed into several subsignals with different frequency range, which can be written aswhere are the input signal and are the decomposed IMFs. denotes the residual and N is the IMF number.
Figure 3 demonstrates the decomposition of a synthetic data series using EMD in t-x domain. The generated signal is the sum of two trigonometric functions (i.e., ). The frequency components of the signal are 0.25 and 0.08 Hz. From Figure 3, we see that EMD successfully separates the signal into two IMFs with different frequency. Due to the advantages of frequency separation, EMD has been adopted in random noise attenuation [22, 23]. The denoising principle is the frequency difference between noise and signal, which random noise mainly locates on the high frequency components. By removing these IMFs, random noise can be attenuated.
(a)
(b)
3.2. f-x Domain EMD Based Dip Filter
Bekara and van der Baan [13] adopt EMD in f-x domain to suppress random noise and steep dip coherent noise. They consider the noise energy is dominant in the high wavenumber portion in the f-x domain. The high wavenumber portion presents the fast oscillation of each frequency slice. Based on the separation between noise and signal, noise can be attenuated by simply removing IMF1 from noisy data. The detailed process is shown as follows:(1)Set the size of the time window.(2)Pick a time window and adopt 1D forward Fourier transform along the time direction.(3)Pick a frequency slice and separate it into real and imaginary parts.(4)Compute IMF1s for real and imaginary parts and subtract them to obtain the filtered parts, respectively.(5)Compose the filtered frequency slice.(6)Repeat (4)-(5) for each frequency slice.(7)Adopt 1D reverse Fourier transform along the time direction.(8)Repeat (2)-(7) for each time window.
The two advantages of f-x EMD is convenience and stability. The f-x EMD is data-driven f‐k filter and does not require the predefined muting zone in f‐k domain, which is easily embedded into field data processing. Moreover, unlike convolutional operator based denoising methods (such as f-x predictive filter), f-x EMD can deal with irregular spatial sampling dataset [13, 24]. For data acquisition of LPR, irregular spatial sampling is inevitable because of the complex terrain, finite time, and expensive cost. Therefore, f-x EMD is a promising tool in LPR data processing.
It should also be noted that the choice of the removed IMFs can be more than one. It is determined by dispersion of high wavenumber components and noise level. When the target noise is located in high wavenumber components or noise level is low, the number of removed IMFs can be small. Conversely, more IMFs should be removed.
Due to the advantage of wavenumber separation, f-x EMD can be used as dip filter [15]. The different IMFs present different dip angle range; i.e., high dip components locate in the low IMFs and low dip components locate in the high IMFs. If we divide the IMF set into several subsets, the dataset is separated by the dip angle. Therefore, we define the dip filter using f-x EMD as follows:where is the filtered frequency slice and is the ith separated IMF. is the ith dip subsets and m is the number of divided subsets. denotes the weighting factors.
Figure 4 demonstrates the results of three types of dip filter (high-pass, mid-pass, and low-pass) working on a plane wave model (reproduced by Bin Hu et al. [20]). In the synthetic data (Figure 4(a)), it contains three dip sets. After the three dip filters, different dip events (Figures 4(b)–4(d)) are well extracted. Table 1 shows the detailed parameters of the three filter. The LPR data are acquired in constant-offset way, whose source can be considered as wave plane. The kinematics characteristics of the main reflection events are similar to those of the terrain. The rocks in lunar regolith are the break point, causing diffractions in LPR data. The kinematics characteristics of these diffractions are hyperbola. In a word, the reflection events are low dip and smooth, whereas the diffraction events are high dip. The diffraction point extraction can be transformed into a problem of steep dip decreasing. Therefore, we should select a simple low-pass dip filter, with , . The key parameter is the number of removed IMFs .
(a)
(b)
(c)
(d)
3.3. Adaptive f-x Domain EMD Based Dip Filter
According to the kinematics of diffractions, the goal of diffraction point extraction is to attenuate as much steep dip components as possible while maintaining the horizontal components in hyperbolic apex. Because of the complexity of rock distribution in lunar regolith, the distribution of different dip components is various. Therefore, it is necessary to develop a method in selecting removed IMFs adaptatively.
In our method, we adopt f-k filtering result as a reference and gradually subtract smaller IMFs. After every subtraction, we calculate the correlation coefficient between subtracted IMF and reference data. The correlation coefficient is calculated as follows:where , denote the ith IMF and the reference data in f-x domain, respectively, and denotes the mean of input data. The IMF according to the maximum of correlation coefficient present the horizontal component; then the number of removed IMF can be defined as where denotes the number of IMF according to the maximum of correlation coefficient. The detailed steps are shown as follows:(1)Set the size of the time window and obtain rough denoised data by f-k filter.(2)Pick a time window and adopt 1D forward Fourier transform along the time direction.(3)Pick a frequency slice and separate it into real and imaginary parts.(4)Compute one IMF for real and imaginary parts and subtract them to obtain the filtered parts, respectively.(5)Compose the filtered frequency slice.(6)Repeat (4)-(5) for each frequency slice.(7)Calculate the correlation coefficient between subtracted IMF and rough denoised data in f-x domain.(8)Repeat (3)-(7) till finding the maximum of correlation coefficient.(9)Adopt 1D reverse Fourier transform along the time direction.(10)Repeat (2)-(7) for each time window.
Figure 5 shows the demonstration of the choice in removing IMFs. We calculated the correlation coefficients (Figure 5(b)) between the first 5 IMFs (Figure 4) and denoised data with f-k filter. We can see that correlation coefficient of IMF3 is largest; i.e., we should remove the first 2 IMFs to extract the horizontal components. We also see the correlation coefficient of IMF1 is larger than other IMFs. The reason is the high wavenumber components residual caused by f-k filter. Note that the muting zone can be defined casually. The dip band divided by f-x EMD is smooth and stable. The high dip components residuals or muting artefacts (red arrows) in reference data have small influence on the calculation of correlation coefficient. Compare with conventional f-x EMD based dip filter; the proposed method is suitable for more complex conditions to highlight horizontal components.
(a)
(b)
4. Regolith Simulations
In order to verify whether this method is suitable for LPR data, we build a complex model (Figure 6). This model considers many factors: random medium, undulating interface, and anomalous body. The modeling method is referenced from [25–27]. FDTD is applied for the simulation of the simple model [28]. According to the actual acquisition parameters of LPR [2], the simulated parameters are shown in Table 2. The forward result is obtained in Figure 6.
(a)
(b)
Figure 7 shows the results of the rocks extraction (left) and the map of corresponding removed IMFs number (right). Note that singularity should be taken out before using the factor to avoid local discontinuities in the extraction results. The comparison of f-k spectrum (Figure 8) demonstrates the effectiveness of the proposed method, and the high-wavenumber components of the simulated data are effectively suppressed. From Figure 9, we see that the dip components are separated, and the remaining horizontal components indicate the location of the rock. Moreover, the IMFs map can reveal the degree of rock enrichment to a certain extent. The reason is that the order of the horizontal components is positively related to the complexity of the dip component, which indicates the intensity of the rock.
(a)
(b)
(a)
(b)
(a)
(b)
Figure 9 shows the interpretation of the integrated model result. According to the adaptive f-x EMD result of the integrated model, we locate the position of each rock block (Figure 9(a)); except for several redundant rocks (shown by red dots) and missing rocks (shown by the green dot), the position of the rock block is matched with the original model (Figure 6(a)). The IMFs map on the one hand helps to divide the area according to the number of rock blocks. On the other hand, it also helped us to get the interface of regolith and bedrock, which is not clear in the radar image. The Chang'E-5 mission will drill and collect the regolith from the moon [29], which requires that there are no rocks below the drilling point; otherwise the drilling machine will be damaged. The adaptive f-x EMD method and its IMFs map can help us select the drilling point, as shown in Figure 8(b) (2m and 18m), which will greatly reduce the probability of machine damage.
5. Results and Discussion
It is proven that the adaptive f-x EMD is suitable for LPR data processing. The processing results are shown in Figure 10. And the corresponding f-k spectrums (Figure 11) also show the effectiveness of high-wavenumber components suppression. We can observe that there are some artefacts in the deep section of the preprocessed data (Figure 2). These artefacts have certain continuity and gentle dip angle, which will affect the recognition of the rock, so we mute the section of deep artefacts. The processing result shows the location of the rock (left) and the degree of the rock enrichment (right).
(a)
(b)
(a)
(b)
In order to further explain the processing results, we extract the rocks (Figure 12) based on the feature of the processed LPR data. After the extraction of each rock, the extracted rocks are superimposed with Figure 12(b). The red area has the largest number of rock blocks, the yellow area is the second, and the blue area is the least, which is consistent with the above simulation results. The dividing lines in Figure 12(b) also help us to stratify the regolith layer and divide the LPR data into three layers.
(a)
(b)
After the above processing, we give an interpretation according to the formation of the lunar regolith and the former results. The LPR data are divided into three layers (Figure 13).
(1) Bedrock (below ~9 m). This layer is basalt, which is the product of the last basalt covering.
(2) Deep regolith (between ~3 m and ~9 m). This regolith layer contains a lot of rock fragments, because, after the formation of paleo-regolith, the small meteorite can not penetrate the upper regolith layer, leading to a lower degree of weathering, a larger grain size, and more rock masses.
(3) Shallow regolith (up ~3 m). This regolith layer has a higher degree of weathering, regardless of the impact of the meteorites or other weathering factors, such as the effect of space particles and electromagnetic radiation, etc., which can all affect the shallow regolith layer, thus decreasing the number and size of the rock blocks.
6. Conclusions
The LPR equipped on Yutu Rover detected the lunar geological structure in the Northern Imbrium. A data preprocessing workflow is designed to solve some types of issues, such as repeated and waste traces and noise. The position of each rock and the contact interface of regolith is still difficult to recognize. On the basis of good preprocessing, we propose an adaptive f-x EMD to locate the rocks in the regolith. The corresponding removed IMF map can reveal the degree of rock enrichment to a certain extent and the contact interface of regolith and the basement rock, which can provide a preference in data interpretation.
Based on the processed LPR CH-2 data, we locate the position of each rock and highlight the contact interface of regolith and the basement rock. The LPR data are divided into three layers: bedrock (below ~9 m), deep regolith (between ~3 m and ~9 m), and shallow regolith (up ~3 m). The bedrock is basalt, which is the product of the last basalt covering. The deep regolith layer contains a lot of rock fragments with a lower degree of weathering. The shallow regolith layer has a higher degree of weathering with a few number of the rocks.
These results provide valuable information regarding our understanding of the modification of the lunar surface and the evolution of the regolith, and the results are also important as a reference for future lunar sample return missions.
Data Availability
Data are hosted by National Astronomical Observatories, Chinese Academy of Sciences. Please visit: http://moon.bao.ac.cn.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors acknowledge Madagascar software for the open-source code and the open-source data in http://moon.bao.ac.cn. This work was supported by the Major Projects of the National Science and Technology of China (Grant No. 2016ZX05026-002-003), the National Natural Science Foundation of China (41374108, 41574097, and 41504083), the National Key Research Development Program (2016YC06001104), and the Graduate Innovation Fund of Jilin University (2018090).