Application of Mathematical Morphological Filtering to Improve the Resolution of Chang ' e-3 Lunar Penetrating Radar Data

As one of the important scientific instruments of lunar exploration, the Lunar Penetrating Radar (LPR) onboard China’s Chang'E-3 (CE-3) provides a unique opportunity to image the lunar subsurface structure. Due to the low-frequency and high-frequency noises of the data, only a few geological structures are visible. In order to better improve the resolution of the data, band-pass filtering and empirical mode decomposition filtering (EMD) methods are usually used, but in this paper, we present a mathematical morphological filtering (MMF) method to reduce the noise. The MMF method uses two structural elements with different scales to extract certain scale-range information from the original signal, at the same time, the noise beyond the scale range of the two different structural elements is suppressed. The application on synthetic signals demonstrates that the morphological filtering method has a better performance in noise suppression compared with band-pass filtering and EMD methods. Then, we apply band-pass filtering, EMD, and MMF methods to the LPR data, and the MMF method also achieves a better result. Furthermore, according to the result by MMF method, three stratigraphic zones are revealed along the rover's route.


Introduction
Lunar Penetrating Radar (LPR) is one of the main payloads mounted on the "Yutu" rover of China's Chang'E-3 (CE-3) spacecraft [1,2].Its objective is to detect geological information of the lunar subsurface such as thickness and distribution of lunar regolith, and the geological structure in the shallow lunar crust.Unlike the Apollo Lunar Sounder Experiment (ALSE) [3] and Lunar Radar Sounder (LRS) [4], Chang'E-3 LPR directly collects data on the lunar surface, which has wider frequency band and higher spatial resolution.Lunar Penetrating Radar (LPR) is a kind of nanosecond imaging radar that is carrier free and operates in the time domain [1].There are two detection channels (Channel-1 and Channel-2) with center frequencies at 60 MHz and 500 MHz, respectively, in the LPR system.Channel-1 (CH1), with meter-level resolution, is used to detect structure of the shallow lunar crust, and Channel-2 (CH2), with a depth resolution of less than 30 cm, is used to map the thickness and structure of lunar regolith.
By processing and analyzing the CH1 and CH2 data, researchers have gained substantial geological information about the lunar subsurface.Su et al. [5] presented the LPR data processing and the preliminary results, which revealed the thickness of regolith with a range of about 4 to 6 m.More than nine subsurface layers had been identified by Xiao et al. [6], and the results suggest that the exploration region has experienced complex geological processes since the Imbrian.In order to analyze the near-surface stratigraphic structure of the regolith in detail, the 500 MHz data was processed, and four major stratigraphic zones were revealed along the exploration line [7].Based on the LPR data, FeO+TiO 2 content and the electrical parameters were estimated [8,9].According to the profile of permittivity, Feng et al. [10] considered that the depth of transition zone between the bedrock and regolith is about 2.5 to 3 m.The processing result of the LPR data from N201 to N203 indicates the thickness of the regolith and the locations and physical parameters of some rocks [2].These results were obtained based on the extraction and analysis of useful echo signals, which is a complicated procedure because the useful signals are often masked by noise.These noises may be the result of multiple factors such as the coupling between antennas, instrument system, electromagnetic interference, etc. [11].In order to improve the quality of radar image and obtain better information of the lunar subsurface structure, many filtering methods have been applied to the LPR data, such as deconvolution processes [12], background removal [10,13], bi-dimensional empirical mode decomposition filtering [2], amplitude compensation [8,10], band-pass filtering [5,11] etc.However, not all the results are satisfactory, so it is still meaningful to improve the resolution of the LPR data by a new filtering method.
Mathematical morphology was presented by Matheron and Serra [14]; it is a nonlinear signal processing method based on mathematical theories such as integral geometry and random set theory.Mathematical morphology has been successfully applied in image processing [15][16][17][18][19] and noise suppression [20][21][22][23].In the processing of geophysical exploration data, Wang et al. [24] first introduced morphological theory into seismic data processing by suppressing outlier impulse noise in seismic data, and the basic waveform of the effective signal and the characteristics of amplitude and frequency of the effective signal were well maintained.Recently, mathematical morphology was applied to fracture detection [25], linear coherent noise attenuation [26], seismic energy dispersion compensation [27], and low-frequency noise attenuation [28] in seismic data.Also, some researchers used it in weak signal detection of the microseismic monitoring [29,30] as well as in noise suppression of magnetotelluric sounding data [31,32].Although much attention has been attracted by mathematical morphology, it is still new in radar signal processing.
This paper first introduces the basic theory of mathematical morphology and the composition of combination morphological filter, then by using the different scale structural elements, the original morphological filter is decomposed into three-kind morphological filters with different scale range.The second scale-range morphological filter can be used to suppress the noise in a certain scale range without damaging the useful signals.The synthetic signal test shows that it has strong applicability and superiority, and is able to suppress the noise and extract the useful signal better than EMD and band-pass filter.Then, the morphological filtering method is applied to the processing of the CE-3's 500 MHz data, and also achieved good results, namely the lunar subsurface structure was clearly imaged.

