Stochastic hyperbola fitting, probabilistic inversion, reverse-time migration and clustering: A novel interpretation toolbox for in-situ planetary radar

Ground-penetrating radar (GPR) is becoming a mainstream tool in planetary exploration, and one of the few in-situ planetary geophysical methods. There are currently three missions (Perseverance, Tianwen-1, Chang’E-4) with GPR-equipped rovers


Introduction
Ground-penetrating radar (GPR) is a near-surface geophysical method with a uniquely diverse set of applications, ranging from glaciology (Schroeder et al., 2020) and environmental science (Sonkamble and Chandra, 2021; Giannakis et al., 2019b) to landmine detection (Daniels, 2005;Giannakis et al., 2016) and non-destructive testing (Giannakis et al., 2021a). GPR is becoming a mainstream tool for planetary exploration -orbiter sounders have been employed both for Martian (Orosei et al., 2018) and Lunar (Kaku et al., 2017) exploration since the early 1970s, and during the last decade in-situ GPR has become a standard part of the scientific payload of planetary rovers (Lai et al., 2020). Giannakis et al. Of the eight active/planned rovers, five (Yutu-2 Ding et al., 2022, Perseverance Eide et al., 2021, Zhurong Zou et al., 2021, Rosalind Franklin Hervé et al., 2020, Chang'E-7 Zou et al., 2020 include GPR as part of their scientific instrumentation. Rover-coupled GPR was also used on the Yutu-1 rover and was part of the payload of the Chang'E-5 lander Su et al., 2022). GPR is, therefore, one of the most widely used scientific instruments on planetary rovers and space missions in general, and is pivotal for future planetary exploration. In particular, radar is critical for exploring Mars, as it allows to unveil in real-time the structure and properties underneath the moving spacecraft. Also, it may be used to assess regions of interest for drilling; interpret the surface layers and the timing of events; distinguish buried materials that may point to past depositional events; and even find reservoirs of water ice, which can be used for In-Situ Resource Utilization (ISRU).
In-situ planetary GPR was first used in 2013 as part of the scientific payload of Yutu-1, the planetary rover of the Chang'E-3 mission on the near side of the Moon (Su et al., 2014;Fa et al., 2015;Yuan et al., 2017;Feng et al., 2017;Lai et al., 2019;Ding et al., 2020). The same antenna systems, mounted on a similar rover named Yutu-2 were also employed in the Chang'E-4 mission in 2019 , the first human object ever landed on the far side of the Moon Tang et al., 2020). A year after, Perseverance was successfully landed at the Jezero crater on Mars, equipped with RIMFAX, a GPR system with 675 MHz central frequency (Hamran et al., 2020). On May 2021, the Chinese mission Tianwen-1 landed on Utopia Planitia on Mars. Zhurong, the rover of the mission, is equipped with two antennas, one low frequency with range 35-75 MHz, and a high frequency system with range 0.8-1.8 GHz (Zhou et al., 2016). Data from the Martian missions have recently become publicly accessible. Consequently, papers for these missions are scarce and focus primarily on preliminary results and observations (Casademont et al., 2022;Hamran et al., 2022). Due to that, the current paper focuses on the radar data from the Lunar Chang'E-4 mission, a well-researched area with numerous papers and various suggested processing and interpretation pipelines (Giannakis et al., 2021b).
Both Yutu-1 and 2 have two antennas with 500 MHz central frequency placed at the bottom of the rovers, and one 60 MHz antenna mounted at the back of the rovers . Due to the close proximity of the un-shielded 60 MHz antennas and the metallic parts of the Yutu rovers, the low frequency data from both missions are corrupted with ringing noise that makes interpretation unreliable (Li et al., 2018;Zhang et al., 2020). On the other hand, good quality data were collected from the 500 MHz antennas in both missions (Lai et al., 2019;Dong et al., 2020a), and especially in Chang'E-4 due to the low electromagnetic losses on the landing site (Dong et al., 2020a) (attributed to the low ilmenite content Giannakis et al., 2021b).
The processing applied to the Yutu-2 radagrams is a typical scheme that was initially developed for utility detection and non-destructive testing (Cassidy, 2009;Daniels, 2004). In particular, cross-coupling and ringing noise are removed using background removal, static noise is reduced using zero-offset filter, and the late reflections are enhanced using time-varying gain (Giannakis et al., 2021b). For the interpretation, linear inversion and conventional migration are applied to increase the overall signal to clutter ratio, and map the subsurface reflectors Li et al., 2020). For both inversion and migration, the electric permittivity is considered homogeneous, and is estimated using hyperbola fitting Li et al., 2020). More advanced approaches use multiple hyperbolas to approximate the 1D permittivity profile of the investigated medium (Dong et al., 2020a;Giannakis et al., 2021b). The estimated permittivity is then used to infer the mechanical and mineralogical properties of Lunar soils (Dong et al., 2020a,b) based on semi-empirical formulas derived from shallow samples collected during the Apollo missions (Olhoeft and Strangway, 1975;Carrier et al., 1991;Hickson et al., 2018). Full-Waveform Inversion -applied to extracted hyperbolas-also shown promising results in Earth applications (Jazayeri et al., 2018;Liu et al., 2018). Its applicability however is compromised by its computational requirements and the need for a realistic numerical equivalent of the antenna-rover system.
The current paper suggests a novel interpretation scheme capable of mapping subsurface targets and estimating the velocity structure of the host medium (Giannakis et al., 2022a). The scheme consists of a stochastic hyperbola fitting; a 1D probabilistic permittivity inversion; a reverse-time migration (RTM) using the finite-difference time-domain (FDTD) method (Leuschen and Plumb, 2001;Lu et al., 2016;Bradford, 2015;Bradford et al., 2018); and unsupervised clustering (Otsu, 1979). The novel stochastic hyperbola fitting addresses the non-unique nature of the problem, and instead of finding a set of parameters that fit the measured arrival times, it estimates a statistical range of parameters that equally fit the measured hyperbola. The statistical range of the bulk permittivities is subsequently used in a novel 1D probabilistic inversion that improves upon (Giannakis et al., 2021b) by estimating not only the permittivity profile but also its uncertainty as well. The estimated 1D permittivity profile is then used as an input to a bespoke RTM coupled with FDTD, which is no longer constrained to homogeneous velocity structures as in Zhang et al. (2020) and Li et al. (2020). This increases the overall accuracy of the migrated radargram, which is subsequently clustered (Otsu, 1979) into two classes, rocks/boulders and background Lunar soil. Due to lack of ground truth in real data from planetary radar, the proposed scheme is validated via a coherent numerical case study, indicating its effectiveness on estimating 1D permittivity profiles and locating distinct subsurface targets such as rocks/boulders. Lastly, the proposed processing pipeline is applied to GPR data from the first two Lunar days of the Chang'E-4 mission. The results reveal a layered structure and a detailed stochastic distribution of rocks/boulders within the top-most 10 m of the Lunar soil.