Basic Theory
The basic idea of mathematical morphology is to quantitatively describe the geometric structure of the target signal by means of sets; that is, to use pre-defined structural elements to locally match or modify the geometric features of the signal, while retaining the main shape features of the target signal, so as to achieve the purpose of extracting useful information [31].Structural elements are the basic elements of mathematical morphology, which can have any shape and scale.
There are four basic operations in mathematical morphology: dilation, erosion, opening, and closing.Taking a one-dimensional discrete signal as an example, let the original signal f (n) be a discrete function defined on P = {0, 1, . . ., N − 1}, and the structural element g(n) be a discrete function defined on Q = {0, 1, . . ., M − 1}, and N >> M, then dilation operation and erosion operation can be defined as Equations ( 1) and (2).
where ⊕ and Θ represent dilation and erosion operations, respectively.Morphological dilation and erosion are essentially maximal and minimal operations within the definition domain of structural elements.The dilation operation represents an expansion process, which is used to fill the depression parts of the unsmooth boundary, and the algorithm can expand the peak and shrink the valley.Meanwhile, the erosion operation represents a contraction process, which is used to remove the bulge parts of the unsmooth boundary, and the algorithm has the ability to shrink the peak and widen the valley [31].Morphological closing and opening operations are formed on the basis of dilation and erosion operations.They are defined as follows: where symbols • and • represent open and closed operations, respectively.Opening operation can eliminate the spikes in the signal and suppress the positive impulse noise; closing operation can compensate the valley bottom and suppress the negative impulse noise [33].
In order to better filter the positive and negative impulse noise in the target signal, the morphological opening-closing and morphological closing-opening filters can be constructed by using the same structural element on the basis of closing and opening operations as follows [31,33]: However, due to the shrinkage of the opening operation, the output of opening-closing filter becomes smaller, and due to the expansion of closing operation, the output of closing-opening filter becomes larger [31,33].In order to achieve better filtering effect, the following filter is constructed through the average value of the above two filters as follows.
where M g f denotes the mathematical morphological filtering (MMF) result of the signal f when the structural element is g.

Selection of Structural Elements
The quality of morphological filtering depends not only on the selected morphological transformation, but also on the selected structural elements.The function of structural elements in morphological operations is similar to that of filter windows or reference templates in general signal processing.Both the shape and scale of structural elements affect the results of morphological operations.Some common shapes of the structural elements are linear, square, diamond, triangle, and sinusoidal.The design of structural elements usually depends on the shape of the signal to be processed.Because the sinusoidal shape is a close match with radar signal waveform, and the shape is simple, symmetrical, and smooth, a sinusoidal structural element is selected in this paper, which is defined as follows: Equation (8) shows that the parameters K and L control the amplitude and length of the structural element.According to the average maximum amplitude of the signal, K is relatively easy to determine.Then, L determines the scale of structural elements; as L increases, the scale of structural element will become larger.Generally speaking, a morphological filter with large scale structural element can extract large contour information of the target signal.A morphological filter with small scale structural element can extract relative detail information of the target signal.
In Figure 1a, Mf1, Mf2, Mf3, Mf4, Mf5, and Mf6 show the results of processing an original signal with Gauss noise using the MMF with structural elements g1, g2, g3, g4, g5, and g6.The K value (Equation ( 8)) of these six structural elements is 0.5, which is determined by the amplitude of the original signal, and their scales are shown in Figure 1b.It can be seen that as the scale of structural element increase, we can get larger contour information and less detail information from the original signal.In addition, we can know that the morphological filter (Equation ( 7)) with a single structural element cannot extract detailed information smaller than the scale of structural element, nor can they extract information within a certain scale range in the middle of the signal.In order to improve the applicability of the morphological filter shown in Equation ( 7) and extract different scale-range information from the target signal, we use the following idea to process the signal using the morphological filtering.
Equation 8 shows that the parameters K and L control the amplitude and length of the structural element.According to the average maximum amplitude of the signal, K is relatively easy to determine.Then, L determines the scale of structural elements; as L increases, the scale of structural element will become larger.Generally speaking, a morphological filter with large scale structural element can extract large contour information of the target signal.A morphological filter with small scale structural element can extract relative detail information of the target signal.
In Figure 1a, Mf1, Mf2, Mf3, Mf4, Mf5, and Mf6 show the results of processing an original signal with Gauss noise using the MMF with structural elements g1, g2, g3, g4, g5, and g6.The K value (Equation 8) of these six structural elements is 0.5, which is determined by the amplitude of the original signal, and their scales are shown in Figure 1b.It can be seen that as the scale of structural element increase, we can get larger contour information and less detail information from the original signal.In addition, we can know that the morphological filter (Equation 7) with a single structural element cannot extract detailed information smaller than the scale of structural element, nor can they extract information within a certain scale range in the middle of the signal.In order to improve the applicability of the morphological filter shown in Equation 7and extract different scale-range information from the target signal, we use the following idea to process the signal using the morphological filtering.When the morphological filter (Equation 7) with the structural element g1 is used to process the target signal f, the target signal f can be decomposed into two scale-range signals f1 and R1, as shown below.
where f1 denotes the MMF filtering result of the signal f when the structural element is g1, which represents the information larger than the scale of structural element g1, while R1 represents the information smaller than the scale of structural element g1.When the signal f1 is processed by the morphological filter (Equation 7) with the structural element g2, it can also be divided into two scalerange signals f2 and R2.
where the scale of g2 is larger than g1.f2 denotes the MMF filtering result of the signal f1 when the structural element is g2, which represents the information larger than the scale of g2, while R2 represents the information smaller than the scale of g2, but larger than the scale of g1.At this point, the original signal f can be made up of three scale-range signals R1, R2, and f2 as follows.
Similarly, when the signal fn-1 is processed by the morphological filter (Equation 7) with the structural element gn, it can also be divided into two scale-range signals fn and Rn.When the morphological filter (Equation ( 7)) with the structural element g 1 is used to process the target signal f, the target signal f can be decomposed into two scale-range signals f 1 and R 1 , as shown below.
where f 1 denotes the MMF filtering result of the signal f when the structural element is g 1 , which represents the information larger than the scale of structural element g 1 , while R 1 represents the information smaller than the scale of structural element g 1 .When the signal f 1 is processed by the morphological filter (Equation ( 7)) with the structural element g 2 , it can also be divided into two scale-range signals f 2 and R 2 .
where the scale of g 2 is larger than g 1 .f 2 denotes the MMF filtering result of the signal f 1 when the structural element is g 2 , which represents the information larger than the scale of g 2 , while R 2 represents the information smaller than the scale of g 2 , but larger than the scale of g 1 .At this point, the original signal f can be made up of three scale-range signals R 1 , R 2 , and f 2 as follows.
Similarly, when the signal f n−1 is processed by the morphological filter (Equation ( 7)) with the structural element g n , it can also be divided into two scale-range signals f n and R n .
where the scale of g n is larger than g n−1 , and n > 2. f n denotes the MMF filtering result of the signal f n−1 when the structural element is g n , which represents the information larger than the scale of structural element g n , whileR n represents the information smaller than the scale of g n , but larger than the scale of g n−1 .At this point, the original signal f can be composed of (n+1) scale-range signals, which can be expressed as follows: where scale (g 1 ) < scale (g 2 ) < scale (g 3 ) . . .< scale (g n−1 ) < scale (g n ).The first scale-range filter is similar to the high-pass filter, which can extract detailed information smaller than the scale of structural element g 1 .The (n+1) th scale-range filter is equivalent to low-pass filter, which can extract long-period information larger than the scale of structural element g n .The second scale-range filter is similar to the band-pass filter, which can extract certain scale-range information in the middle of the signal.However, these three filters differ from the frequency filtering due to the fact that the scale ranges do not correspond to the frequency bands.We can use different scale-range filters to suppress the noise beyond the scale range to separate the useful signal and the unnecessary information.The LPR CH2 data contains some unnecessary detailed information and some unnecessary rough outline information, and these two kinds of information are often removed by the band-pass filtering and EMD methods, so as to clearly show the subsurface lunar structure [11,13].However, the results are not attractive enough, so we try to use the second scale-range filter in Equation ( 16) to solve the problem.Therefore, it is important for us to find the appropriate two different scale structural elements of the second scale-range filter; in other words, we only need to find two different "L"s in Equation ( 8).According to our experience, the two L's values can be selected by our desired frequencies, since the relationship between L and frequency approximately obeys a function similar to power function, when the sampling frequency is constant.

Synthetic signal test
For this part, we applied the second scale-range morphological filter to the synthetic signal to test its effectiveness on noise suppression and useful information extraction, and made comparisons with band-pass filter and empirical mode decomposition (EMD) filter.The synthesized signal (signal 3) was composed of signal 1 and signal 2, as shown in Figure 2a.The signal 1 is the low-frequency and high-frequency noises, and signal 2 is the original signal with the main frequency of 500 MHZ.The amplitude spectra of signal 3 and signal 2 are shown in black and green curves in Figure 2b, respectively.It can be seen that there are obvious high-frequency and low-frequency interferences in signal 3 compared with the amplitude spectrum of signal 2.
Firstly, signal 3 was processed by morphological filtering method, the selected K value of the structural elements was 0.5, and the selected L's values of the structural elements were 7 and 10, respectively.Signal 4 is the processing result shown in Figure 2a and its amplitude spectrum is displayed in Figure 2b (red curve).It can be seen that the second scale-range morphological filter achieved a good result.In order to show the result of morphological filtering more clearly, we compared the normalized results of signal 2 and signal 4, respectively, as shown in Figure 2c.Although the morphological filtering result had some residual noise information as highlighted by the red arrow, they still matched signal 2 very well.
In order to show the superiority of morphological filtering in suppressing both high frequency and low frequency information, we also used band-pass filtering and EMD methods to process signal 3, respectively.Three band-pass filters (bp filter 1, bp filter 2, and bp filter 3) with different frequency ranges were used to process signal 3.The processing results and the corresponding amplitude spectrum are shown in Figure 3a,b, respectively.Compared with bp filter 1 and bp filter 2, bp filter 3 (300, 450, 600, 800 MHz) acquired a better result.However, there was still residual low-frequency and high-frequency information in the processing results, and compared with the original signal (signal 2), there were still significant differences.The results obtained by EMD method were not satisfactory either, because almost no good wavelet signal was in good agreement with the original signal.We chose Imf 1 as the filtering result, since Imf 1 contained the waveform information of the original signal, but its waveform and amplitude spectra still showed that it contains relatively much noise (Figure 4a,b).In order to quantitatively evaluate the results by different filtering methods, the S/Ns (signal-to-noise ratios) were calculated.The original S/N was -9.38 dB, and the S/Ns of the denoised results using bp filter 1, bp filter 2, bp filter 3, EMD method, and MMF were −7.99, −5.24, −0.65, −11.94, and 1.73 dB, respectively.It can be observed that the result by MMF had higher S/N.Therefore, compared with band-pass filtering and EMD, morphological filtering had more advantages and could better extract useful information from the signals.In order to show the superiority of morphological filtering in suppressing both high frequency and low frequency information, we also used band-pass filtering and EMD methods to process signal 3, respectively.Three band-pass filters (bp filter 1, bp filter 2, and bp filter 3) with different frequency ranges were used to process signal 3.The processing results and the corresponding amplitude spectrum are shown in Figure 3a,b, respectively.Compared with bp filter 1 and bp filter 2, bp filter 3 (300, 450, 600, 800 MHz) acquired a better result.However, there was still residual low-frequency and high-frequency information in the processing results, and compared with the original signal (signal 2), there were still significant differences.The results obtained by EMD method were not satisfactory either, because almost no good wavelet signal was in good agreement with the original signal.We chose Imf 1 as the filtering result, since Imf 1 contained the waveform information of the original signal, but its waveform and amplitude spectra still showed that it contains relatively much noise (Figure 4a,b).In order to quantitatively evaluate the results by different filtering methods, the S/Ns (signal-to-noise ratios) were calculated.The original S/N was -9.38 dB, and the S/Ns of the denoised results using bp filter 1, bp filter 2, bp filter 3, EMD method, and MMF were −7.99, −5.24, −0.65, −11.94, and 1.73 dB, respectively.It can be observed that the result by MMF had higher S/N.Therefore, compared with band-pass filtering and EMD, morphological filtering had more advantages and could better extract useful information from the signals.In order to show the superiority of morphological filtering in suppressing both high frequency and low frequency information, we also used band-pass filtering and EMD methods to process signal 3, respectively.Three band-pass filters (bp filter 1, bp filter 2, and bp filter 3) with different frequency ranges were used to process signal 3.The processing results and the corresponding amplitude spectrum are shown in Figure 3a,b, respectively.Compared with bp filter 1 and bp filter 2, bp filter 3 (300, 450, 600, 800 MHz) acquired a better result.However, there was still residual low-frequency and high-frequency information in the processing results, and compared with the original signal (signal 2), there were still significant differences.The results obtained by EMD method were not satisfactory either, because almost no good wavelet signal was in good agreement with the original signal.We chose Imf 1 as the filtering result, since Imf 1 contained the waveform information of the original signal, but its waveform and amplitude spectra still showed that it contains relatively much noise (Figure 4a,b).In order to quantitatively evaluate the results by different filtering methods, the S/Ns (signal-to-noise ratios) were calculated.The original S/N was -9.38 dB, and the S/Ns of the denoised results using bp filter 1, bp filter 2, bp filter 3, EMD method, and MMF were −7.99, −5.24, −0.65, −11.94, and 1.73 dB, respectively.It can be observed that the result by MMF had higher S/N.Therefore, compared with band-pass filtering and EMD, morphological filtering had more advantages and could better extract useful information from the signals.