Methodology
The novel interpretation pipeline consists of four distinct and sequential steps: (A) stochastic hyperbola fitting; (B) probabilistic 1D velocity inversion; (C) RTM using FDTD; and (D) post-migration processing and clustering. Each of these steps will be described in detail in the following subsections. Numerical examples will also be given to showcase the effectiveness and applicability of the proposed scheme in cases that resemble Lunar soils where multiple targets are buried in layered media with gradational variations between layers.

Stochastic hyperbola fitting
Conventional hyperbola fitting assumes that GPR operates on top of a point target buried in a non-magnetic half-space with uniform permittivity distribution (Daniels, 2004). This configuration has been widely applied in both Chang'E-3 and E-4 missions Fa, 2020;Dong et al., 2020a,b;Lai et al., 2019) for estimating the bulk permittivity of the Lunar soil. From the estimated permittivity, the velocity of the medium is inferred, which is then used as an input in linear inversion or migration . Although there are automatic ways for picking hyperbolas (Lei et al., 2019), nonetheless, these approaches are applicable to clear hyperbolas in applications with high signal to clutter ratio such as concrete inspection and utility detection. In planetary radar, hyperbolas are not easy to pick, and an experienced GPR user is needed to reduce subjectivity. The manually picked arrival times are the raw data for the subsequent processing steps, and therefore special care should be taken to correctly pick and map them. Moreover, the selected hyperbola should be selected throughout the investigated range of depth to reduce the uncertainty in the final estimations.
As stated in Giannakis et al. (2021b), the underlying assumptions of hyperbola fitting constrain its applicability, especially in areas with I. Giannakis et al. Fig. 1. The generic scenario under consideration in stochastic hyperbola fitting. A common-offset GPR system operates on top of a spherical target with radius buried at depth inside a layered medium with varying permittivity with respect to depth ( ). In planetary environments, the conductivity ( ) is low, with negligible affects to the overall electromagnetic velocity (and consequently to the arrival times and hyperbola fitting).
non-uniform permittivity, and in the presence of relatively large targets. To address that, hyperbola fitting can be easily modified to include vertical variations of the permittivity (Giannakis et al., 2021b), and targets with arbitrary sizes (see Fig. 1). Despite these modifications, hyperbola fitting still suffers from non-uniqueness that greatly compromises the accuracy of the resulting permittivity estimation (Giannakis et al., 2022a). In this section we propose a novel scheme that takes non-uniqueness into account and calculates the probability density of permittivity and depth instead of specific values based on the point-target assumption (Giannakis et al., 2022a).
From Fig. 1, for a small offset, the arrival times of the reflections from a spherical target, with relatively small radius, buried in a half-space with smooth vertical permittivity variations ( ) can be approximated by: where is the arrival time (s); is the position of the centre of the bi-static antenna (m); 2 is the offset between the transmitter and the receiver (m); is the cover depth of the target (m); 0 is the projection of the target to the surface (m); is the square of the averaged squared root of the permittivity ( ) over ∈ [0, ] (Giannakis et al., 2021b); is the radius of the buried target (m), 0 = 2.99 ⋅ 10 8 m∕s is the velocity of light in free space; and and are the distances (m) from the centre of the target to the transmitter and receiver respectively. Eqs. (1)-(4) can sufficiently represent the arrival times of non-shallow point-like targets, regardless of the radiation pattern of the antenna system. For very shallow targets, the first arrivals originate from the surface waves. Therefore Eqs. (1)-(4) should be used with cautious for shallow targets. Moreover, in the presence of a high permittivity surface layer, significant refractions are expected between the free-space and the top-soil surface. In this case, hyperbola fitting should be modified to address this and incorporate the refractions in the calculations .
Hyperbola fitting can be expressed as an optimisation problem, in which the parameters { 0 , , , } are tuned in order to minimise the norm-2 error between number of measured = [ 1 , 2 , 3 ... ] ∈ R and simulated = [ 1 , 2 , 3 ... ] ∈ R arrival times: The minimisation in (5) is realised using a hybrid approach, where the results from particle swarm optimisation (PSO) (Kennedy and Eberhart, 1995) are used as initial points to a convex optimiser, in this case the simplex method. Via this approach we overcome potential local minima, and ensure that (5) will always converge. PSO and simplex method are arbitrarily chosen, similar results are achieved using any global optimiser (genetic algorithms, ant colony optimisation) coupled with any convex optimiser (non-linear least squares, gradient descent etc.). Notice that the apex of the hyperbola 0 is part of the unknown parameters to be tuned, since we investigate the whole profile instead of isolated sections. The minimisation in (5) has unique solutions only in clinical noise-free environments (Giannakis et al., 2022a,b). In the presence of noise the problem becomes ill posed with non-unique solutions, and a range of { 0 , , , } that gives rise to an equally good fit (Giannakis et al., 2019b(Giannakis et al., , 2022aMertens et al., 2016). This is clearly illustrated in Giannakis et al. (2022b) and also shown in Fig. 2, where (5) is applied to a typical scenario encountered in in-situ planetary radar. A target with = 0.2 m is buried at = 1.5 m in a homogeneous half-space with = 5. The arrival times of the simulated hyperbola were calculated using Eqs. (1)-(4) for homogeneous half-space. In clinical noise-free data, the minimisation in (5) always converges to the correct { 0 , , , }. In the presence of just 2% of Gaussian noise, the problem becomes unstable. Executing (5) 120 times, with different Gaussian noise each time, the results vary as illustrated in Fig. 2. This is more apparent in Fig. 3 where a noisy set of data from the scenario mentioned above, is fitted equally well with numerous hyperbolas, with permittivity ranging from = 4.6 − 7.
To overcome non-uniqueness, conventional hyperbola fitting as applied in planetary science Fa, 2020;Dong et al., 2020a,b;Lai et al., 2019), assumes a point target = 0. This is a simplification, since the rock/boulders on the Moon are expected to have varying sizes from = 0 − 5 m ( Bart and Melosh, 2010;Li et al., 2017). This can result in errors when estimating the bulk permittivity, as shown in Fig. 3, where the permittivity can vary from = 4.6 − 7 depending on the size of the target.
To overcome the inherent non-uniqueness of hyperbola fitting, we propose a novel technique that we call stochastic hyperbola fitting (SHF). In SHF, instead of calculating a specific set of { 0 , , } based on a given (conventional hyperbola fitting), we estimate the probability of { 0 , , , } to be true subject to the measured hyperbola, without making any assumptions regarding the radii of the targets. SHF first executes (5), and subsequently calculates the mean = [ − ] and the standard deviation = √ [( − − ]) 2 of the differences between the real ( ) and the synthetic ( ) observations. Then we assume that the measurements are corrupted with Gaussian noise, and that can be expressed as = + ( , 2 ). Based on that, we can derive a sufficient number of different sets of noisy measurements , and then execute the minimisation in (5) multiple times to estimate the range of { 0 , , , } that can satisfy (4) in the presence of different sets of noise  ( , 2 ). Regarding 0 , the uncertainty of the estimation is low due to the fact that 0 can be trivially calculated from the apex of the hyperbola. Furthermore, is irrelevant, since SHF is used to estimate the bulk permittivity at a given depth regardless of . Moreover, the subsequent processing steps are radius-agnostic and they do not require any prior information regarding the radii of the targets. Lastly, as stated in Giannakis et al. (2022b), the uncertainty of radius estimation using hyperbola fitting is high, making any estimation (even probabilistic) unreliable. For the reasons above, 0 and are omitted from the final output, which is now the multivariate kernel density estimation (KDE) (via diffusion Botev et al., 2010) with respect to and . Fig. 4 shows the KDE with respect to and using SHF on a hyperbola measured for = 0.2 m, = 1.5 m and = 6. The measurements are taken every 5 cm and the measuring line is 5 m long symmetrically around the centre of the target. The data are corrupted I. Giannakis et al.  with 1% Gaussian noise. The actual KDE is calculated by executing the minimisation in (5) a sufficient number of times using different noisy measurements with the same standard deviation. This is possible because both the noise-free hyperbola and the noise level are known. In real measurements we cannot separate the noise-free measurements and the noise, and therefore we can only estimate the KDE using SHF via the procedure described in the previous paragraphs. From Fig. 4 it is evident that SHF manages to approximate the multivariate KDE and provide a range of and , instead of specific uncertain values based on the clinical assumptions of conventional hyperbola fitting (i.e. = 0). The computational time needed for fitting one hyperbola using SHF is approximately 1-2 min using 1.6 GHz Intel Core i5.