Lunar Penetrating Radar data processing
For this part, we applied band-pass filtering, EMD, and morphological filtering methods to the LPR Channel-2 data, and made comparisons on their results.Channel-2 is installed on the bottom board of the lunar rover with a distance of about 300 mm from the lunar surface [1,12].Channel-2 has two receiving antennas, Channel-2A and Channel-2B, their working frequency ranges from 250 MHz to 750 MHz, and the main frequency is 500 MHz.The time window of collecting LPR data is 640 ns, and the sampling interval is 0.3125 ns.The route of lunar rover is shown in Figure 5.In this paper, we selected Channel-2B data of the route N104-N105 to study, and the distance of the N104-N105 route is 10 m.
The raw LPR data are archived and distributed by the National Astronomical Observatories, Chinese Academy of Sciences (NAOC) (http://moon.bao.ac.cn/).The radar data are divided into three levels including Level 0, Level 1, and Level 2. The initial data we chose in this paper was Level 2 data without band-pass filtering.Before using the band-pass filtering, EMD, and morphological filtering, some operations needed to be done due to the unique characteristics of LPR data acquisition.LPR collects data in unique ways, the motion of the Yutu rover is non-uniform, and the coordinates are recorded once a minute, thus there is a situation of repeated data acquisition at the same location.So, in order to avoid data discontinuity, it is necessary to eliminate the repeated data [12].According to the calibration tests that were performed by the Institute of Electronics, Chinese Academy of Sciences, the time delay of the radar echo signal transmitted from the lunar surface is 28.203 ns for Channel 2; so, a 28.203-ns time delay correction was needed [11].Radar wave propagates in the stratum and its amplitude decreased rapidly with the increase of depth, so we used the automatic gain control (AGC) method to balance the weak and strong amplitude of each trace data.In order to remove direct wave and highlight local anomaly information, we used a median filter to remove the background.Because most of the available information of the 500 MHz LPR data was concentrated in less than 100 nanoseconds, we removed the data after 100 ns.After the above data processing steps, the result and its amplitude spectrum are shown in Figure 6.It can be seen that there are obvious low-frequency (<150 MHz) and high-frequency (>600 MHz) noise in the data.Then, we applied bandpass filtering and EMD methods to the Channel-2B LPR profile data of the route N104-N105, respectively.The filtering results are shown in Figure 7.

Lunar Penetrating Radar Data Processing
For this part, we applied band-pass filtering, EMD, and morphological filtering methods to the LPR Channel-2 data, and made comparisons on their results.Channel-2 is installed on the bottom board of the lunar rover with a distance of about 300 mm from the lunar surface [1,12].Channel-2 has two receiving antennas, Channel-2A and Channel-2B, their working frequency ranges from 250 MHz to 750 MHz, and the main frequency is 500 MHz.The time window of collecting LPR data is 640 ns, and the sampling interval is 0.3125 ns.The route of lunar rover is shown in Figure 5.In this paper, we selected Channel-2B data of the route N104-N105 to study, and the distance of the N104-N105 route is 10 m.
The raw LPR data are archived and distributed by the National Astronomical Observatories, Chinese Academy of Sciences (NAOC) (http://moon.bao.ac.cn/).The radar data are divided into three levels including Level 0, Level 1, and Level 2. The initial data we chose in this paper was Level 2 data without band-pass filtering.Before using the band-pass filtering, EMD, and morphological filtering, some operations needed to be done due to the unique characteristics of LPR data acquisition.The LPR collects data in unique ways, the motion of the Yutu rover is non-uniform, and the coordinates are recorded once a minute, thus there is a situation of repeated data acquisition at the same location.So, in order to avoid data discontinuity, it is necessary to eliminate the repeated data [12].According to the calibration tests that were performed by the Institute of Electronics, Chinese Academy of Sciences, the time delay of the radar echo signal transmitted from the lunar surface is 28.203 ns for Channel 2; so, a 28.203-ns time delay correction was needed [11].Radar wave propagates in the stratum and its amplitude decreased rapidly with the increase of depth, so we used the automatic gain control (AGC) method to balance the weak and strong amplitude of each trace data.In order to remove direct wave and highlight local anomaly information, we used a median filter to remove the background.Because most of the available information of the 500 MHz LPR data was concentrated in less than 100 nanoseconds, we removed the data after 100 ns.After the above data processing steps, the result and its amplitude spectrum are shown in Figure 6.It can be seen that there are obvious low-frequency (<150 MHz) and high-frequency (>600 MHz) noise in the data.Then, we applied band-pass filtering and EMD methods to the Channel-2B LPR profile data of the route N104-N105, respectively.The filtering results are shown in Figure 7.The band-pass filtering is widely used to tackle low-frequency and high-frequency noises.The real LPR data is tackled by this method, and the result shows that this method has a good performance to reject the unnecessary information, and the radargram structure become much clearer.However, there was still some noise in the result, as highlighted by the red arrow and circle (Figure 7a,b).Figure 7c,d are the filtering results by EMD method.It can be seen from the amplitude spectrum that this method can also remove the noises effectively, but the radargram displays an unsatisfactory result because there are many strips of interference in radar images that the echo signals characteristics cannot be clearly presented in some locations.
Before morphological filtering, the most important step is to determine the K and L of the two different structural elements.The K value is easy to select by the amplitude of preprocessing resulthere it is 2.5-while L is usually determined by the amplitude spectrum of the processing result.In order to conveniently select the L of the structural element, here an empirical relationship between L and frequency of the LPR data approximately has been obtained, which can be defined as: Figure 8a shows the empirical Equation 17 acquired by the experimental data; however, it is worth noting that if the sampling frequency changes, the empirical relationship will change.Here, Equation 17 is only suitable for the LPR data with the sampling interval 0.3125 ns.From Figure 6b, the frequency band is between 150 and 600 MHz; therefore, when the required L values of the structural elements of the second scale-range filter are 3 and 7, the proposed morphological filtering is used to extract useful information.The two structural elements are shown in Figure 8b.The results of the data processing using the second scale-range morphological filter in Equation 16are shown in Figure 9a,b.It can be seen that the morphological filtering method demonstrates a better result.Not only were the low-frequency and high-frequency noises effectively suppressed, but also the structure  The band-pass filtering is widely used to tackle low-frequency and high-frequency noises.The real LPR data is tackled by this method, and the result shows that this method has a good performance to reject the unnecessary information, and the radargram structure become much clearer.However, there was still some noise in the result, as highlighted by the red arrow and circle (Figure 7a,b).Figure 7c,d are the filtering results by EMD method.It can be seen from the amplitude spectrum that this method can also remove the noises effectively, but the radargram displays an unsatisfactory result because there are many strips of interference in radar images that the echo signals characteristics cannot be clearly presented in some locations.
Before morphological filtering, the most important step is to determine the K and L of the two different structural elements.The K value is easy to select by the amplitude of preprocessing result-here it is 2.5-while L is usually determined by the amplitude spectrum of the processing result.In order to conveniently select the L of the structural element, here an empirical relationship between L and frequency of the LPR data approximately has been obtained, which can be defined as: Figure 8a shows the empirical Equation ( 17) acquired by the experimental data; however, it is worth noting that if the sampling frequency changes, the empirical relationship will change.Here, Equation ( 17) is only suitable for the LPR data with the sampling interval 0.3125 ns.From Figure 6b, the frequency band is between 150 and 600 MHz; therefore, when the required L values of the structural elements of the second scale-range filter are 3 and 7, the proposed morphological filtering is used to extract useful information.The two structural elements are shown in Figure 8b.The results of the data processing using the second scale-range morphological filter in Equation ( 16) are shown in Figure 9a,b.It can be seen that the morphological filtering method demonstrates a better result.Not only were the low-frequency and high-frequency noises effectively suppressed, but also the structure of the radargram and the echo signals characteristics were imaged clearer, which provides much useful information for the geological interpretation of the lunar subsurface.

Stratigraphic division
The geological structure of the subsurface of the Moon provides valuable information on lunar evolution, and LPR is an important tool to reveal the lunar subsurface structure.For this part, according to the characteristics of echo signals displayed by morphological filtering result, a threelayer structure is separated by dotted lines (Figure 10).Meanwhile, some hyperbolas are selected, because every hyperbola could represent a buried rock or a void, and we can use the hyperbolic fitting method presented by Feng et al. [10] to estimate the velocity and depth.Since the middle part has similar characteristics of the echo signals, here we marked the upper and lower hyperbolas, and their corresponding estimated velocities and depths are displayed in Table 1, which is consistent with the result in [10].Then, combining with the forming mechanism of the regolith in the study area, a further interpretation of the LPR data was revealed.
(1) The thickness of the first layer varies slightly.The estimated maximum thickness of the layer is about less than 1 m, which is consistent with the estimated regolith thickness on the ejecta [7,11].The reflection information of this layer is little and the signal amplitude is relatively weak.We interpret this layer as fine-grained regolith.Even though more layers can be recognized based on the samples in the core tube, their thicknesses are typically of the order several centimeters, which are much smaller than the LPR range resolution [7], therefore it is not easy to identify more layers solely based on LPR image.
(2) The second layer is a region with bright radar echoes with larger amplitude, which contains numerous chaotic, irregular layers and hyperbolic curves.From Table 1, we can see that the velocities in the upper part are higher than in the lower part, which may be due to the increasing rock content

Stratigraphic division
The geological structure of the subsurface of the Moon provides valuable information on lunar evolution, and LPR is an important tool to reveal the lunar subsurface structure.For this part, according to the characteristics of echo signals displayed by morphological filtering result, a threelayer structure is separated by dotted lines (Figure 10).Meanwhile, some hyperbolas are selected, because every hyperbola could represent a buried rock or a void, and we can use the hyperbolic fitting method presented by Feng et al. [10] to estimate the velocity and depth.Since the middle part has similar characteristics of the echo signals, here we marked the upper and lower hyperbolas, and their corresponding estimated velocities and depths are displayed in Table 1, which is consistent with the result in [10].Then, combining with the forming mechanism of the regolith in the study area, a further interpretation of the LPR data was revealed.
(1) The thickness of the first layer varies slightly.The estimated maximum thickness of the layer is about less than 1 m, which is consistent with the estimated regolith thickness on the ejecta [7,11].The reflection information of this layer is little and the signal amplitude is relatively weak.We interpret this layer as fine-grained regolith.Even though more layers can be recognized based on the samples in the core tube, their thicknesses are typically of the order several centimeters, which are much smaller than the LPR range resolution [7], therefore it is not easy to identify more layers solely based on LPR image.
(2) The second layer is a region with bright radar echoes with larger amplitude, which contains numerous chaotic, irregular layers and hyperbolic curves.From Table 1, we can see that the velocities in the upper part are higher than in the lower part, which may be due to the increasing rock content

Stratigraphic Division
The geological structure of the subsurface of the Moon provides valuable information on lunar evolution, and LPR is an important tool to reveal the lunar subsurface structure.For this part, according to the characteristics of echo signals displayed by morphological filtering result, a three-layer structure is separated by dotted lines (Figure 10).Meanwhile, some hyperbolas are selected, because every hyperbola could represent a buried rock or a void, and we can use the hyperbolic fitting method presented by Feng et al. [10] to estimate the velocity and depth.Since the middle part has similar characteristics of the echo signals, here we marked the upper and lower hyperbolas, and their corresponding estimated velocities and depths are displayed in Table 1, which is consistent with the result in [10].Then, combining with the forming mechanism of the regolith in the study area, a further interpretation of the LPR data was revealed.
(1) The thickness of the first layer varies slightly.The estimated maximum thickness of the layer is about less than 1 m, which is consistent with the estimated regolith thickness on the ejecta [7,11].The reflection information of this layer is little and the signal amplitude is relatively weak.We interpret this layer as fine-grained regolith.Even though more layers can be recognized based on the samples in the core tube, their thicknesses are typically of the order several centimeters, which are much smaller than the LPR range resolution [7], therefore it is not easy to identify more layers solely based on LPR image.
(2) The second layer is a region with bright radar echoes with larger amplitude, which contains numerous chaotic, irregular layers and hyperbolic curves.From Table 1, we can see that the velocities in the upper part are higher than in the lower part, which may be due to the increasing rock content in the regolith.This layer is interpreted as ejecta containing many fragmented rocks, and the rocks might be produced during the formation of the CE-3 crater.
(3) The third layer is masked by a lot of salt-and-pepper noises, as shown in Figure 6a, which might be the result from overgained control of the LPR data.Through the morphological filtering, this noise can be removed, and the acquired characteristic of the radar echo signal is consistent with the results in [7,11].The radar echo signal displayed in this layer is weak, and only a few hyperbolic curves can be seen in this layer.Also, based on the velocities of the hyperbolas P4 and P5, √ ε = c/v is used to calculate the dielectric permittivity, and their corresponding dielectric permittivities are 7.7 and 9.2, respectively, which are close to 8. Previous research on Apollo samples showed that the permittivity of solid basalt is about 8, so we consider that this layer is basalt base.
in the regolith.This layer is interpreted as ejecta containing many fragmented rocks, and the rocks might be produced during the formation of the CE-3 crater.
(3) The third layer is masked by a lot of salt-and-pepper noises, as shown in Figure 6a, which might be the result from overgained control of the LPR data.Through the morphological filtering, this noise can be removed, and the acquired characteristic of the radar echo signal is consistent with the results in [7,11].The radar echo signal displayed in this layer is weak, and only a few hyperbolic curves can be seen in this layer.Also, based on the velocities of the hyperbolas P4 and P5, / c v ε = is used to calculate the dielectric permittivity, and their corresponding dielectric permittivities are 7.7 and 9.2, respectively, which are close to 8. Previous research on Apollo samples showed that the permittivity of solid basalt is about 8, so we consider that this layer is basalt base.

Discussion
In the above sections, the synthetic signal and the LPR data tests demonstrate the advantages of the proposed morphological filtering method to reject the low-frequency and high-frequency noise, compared with band-pass filtering and EMD methods.The algorithm of morphological filtering is easier to understand and operate, as it uses two different structural elements to extract the certain scale-range information from the original signal, and the result obtained by morphological filtering method has higher signal-to-noise ratio.Although the traditional band-pass filtering method is able to preserve the information of a certain frequency range in the signal, there are still many residual low-frequency and high-frequency noises in the result.The EMD method can decompose a signal into a group of sub-signals adaptively, and we can select one sub-signal or reconstruct several decomposed signals into one signal as the filtering result.However, there is no uniform standard for sub-signal selection or component reconstruction, so the acquired result varies from person to person and is not compelling enough.Though the proposed morphological filtering method has better effectiveness in noise suppression, the L's selection of the proper structural elements is not easy but very crucial.It is still necessary to present a more convenient approach to confirm the L of structural element comparing with the proposed method of empirical relation, since the empirical relation will change with the data sampling frequency.
Based on the result obtained by the morphological filtering, a three-layer subsurface structure at the CE-3 landing site can be revealed; meanwhile, the EM propagation velocity in different layers of the subsurface is estimated.According to the depth of the hyperbola P1 and arrival time, the first layer shows a high velocity, and its corresponding dielectric permittivity is about 3, which is

Discussion
In the above sections, the synthetic signal and the LPR data tests demonstrate the advantages of the proposed morphological filtering method to reject the low-frequency and high-frequency noise, compared with band-pass filtering and EMD methods.The algorithm of morphological filtering is easier to understand and operate, as it uses two different structural elements to extract the certain scale-range information from the original signal, and the result obtained by morphological filtering method has higher signal-to-noise ratio.Although the traditional band-pass filtering method is able to preserve the information of a certain frequency range in the signal, there are still many residual low-frequency and high-frequency noises in the result.The EMD method can decompose a signal into a group of sub-signals adaptively, and we can select one sub-signal or reconstruct several decomposed signals into one signal as the filtering result.However, there is no uniform standard for sub-signal selection or component reconstruction, so the acquired result varies from person to person and is not compelling enough.Though the proposed morphological filtering method has better effectiveness in noise suppression, the L's selection of the proper structural elements is not easy but very crucial.It is still necessary to present a more convenient approach to confirm the L of structural element comparing with the proposed method of empirical relation, since the empirical relation will change with the data sampling frequency.
Based on the result obtained by the morphological filtering, a three-layer subsurface structure at the CE-3 landing site can be revealed; meanwhile, the EM propagation velocity in different layers of the subsurface is estimated.According to the depth of the hyperbola P1 and arrival time, the first layer shows a high velocity, and its corresponding dielectric permittivity is about 3, which is consistent with that of Apollo regolith samples.In the second layer, the velocities in the upper part are higher than the lower part, which may be due to the higher rock content in the lower area.There is no obvious reflective interface between the regolith and basalt, and almost no hyperbola can be observed in the third layer.We suggest that the velocity of the third layer approximates to the velocity of hyperbola P5.These velocities can be utilized to do the migration, so that more details of the subsurface structure can be revealed.

Conclusions
In order to better suppress low-frequency and high-frequency noises, we proposed a morphological filtering method in this paper.This method uses two structural elements with different scales to extract a certain scale-range information from the original signal.The synthetic signal experiment shows that this method can extract useful information by suppressing information beyond the scale range of the two different structural elements.Compared with band-pass filtering and EMD methods, this method has better performance in noise suppression.Finally, we applied band-pass filtering, EMD, and morphological filtering methods to LPR CH2 data.The morphological filtering method also achieves a better result, where the low-frequency and high-frequency information covering the useful signal are obviously rejected much more, and the resolution of radar data profile is significantly improved.In addition, according to the echo signal characteristics of morphological filtering result, three stratigraphic zones are depicted along the exploration route, which provides valuable information to understand the evolution of the lunar subsurface regolith.

Figure 1 .
Figure 1.(a) The results of processing a signal using the mathematical morphological filtering (MMF) with different structural elements; (b) the corresponding structural elements of the results in Figure 1a.

Figure 1 .
Figure 1.(a) The results of processing a signal using the mathematical morphological filtering (MMF) with different structural elements; (b) the corresponding structural elements of the results in Figure 1a.

Figure 2 .
Figure 2. (a) The signal 3 is the synthesized signal which is composed of signal 1 and signal 2, signal 4 is the processing result of signal 3 using MMF; (b) the amplitude spectra of signal 2, signal 3 and signal 4; (c) the normalized results of signal 2 and signal 4.

Figure 2 .
Figure 2. (a) The signal 3 is the synthesized signal which is composed of signal 1 and signal 2, signal 4 is the processing result of signal 3 using MMF; (b) the amplitude spectra of signal 2, signal 3 and signal 4; (c) the normalized results of signal 2 and signal 4.

Figure 2 .
Figure 2. (a) The signal 3 is the synthesized signal which is composed of signal 1 and signal 2, signal 4 is the processing result of signal 3 using MMF; (b) the amplitude spectra of signal 2, signal 3 and signal 4; (c) the normalized results of signal 2 and signal 4.

Figure 3 .Figure 3 .
Figure 3. (a) The denoised results of the synthetic signal using band-pass filters; (b) the corresponding amplitude spectra of the results in Figure 3a.

Figure 4 .
Figure 4. (a) The denoised results of the synthetic signal using empirical mode decomposition filtering (EMD); (b) the corresponding amplitude spectrum of Imf 1.

Figure 4 .
Figure 4. (a) The denoised results of the synthetic signal using empirical mode decomposition filtering (EMD); (b) the corresponding amplitude spectrum of Imf 1.

Figure 5 .
Figure 5.The route of Yutu lunar rover, the distance from N104 to N105 is the research route (10 m).

Figure 6 .Figure 5 .Figure 5 .
Figure 6.(a) The preprocessing result of the Channel-2B LPR data of the route N104-N105; (b) the amplitude spectrum of the preprocessing result in Figure 6a.

Figure 6 .
Figure 6.(a) The preprocessing result of the Channel-2B LPR data of the route N104-N105; (b) the amplitude spectrum of the preprocessing result in Figure 6a.

Figure 6 .
Figure 6.(a) The preprocessing result of the Channel-2B LPR data of the route N104-N105; (b) the amplitude spectrum of the preprocessing result in Figure 6a.

Figure 7 .
Figure 7. (a) The filtering result using band-pass filtering method; (b) the amplitude spectra of the results in Fig. 6a and 7a; (c) the filtering result using EMD method; (d) the amplitude spectra of the results in Fig. 6a and 7c.

Figure 7 .
Figure 7. (a) The filtering result using band-pass filtering method; (b) the amplitude spectra of the results in Figures 6a and 7a; (c) the filtering result using EMD method; (d) the amplitude spectra of the results in Figures 6a and 7c.
Figure 7. (a) The filtering result using band-pass filtering method; (b) the amplitude spectra of the results in Figures 6a and 7a; (c) the filtering result using EMD method; (d) the amplitude spectra of the results in Figures 6a and 7c.
and the echo signals characteristics were imaged clearer, which provides much useful information for the geological interpretation of the lunar subsurface.

Figure 8 .
Figure 8.(a) Experimental data for fitting Equation 17; (b) the two structural elements of the second scale-range filter in Equation 16.

Figure 9 .
Figure 9. (a) The result by morphological filtering; (b) the amplitude spectra of the results in Figures 6a and 9a.

Figure 8 .
Figure 8.(a) Experimental data for fitting Equation (17); (b) the two structural elements of the second scale-range filter in Equation (16).

Figure 8 .
Figure 8.(a) Experimental data for fitting Equation 17; (b) the two structural elements of the second scale-range filter in Equation 16.

Figure 9 .
Figure 9. (a) The result by morphological filtering; (b) the amplitude spectra of the results in Figures 6a and 9a.
Figure 9. (a) The result by morphological filtering; (b) the amplitude spectra of the results in Figures 6a and 9a.

Figure 9 .
Figure 9. (a) The result by morphological filtering; (b) the amplitude spectra of the results in Figures 6a and 9a.

Figure 10 .
Figure 10.Interpretation of the lunar penetrating radar (LPR) data from N104 to N105

Figure 10 .
Figure 10.Interpretation of the lunar penetrating radar (LPR) data from N104 to N105.

Table 1 .
The estimated velocities and depths of the hyperbolas

Table 1 .
The estimated velocities and depths of the hyperbolas.