Probabilistic 1D permittivity inversion
The main mechanism of deposition in Lunar soils is via cratering, which results in layered ejecta-blankets around impact craters (Melosh, 1989). Consequently, the proposed probabilistic inversion is based on the assumption that the permittivity in Lunar soils, for small scales, varies only with respect to depth (i.e. 1D model). Similar to any 1D geophysical method (vertical electrical sounding, 1D magnetotellurics etc.), if the investigated Lunar soil cannot be approximated by a 1D profile, the validity and reliability of the proposed scheme will be compromised.
SHF estimates the KDE with respect to and , which can be expressed as ( , | ) i.e. the conditional probability of ( , ) subject to the measured hyperbola . KDE can also be used to estimate the uncertainty for 0 and , nonetheless (as described in 2.1) these are omitted from the final KDE, primarily due to the fact that the subsequent processing steps only require knowledge regarding the bulk permittivity with respect to depth. The novel 1D probabilistic inversion utilises multiple KDEs from different targets in order to infer the 1D permittivity distribution in layered media, where the permittivity varies with respect to depth ( ). Similar to Giannakis et al. (2021b), the novel 1D probabilistic permittivity inversion can handle smooth velocity variations in contrast to conventional approaches based on Dix conversion (Giannakis et al., 2021b;Dix, 1955).  and . The parameters of the investigated case study are = 0.2 m, = 1.5 m and = 6. The arrival times are corrupted with 1% Gaussian noise. Up: Actual normalised KDE. Down: Estimated normalised KDE using SHF. SHF estimates the range of permittivities that can equally fit the measured data, subject to noise and the inherent non-uniqueness of the problem. Notice that in contrast to conventional hyperbola fitting, SHF makes no assumptions regarding the radius of the target.
As with any inversion scheme, before anything, we first need to parameterise our model, in this case the permittivity profile ( ). There are different ways to parameterise models prior to inversion, from typical grid-based discretisation to generic mathematical models. The novel 1D probabilistic inversion that we propose assumes that the permittivity can be described as a spline interpolation between equidistant depths. This model allows for smooth gradational layers; reduces the number of the parameters that need to be optimised; and greatly decreases the overall computational time.
Therefore, the first step of the probabilistic inversion is to discretise the permittivity in equidistant depths { 1 , 2 ... }. Spline interpolation is then applied to the discretised points resulting in a continuous and analytic representation of permittivity with respect to depth ( , { 1 , 2 ... }). A large allows for more complex profiles but increases the chances of over-fitting, while a small generates conservative low resolution models that can potentially under-fit the measurements.
The 1D probabilistic inversion tries to find the best set of { 1 , 2 ... } that minimises the error between the measured bulk permittivities using SHF and the bulk permittivity of ( , where { 1 , 2 ... } are the discretised permittivities used in spline interpolation; ( , { 1 , 2 ... }) is the analytic formula for the permittivity derived from spline interpolation; , and are random real numbers selected based on the KDE ( , , | ) estimated using SHF for the th hyperbola; and is the number of hyperbolas used in the inversion. In other words, we are trying to find the optimum set of { 1 , 2 ... }, which when used in spline interpolation gives rise to a continuous function ( , { 1 , 2 ... }) that minimises the squared error between the actual and the estimated bulk permittivities using SHF. Similar to (5), a hybrid optimisation approach is employed to solve (6), where the results from PSO (Kennedy and Eberhart, 1995) are used as initial points to the simplex method.
In (6) a random set of { , , } is used based on the KDE ( , , | ) calculated using SHF. Executing (6) multiple times, with a different set of { , , } each time (based on ( , , | )), will result in a set of permittivity profiles. All of them will equally satisfy our data, subject to the uncertainty due to noise and the inherent non-uniqueness of hyperbola fitting. Via this approach we get an insight into the uncertainty of the results and provide a range of acceptable solutions. This is in contrast to Giannakis et al. (2021b) where a single permittivity profile is derived based on the unrealistic assumption of conventional hyperbola fitting i.e. that the subsurface rocks/boulders are point targets ( = 0) .
A case study is examined next to investigate the accuracy of the proposed probabilistic inversion. The model is shown in Fig. 5. The open-source electromagnetic simulator gprMax  was used for the simulations. The dimensions of the model are 2 × 1 Fig. 5. The investigated scaled numerical case study. Ten targets are buried in a non-magnetic and non-conductive medium with varying permittivity with respect to depth. The shapes and sizes of the targets are not consistent, introducing uncertainty and non-uniqueness to hyperbola fitting. Fig. 6. The processed radagram from the case study shown in Fig. 5. Notice that the reflections from the layers are absent, due to the smooth variation of permittivity that decreases the reflection coefficients of the layers. m, with 2.5 mm spatial step. The time-step is 0.99 times the Courant limit (Taflove and Hagness, 2000). The host medium is non-magnetic and non-conductive with varying permittivity with respect to depth. Ten targets, nine cylinders and one box, are buried at different depths. The radii of the cylinders vary in order to add to the overall uncertainty of the estimated bulk permittivities. The model is excited by a monostatic Hertzian dipole with a Ricker pulse with 1 GHz central frequency . The measurements are taken every 2 cm from left to right (see Fig. 5). The resulting radagram is shown in Fig. 6. Only time-zero correction and time-varying gain were applied to the raw numerical data. Notice that there are no reflections from the horizontal layers due to the gradational transition between them (Diamanti et al., 2014). This means, that layers with smooth interfaces might go undetected using conventional reflection-based GPR interpretation techniques (Giannakis et al., 2021b). This is very important in planetary radar, since gradational layers are expected in ejecta blankets, especially at the transition between the weathered top soil and the underlying ejecta.

I. Giannakis et al.
The hyperbolas were fitted using SHF, and the KDE ( , , | ) was derived for every target. Subsequently, the 1D probabilistic inversion was executed 200 times, and each time different sets of { , } were chosen subject to their KDE ( , , | ). The 200 different 1D profiles were then used to estimate the range (95% confidence interval) of acceptable solutions, the results are shown in Fig. 7. Despite using hyperbolas from non-clinical targets with various shapes and sizes, the real permittivity profile is within the estimated range. The mean permittivity profile (red line in Fig. 7) is then used to calculate the bulk permittivity with respect to depth, which is then compared with the estimated distribution of bulk permittivities using SHF. From Fig. 8 it is apparent that the estimated permittivity profile ( ) results in a bulk permittivity variation that is aligned with the measurements, indicating the reliability of the estimated permittivity profile. The computational time needed for the investigated numerical study was approximately 5 min using 1.6 GHz Intel Core i5.

Reverse-time migration using FDTD
RTM coupled with FDTD can effectively focus the signal subject to any arbitrary permittivity distribution (Giannakis et al., 2020). This is an improvement compared to conventional migration that is typically applied in planetary radar, where the medium is assumed to be a homogeneous half-space; a clinical assumption that most often deviates from the truth. To our knowledge, there is no commercial software that  can perform RTM for GPR using numerical solvers such as FDTD. Due to that, we have developed our own bespoke FDTD solver (Taflove and Hagness, 2000) coupled with RTM as described in Leuschen and Plumb (2001). Similar to Giannakis et al. (2022a), a TM-FDTD with second order accuracy in both space and time is used in the current paper.
RTM requires the execution of FDTD, a time-consuming and complex algorithm with high computational requirements. Moreover, RTM uses 4 times the actual permittivity as input (to account for the twoway travel time in RTM), which leads to large numerical errors (Taflove and Hagness, 2000). In order to mitigate that, small discretisation steps need to be selected (Taflove and Hagness, 2000), which further I. Giannakis et al. Fig. 9. Applying Kirchhoff migration to the radagram shown in Fig. 6. The permittivity is assumed to be homogeneous and it is equal to = 3, 5, 7. Due to the layered structure of the medium (see Fig. 6), using low permittivity over-migrates deeper targets while using high permittivity under-migrates the shallow ones.
increases the overall computational requirements of RTM, making it unattainable for large radagrams measured in high permittivity media. Nonetheless, for Lunar applications, the investigated media are low dielectric targets (depending on their ilmenite content), with no high dielectric contrasts nor high permittivity regions (Dong et al., 2020a) that could compromise the accuracy and overall computational efficiency of RTM. Fig. 9 illustrates the migrated radagram -from here on denoted as ( , ) -of Fig. 6, using Kirchhoff migration with three different uniform permittivity values = 3, 5, 7. It is apparent that for < 5 the shallow targets are focused correctly while deeper targets are overmigrated. Increasing the permittivity results in under-migrating shallow targets while correctly migrating deeper ones. This is due to the smooth variation of permittivity with depth (as shown in Fig. 5) that deviates from the underlying assumption of conventional migration i.e. that the host medium is a homogeneous half-space. In Fig. 10, our bespoke RTM coupled with FDTD is applied to the same radagram (Fig. 6) subject to the 1D permittivity profile estimated using the probabilistic inversion (see Fig. 7). It is apparent, that using a good approximation of the actual permittivity profile, results in a more detailed migration, where all targets are equally focused regardless of their depth. The computational time for RTM varies greatly with the dimensions of the problem. In the absence of sufficient computational resources, large measuring lines should be divided and migrated separately.

Unsupervised clustering
The last stage of the novel processing pipeline consists of a series of post-migration steps followed by an un-supervised clustering. The main goal of this stage is to separate the rocks/boulders from the Lunar soil and map the rock distribution of the investigated area.
Firstly, the absolute value of the Hilbert transform is calculated for every trace, in order to extract its envelope ( , ) = |̂( , )| (Daniels, 2004), wherê( , ) is the Hilbert transform of ( , ). Subsequently, ( , ) is raised to the power ( , ) = ( , ) 1.5 to enhance the signal over migration artefacts and increase the compactness of the results. Raising ( , ) to a higher power could potentially suppress meaningful signal alongside artefacts and noise. Notice that in Lunar I. Giannakis et al. soil, the targets (rocks and boulders) are weak scatterers due to their low dielectric contrast with the background soil, resulting in weak secondary reflections that could otherwise give rise to artefacts in the migrated image. The processed image ( , ) is then convolved with a 2D Gaussian filter ( , ) = ( , ) * ( , ) to remove spikes and leftover artefacts from the Hilbert transform. Lastly, a threshold is selected that maximises the discriminant measure or the separability between the two classes as described in Otsu (1979). Based on this unsupervised threshold (Otsu, 1979), the image is clustered ( , ) into two classes i.e. targets and background medium. The clustering method described in Otsu (1979), otherwise known as Otsus' threshold, has never been reported in the GPR literature, but has been successfully applied in numerous two-class computer vision problems trying to separate the foreground from the background (Lee et al., 1990;Sezgin and Sankur, 2004). Otsu's method is equivalent to K-means clustering for image segmentation (Liu and Yu, 2009), and it produces a conservative threshold that minimises intra-class intensity variance . Fig. 10(A-E) illustrates the aforementioned steps, where the initial migrated image ( , ) (Fig. 10A) is transformed to the clustered ( , ) (Fig. 10E) . Every step in Fig. 10 is fully automatic and requires minimal computational resources. Comparing the results to the synthetic testcase problem (Fig. 5), it is apparent that the proposed processing framework can reliably estimate both the 1D velocity profile and the distribution of subsurface targets. Notice that although the presence and coordinates of the targets are correctly derived, nonetheless, the estimated sizes are not in good agreement with the actual sizes shown in Fig. 5. This is a known issue in migration/linear-inversion problems (Giannakis et al., 2021a(Giannakis et al., , 2019a. To accurately and reliably estimate the size of the investigated targets, more advanced approaches like FWI and machine learning should be considered (Giannakis et al., 2021a(Giannakis et al., , 2019a. To that extent, the proposed processing scheme, similar to all the previous approaches suggested for planetary radar, cannot be used for estimating the size of subsurface targets. .

Results from the Chang'E-4 radar data
The suggested interpretation pipeline is now applied to the data collected at the Von Kármán (VK) crater by the channel 2 (500 MHz) of the Yutu-2 rover. The underlying 1D assumption of the proposed scheme is not violated, since a plethora of evidence suggest that the investigated area is a layered medium (Huang et al., 2018;Li et al., 2020;Zhang et al., 2020;Dong et al., 2020a). The data consist of ∼ 106 m radagram collected over an irregular path during the first two Lunar days of the Chang'E-4 mission . Similar to Giannakis et al. (2021b), we focus on the first 150 ns of the scan, to find any hidden gradational layers within the Lunar regolith that were not detected using conventional GPR interpretation. Based on the lack I. Giannakis et al. of reflections within the first 150 ns, this layer was considered to be the regolith i.e. a relative homogeneous weathered layered. Nonetheless, as it is shown in Giannakis et al. (2021b) and in Fig. 6, layers with smooth gradational boundaries are transparent to electromagnetic waves, with practically zero reflection coefficients and no visible reflections. Consequently, lack of reflections does not exclude the presence of layers within the Lunar soil.
Typical processing was applied to the radargram prior to our scheme i.e. zero-time correction, dewow, exponential gain, and background removal Giannakis et al., 2021b). Left-over artefacts from background removal can corrupt the results especially in the shallow layers, where cross-coupling and ringing noise are most dominant. Nonetheless, that is not the case here as shown in Fig. 11. Subsequently, SHF was applied to 14 manually picked hyperbolas along the scan. The processed radagram and the selected hyperbolas are illustrated in Fig. 11. These hyperbolas were chosen because they are regularly distributed along the measuring line, and they extend throughout the investigated depth.
The 1D probabilistic inversion was executed using = 5 equidistant discretisation points from = 0−11 m. As discussed in Section 2.2, the value of dictates the complexity of the resulting permittivity profile. A large allows for more complex layered structures while a small results to simple conservative models. The choice of = 5 was done via a trial and error process where it was observed that increasing the discretisation points > 5 had minor effects on the estimated mean permittivity profile, while reducing < 5 led to simplified permittivity structures that under-fit the KDEs. The results are shown in Fig. 12.
In the presence of more hyperbolas the uncertainty range in Fig. 12 will be reduced. The KDEs ( , , | ) for the 14 hyperbolas, and the fitted bulk permittivity for the mean estimated profile (red line in Fig. 12) are shown in Fig. 13, which is also in good agreement with recently published data on the bulk permittivity at the landing site (Feng, Jianqing et al., 2022).
From Fig. 12, it is apparent that there is a low permittivity layer down to ≈ 3 m, followed by a relatively high permittivity layer until ≈ 6 − 7 m, followed by another low permittivity layer until ≈ 11 m. At ≈ 11 m, there are indications supporting the presence of another high permittivity layer, which is consistent with the strong radar reflections observed around 150 ns as shown in Li et al. (2020), Zhang et al. (2020), and remote sensing reflectance data that suggest the presence of a layer at around ≈ 8 − 13 m (Huang et al., 2018). The Fig. 13. The distribution of the bulk permittivities with respect to depth for the 14 hyperbolas shown in Fig. 11. The solid line is the bulk permittivity of the mean permittivity profile shown in Fig. 12. The scattered bulk permittivity is due to the ill-posed nature of hyperbola fitting, where in the presence of minimum noise it gives rise to non-unique solutions that satisfy equally well the measured hyperbola. permittivity profile shown in Fig. 12 is also in good agreement with Giannakis et al. (2021b) and Giannakis et al. (2022a) that estimated a similar layered structure in this area. The main difference between the current approach with Giannakis et al. (2021b), is the fact that the current processing pipeline takes into account the ill-posed nature of hyperbola fitting by mapping the uncertainty range, and plotting the plethora of permittivity profiles that can equally satisfy the measured hyperbolas.
The mean permittivity profile shown in Fig. 12, is then used as an input to our bespoke RTM coupled with TM-FDTD. The spatial step of FDTD grid was = = 1 cm, and the time-step was 0.99 times the Courant limit (Taflove and Hagness, 2000). The migrated image ( , ) is then transformed to ( , ) by taking the absolute value of the Hilbert transform of ( , ). Subsequently, ( , ) is raised to the power of 1.5 and convolved using a Gaussian filter resulting in ( , ). Lastly, ( , ) is clustered in two classes (rocks/boulders and host medium) using Otsus' method, a threshold that maximises the separability between the two classes (Otsu, 1979).
The post-processed migration processing steps and the final clustered image ( , ) are shown in Fig. 14. It is apparent that there is a stochastic distribution of rocks/boulders with some areas exhibiting very low rock abundance at ≈ 50 and ≈ 85 meters. Such areas could be potential targets for future drilling missions (similar to Chang'E-5 (Qian et al., 2021;Su et al., 2022)), where an optimised area with the least amount of rocks/boulders could be selected in order to reduce the energy requirements by avoiding drilling through potential large boulders. This approach can also be followed to detect buried discontinuous patches of water ice, permafrost or brines in the top meters of the soil.
Notice that the distribution of rocks follows a layered structure that is approximately aligned with the layered structure estimated by the 1D probabilistic inversion (Fig. 12). This is more apparent in Fig. 15. For every depth instance, all the pixels in Fig. 14 (E) classified as ''rock/boulder'' are summed and then normalised with respect to the depth with the maximum rock abundance. The layered structure of the rock distribution is aligned with the results from the 1D probabilistic inversion (Fig. 12), pointing to layered structure in the previously-assumed homogeneous regolith. Lunar rocks (primarily basalts, anorthosites and breccia) have higher permittivity than Lunar I. Giannakis et al. Fig. 14. The post migration processing pipeline as applied to the processed radagram for the first two Lunar days of the Chang'E-4 mission (Fig. 11). Green and blue colours in ( , ) (E) correspond to rock/boulders and host medium respectively. RTM coupled with TM-FDTD is used to migrate the data subject to the mean estimated permittivity profile shown in Fig. 12. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.) soil (Olhoeft and Strangway, 1975). Consequently, the presence of rocks within a layer can raise its overall bulk permittivity. Therefore, the changes in permittivity with depth as observed in Fig. 12 might be due to the changes in rock abundance rather than dielectric differences in the Lunar soil. In other words, the first 10 m of the landing site could be composed of varying rock distributions from different ejecta buried in weathered soils with similar low dielectric properties.

Conclusions
A novel set of interpretation tools for in-situ planetary radar is described. The proposed framework consists of a novel stochastic hyperbola fitting that estimates the range of permittivities that fit a given hyperbola. Stochastic hyperbola fitting considers the inherent non-uniqueness of the problem when the subsurface targets have varying sizes. Subsequently, a novel probabilistic inversion is applied that utilises multiple hyperbolas to infer the 1D permittivity profile of the area. The estimated permittivity profile is then used as input to our bespoke reverse-time migration coupled with finite-differences timedomain method, capable of migrating the processed radagram subject to any arbitrary permittivity structure. Lastly, a post-migration processing is applied, which effectively clusters the data into rocks/boulders and Lunar soil. The proposed scheme has been numerically validated and subsequently applied to the radar data from the first two Lunar days of the Chang'E-4 mission. A layered structure is detected within the first 10 m of the regolith, with a stochastic rock distribution, exhibiting areas with low rock densities, that can be potential targets for future drilling missions. This approach can be equivalently applied to radar data from Zhurong and Perseverance, which are currently exploring the surface of Mars.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request I. Giannakis et al. Fig. 15. The normalised rock abundance with respect to depth as estimated from ( , ) in Fig. 14 (E). It is apparent that there are 3 clearly seen layers within the first 10 m of the regolith. Notice that the results from the migration are in good agreement with the independent results obtained using the 1D inversion (see Fig. 13).