© The Authors 2023

1 Introduction

Radar tomography constitutes an ill-posed and ill-conditioned inverse problem in which the internal permittivity is to be found by observing the electromagnetic field transmitted and measured by the radar and penetrating through the target. That is, the solution is non-unique, and only a small amount of measurement or modelling errors can have a significant effect on the final reconstruction. Our aim is to advance the development of an inversion methodology to be used in potential planetary scientific tomography investigations (Asphaug et al. 2002, 2008) targeting small Solar System bodies (SSSBs) such as asteroids and comets, which have complex shapes and whose interior structure is largely unknown. Discovering the interior structure is an important scientific objective of Solar System research in order to understand its history and evolution. The sounding radar is, in general, one of the most promising techniques to directly yield information about the interior of SSSBs (Herique et al. 2018). The first tomographic radar measurements for an SSSB were performed as a part of the European Space Agency's (ESA's) Rosetta mission, which included the COmet Nucleus Sounding Experiment by Radiowave Transmission (CONSERT; Kofman et al. 2007, 2015; Herique et al. 2019a) as a part of its science programme. In addition to the complex scattering phenomena, specific challenges related to inverting in situ measurements including the incompleteness of the data with respect to both point and spectral coverage as well as the small carrier wavelength in comparison to the target body size. To understand how these features can be taken appropriately into account in the tomographic reconstruction process, it is necessary to perform experimental laboratory studies with analogue models.

Our concentration is, in particular, on asteroids, whose compositional structure (Britt & Lebofsky 1997) can be divided into several taxonomies based on their morphology or mineralogy. In the density studies, asteroids have been observed to contain large amounts of porosity (Carry 2012), the distribution of which is unknown, allowing the potential existence of internal macro-scale inhomogeneities such as cracks and voids (Sorsa et al. 2019). Impact simulation studies suggest that the density of the interior grows towards the centre of the body (Jutzi & Benz 2017). That is, the deeper interior structure might be covered by a less dense mantle which might be composed, for example, by powdery minerals. Some asteroids are also known to contain significant amounts of metals that might form areas impenetrable by a radar signal. Missions to discover the interior structure of asteroids include both sample return missions such as the Japanese Space Agency's (JAXA) Hayabusa-2 mission and the National Aeronautics and Space Administration's (NASA) OSIRIS-REx, as well as radar missions such as ESA's future confirmed HERA mission, whose Juventas Radar (JuRa; Herique et al. 2019b) will perform low-frequency penetrating radar measurements to reconstruct the interior structure of the Dimorphos asteroid. With the successful outcome of the DART mission (Rivkin et al. 2021) according to initial reports from October 2022 showing that both the orbital speed and the trajectory of the body were affected and a crater formed on the impacted target Dimorphos, the meeting of HERA with the same body in December 2026 and the following science experiments involving also radar measurements will be of a special interest.

Our work, structured in two joint papers (this one and Dufaure et al. 2023), concerns the investigation of an experimental tomographic reconstruction of high-contrast details inside a structurally complex target body, using experimental data. The present paper investigates the problem with a time domain method and Dufaure et al. (2023) studies the problem in the frequency domain. Our focus is on inverting quasi-monostatic back-scattering measurements incorporating advanced full-wave forward simulation (Eyraud et al. 2020; Sorsa et al. 2020) results and data obtained from multiple measurement points. In particular, multiple scattering and multipath effect suppression (MES) is considered to reduce artefacts in the reconstruction of the internal permittivity. As a potential MES technique, we examined (i) geometrical optics-based pathlength thresholding, which allows the estimation of the two-way traveltime inside the body to filter out intense reflections originating from the wall opposite to the measurement point, and (ii) a peak detection method to investigate whether a data-driven approach could be used to exclude multipath effects. In the inversion stage, we used a marginalisation technique to suppress the random effects due to modelling errors by splitting a larger point several sparse sets of measurements. Marginalisation can be motivated as a way to reduce the effect of modelling errors that might follow from the small wavelength compared to the target size (Deng et al. 2021a) or the inaccuracies of the volumetric model, leading to phase differences in the signal and, thereby, affecting the final reconstruction. As a marginalisation technique we applied independence sampling (Liu 2008) proposed in a recent 2D study (Yusuf et al. 2022).

We analysed results obtained with a laboratory analogue (Sorsa et al. 2021a) that is based on the shape of asteroid 25143 Itokawa, the target of the first Hayabusa mission, and contains a void and surface layer in its interior. As an inversion technique, we used back-propagation as proposed in Sorsa et al. (2021b). The results suggest that a reconstruction from multiple points can be successfully combined and that high-contrast details inside a complex body can be localised, while demonstrating the advantage of the approaches (i) and (ii) in MES, especially in filtering the back-wall echo.

2 Materials and methods

2.1 The tomographic target

In the measurements, we used a target analogue manufactured from a dielectric plastic filament by 3D printing and matching to the shape of the asteroid 25143 Itokawa (Gaskell et al. 2008) and an interior structure containing a mantle and a deep interior void (Fig. 1). The same target object structure was used in an earlier numerical study investigating tomographic inversion of low-frequency signals (Sorsa et al. 2019) and in an analysis of signal propagation and back-propagation in the case of a single-point measurement comparing laboratory measurements with simulations (Eyraud et al. 2020; Sorsa et al. 2021b). Our basic assumption when building the simple model was based on the results of impact simulations (Jutzi & Benz 2017) showing that sub-catastrophic impacts of elongated bodies can lead to bi-lobed structures where the porosity of the body decreases towards the centre of the body. As our model is a rubble-pile asteroid Itokawa, we assume that the overall bulk porosity (up to 40% as reported by Fujiwara et al. 2006), consists of overall material porosity in addition to cavities and/or voids inside the rubble-pile structure.

To obtain a 3D-printable structure, the edges of the conforming the finite element mesh constituting the numerical asteroid model were inflated to yield a permittivity-controlled wireframe according to the procedure detailed and analysed in Sorsa et al. (2021a). The target consists of an interior compartment (ε_r = 3.40 + j0.04), deep interior void (ε_r = 1.00) and a surface layer (ε_r = 2.56 + j0.02; Table 1) constituting the detailed model (DM) in this investigation. The electric permittivity of different compartments within the model is controlled by varying the edge width of the wireframe based on classical mixing formulae for the effective permittivity of a mixture. The wireframe was manufactured using the commercially available low-loss Preperm ABS450 filament (Premix Oy, Finland) whose relative electric permittivity at 2.4 GHz frequency is 4.50 + j0.019. The analogue is a low-loss body, which is a reasonable assumption for an asteroid (Kofman 2012).

Fig. 1 Analogue based on the shape model of the asteroid Itokawa, whose larger lobe is called the body and the smaller one the head. Top: target detailed model (DM) analogue finite element structure showing the surface layer and the deep interior void. Bottom: wireframe surface structure of the finite element mesh in which the edges have been substituted with prisms to enable 3D printing of the object.

2.2 Microwave radar laboratory measurements

The experimental measurements on the tomographic targets were carried out in the anechoic chamber of the Centre Commun de Recherches en Micro-ondes (CCRM) in Marseille, France (Geffrin et al. 2009). The experiment was carried out in an igloo-type spherical configuration of 2372 quasi-monostatic measurements around the object covering the full range of the azimuth circle, the angles from 26° to 154° in the elevation direction, and the angular steps according to the Nyquist criterion with respect to the maximum frequency (Bucci & Isernia 1997). The elevation angle between a transmitter and a receiver at each position was fixed at 12 degrees.

The antennas were placed at a 1.794 metre radius from the centre of the target. The polarisation of the antennas was linear along the elevation unit vector for both the transmitter and the receiver. The covered frequency range was between 2 and 18 GHz in 0.05 GHz increments. The measurements were made frequency-per -frequency, resembling a setup in a continuous wave radar. The largest dimension of the laboratory target (20.5 cm) and the available microwave frequency range correspond to a low-frequency radar measurement with the respective centre frequency and bandwidth of 4.92 MHz and 2.20 MHz on a real 535 m asteroid in space (Table 1).

2.3 Modelling and simulations

2.3.1 The direct problem: Full wave modelling

In addition to the laboratory measurements, a numerically simulated data set was created for the same target with a full-wave approach using the finite-element time-domain (FETD) method. Similarly to our earlier studies (Sorsa et al. 2020; Eyraud et al. 2020), the radar signal in time was modelled with an isotropic Blackman-Harris pulse with a 12.9 GHz centre frequency and a 5.70 GHz bandwidth. The necessity of the full wavefield modelling follows from its ability to distinguish scattering from different parts of a complex and high-contrast target body, also covering the higher order scattering effects of multipath signal propagation.

2.3.2 The inverse problem

The inverse problem is solved by a simple back-propagation analysis relying on the first-order Born approximation (BA) of the demodulated full wavefield. The BA is found with respect to a numerical domain which has the exterior shape of the actual target, but a homogeneous permittivity distribution matching with the constant interior value of the DM. Therefore, the numerical computations were run for two target structures, the DM shown in Fig. 1 and a homogeneous model (HM), the latter of which has the target asteroid shape with a constant interior permittivity ε_r = 3.40 + j0.04 matching the permittivity of the interior compartment of the DM. The HM interior, however, lacks the structural details (mantle, void) present in the DM.

The first-order BA (Sorsa et al. 2021b) provides a linearised mapping: $$Lx+n=y.$$(1)

There, $L\in {\text{ℝ}}^{I_I\times I_F}$ is the BA of the differentiated signal, where I_I and I_F are the numbers of elements in inversion and forward meshes, respectively. The variable $x\in {\text{ℝ}}^{I_F\times 1}$ denotes a discretised estimate of the perturbation between the actual permittivity distribution and its background value, n is a measurement noise vector, and $y\in {\text{ℝ}}^{I_I\times 1}$ is a vector containing the full-wave measurements. A back-propagated estimate for a given data can be obtained as $$x^{†}=L^{T}y.$$(2)

2.4 Marginalising phase errors by sampling and averaging

The measured data sets y₁, y₂,…, y_n, where n is the total number of independent data sets, can be interpreted to correspond to data set k specific data vectors ${\tilde{y}}_k$ according to the formula $$y_k={\tilde{y}}_k+n+{\Delta }_k,$$(3)

where n is a noise measurement term and Δ_k denotes a modelling error, which, based on the earlier investigation (Eyraud et al. 2020), can be assumed greater than the inaccuracy of the actual measurement. Hence, |n| ≪ |Δ_k| for all k = 1,2,…, n. The resulting signal oscillation artefacts induced by the modelling error and depicted schematically in Fig. 2 are difficult to predict due to the relatively small wavelength of measurement with respect to the target size and the non-linearity of the wave equation with respect to the unknown permittivity. This is especially so because they can occur even if the differences between the actual and estimated permittivity distributions in the numerical model were small. Due to such modelling errors, predicting the mutual correlation between two data sets measured at different points is likely to be difficult without accurate a priori knowledge of the global structure of the permittivity estimate x.

Because coarse distributions are generally less sensitive to this noise and can be found with fewer data points and lower frequency data than fine details in linear or linearised inverse problems (Pursiainen 2008), the coarse global permittivity structure x can be approximated based on a sparse distribution of measurement points. As a measure of sparsity, we considered the minimum angular aperture, α > 0, between two different measurement points in the set. All possible data vector realisations with sparsity are denoted by 𝒟_α. To denote the kth estimator of x, we used $$x_k^{*}=F\left(y_k\right).$$(4)

The function F: 𝒟_α → 𝒮 from the set 𝒟_α to the set 𝒮 of candidate solutions is assumed to be determined by an inversion technique which by itself incorporates global a priori information into the process of recovering x, which is crucial due to the ill-posedness of the inverse problem. Interpreting the data vectors as mutually independent (or weakly correlated) and identically distributed (i.i.d.) random variables for k = 1,2,…,n, the sample mean $\frac{1}{n}{\displaystyle {\sum }_{k=1}^n{x}_k\to x^{*}}$ approximates a normally distributed random variable with mean x^* and a standard deviation proportional to $1/\sqrt{n}$ (Chung & Zhong 2001; Liu 2008).

Furthermore, $$x^{*}\approx {\displaystyle \sum _{k=1}^{n}F\left(x_k^{*}\right)\approx {\displaystyle {\int }_{{𝒟}_{\alpha }}F}\left(y\right)\text{d}y.}$$(5)

In the present study, incorporating full-wave data from multiple measurement points poses a challenge due to the high centre frequency of the signal wave. This means that the phases of the measured and simulated data might not match, since the phase accuracy of the simulation decreases when the frequency increases, inducing mutual position-based oscillating differences between the measured and simulated signals as demonstrated in Fig. 2. These oscillations can vary between measurement points and are likely to intensify when time increases, as longer propagation times correspond to higher order scattering effects and multipath signal propagation. A similar approach has been shown to improve reconstruction quality in 2D (Yusuf et al. 2022).

2.5 Sparse measurement configuration point sets

The large laboratory measurement point set was down-sampled to obtain a computationally convenient number of points for the purpose of present full-wave modelling and inversion analysis. We investigated two sparse measurement configurations, one containing 92 points and the other 363 points out of the 2372 points measured in the laboratory. The two configurations were further divided into 10 and 40 sub-configurations consisting very sparse measurement point clouds (Fig. 3). Table 2 describes the four configurations that are labelled based on the number of points in each sub-configuration and the number of sub-configurations in the total data set, resulting in the data sets 92/1, 363/1, 10/10, and 40/10. In each configuration and sub-configuration, the neighbouring points are separated by a minimal angle given in Table 2 to ensure an even spread of points around the target. The resulting measurement point configurations indicating the transmitter and receiver locations and orientations are illustrated in Fig. 3 (left and middle), and an example of one of the sparse sub-configurations is given on the right.

2.6 Signal filtering based on pathlength and geometrical optics

Strong multipath effects and scattering effects, which can both affect the accuracy of the inversion scheme proposed here, were observed for the longer propagation paths in Eyraud et al. (2020) and Sorsa et al. (2021b). To filter out the effect from the inverse estimates, we estimate the travel time of the first echo from the wall opposite to the measurement point, which approximately marks the start of the higher order scattering interval in the temporal domain (Sorsa et al. 2021b). To accomplish this, we used a simple approach to estimate the signal path length based on geometrical optics, first finding the closest point between the point of transmission and the exterior surface of the target, and then the propagation path inside the target from the front wall to the back wall based on the Snell–Descartes law of the first refraction. The refractive index of the target $n=\sqrt{{\epsilon }_{\text{r}}^{\prime }}$ is calculated based on the real part ${\epsilon }_{\text{r}}^{\prime }$ of the permittivity value of the HM. The velocity of the wave inside the target is obtained as ${\text{C}}_0/\sqrt{{\epsilon }_{\text{r}}^{\prime }}$, where c₀ is the wave velocity in free space. The time interval of the measurement data utilised in the inversion process is denoted as $$T_{\text{inv}}=\gamma T_{\text{path}},$$(6)

where T_path is the two-way travel-time estimate obtained with the geometric optics-based pathlength estimate and γ > 0 is a threshold parameter giving the fraction of the total path length to back wall that is included in the signal used in the inversion stage.

Fig. 2 Schematic illustrations on the choice of measurement configuration and phase difference. (a): a single transmitter (blue) and receiver (red) pair orbiting a non-convex asteroid body. (b): two signal paths (solid and dashed), whose difference might be significant with respect to the carrier wavelength and thereby its phase. (c): outline of a demodulated measured signal y (red) and numerically simulated (dashed blue) curve $\tilde{y}$ that is in phase with respect to the measurement. (d): simulated signal partly in phase (left) and partly out of phase (right), while otherwise being similar to the signal. Despite their difference with respect to the phase, the signals in (c) and (d) can correspond to nearby measurement points. Consequently, sparse measurement configurations might be advantageous in reconstructing coarse details (oval in (a)), as incorporating data from multiple measurement points to the reconstruction might introduce artefacts in a reconstructed distribution (Yusuf et al. 2022) when found using the Born approximation of the simulated full-wave signal (striped detail in (b)). The density of the transmitter and measurement points in each configuration is depicted by the sectors of the encircling circle.

Fig. 3 Signal configurations used in the study. (a): Transmitter-receiver pairs of the 92 measurement points which are used in the 92/1 and 10/10 data sets. (b): the 363 transmitter-receiver pairs forming the 363/1 and 10/40 sets. (c): example of a sparse 10-point-pair sub-configuration. The different colours in (a) and (b) indicate the sub-systems into which each measurement point belongs to. The target has been magnified for illustration purposes.

2.7 Investigated reconstruction cases I-V

The permittivity reconstructions are evaluated for five different cases I–V (Table 3), which are different with respect to the timedomain filtering of the signal to suppress the multipath effects and back-wall scattering. The first of these, (I), uses full data including all multipath effects carried by the full wave. In cases II–IV, the full-wave signal is thresholded by Eq. (6) to investigate the effect multipath effect suppression has on the reconstruction. The threshold parameter γ is 1.1, 0.95, and 0.75 for (II), (III), and (IV), respectively. That is, the maximal signal propagation time is 110, 95, and 75% of the two-way travel time between the transmitting antenna and the back wall. Case V thresholds the signal at 80% of the maximum peak, which is assumed to be the data-driven way of determining the back-wall echo.

3 Results

3.1 Propagation inside the target

Wave propagation inside the target is visualised in Fig. 4. The wavefront transmitted from a point at the body end of the target object will enter the target at two different interfaces: first at the body end of the target, and secondly at the head end of the target due to the concavity in the neck region. The wave that travels inside the target from the body region is significantly slower than that propagating in the free space surrounding the target. When the wave enters the target in the head region, it is further away from the receiver in comparison to the part of the wave which enters the target first, causing ambiguity in identifying the origin of the time-domain signal at the receiver. Furthermore, the two wavefronts inside the target eventually combine due to scattering and further complicate the simple interpretation of the signal curves.

3.2 Time domain filtering based on path length and geometrical optics

The path length to the back wall of the target was identified using two methods: (1) geometrical optics with the threshold parameter γ = 1 and (2) data-driven identification of the first occurrence of 80% of the maximum peak in the data. Figure 5 shows the topographical maps of the path lengths to the back wall based on these two methods and traces two of the shortest and the longest paths found in the system. The topographies are shown for both the 92- and 363-point systems. The general distributions of the path lengths do not change between the two system sizes, but, as can be expected with a higher number of measurement points, the path-length distribution of the larger point system is smoother. In the geometrical optics-based path-length determination, the longest path lengths are found in the head and body regions close to the xy plane, where z is close to zero. In these areas, the wave enters the target fairly perpendicular to the front surface and is hence refracted so that it travels inside the body for a long time. The shortest times are in the middle of the body measured from the top and bottom of the z axis. This is expectantly due to the shape of the target, as the y and z dimensions of the body are approximately half that of the x direction.

The data-driven way of identifying 80% of maximum signal peaks of the HM data produces the longest paths in the head and the body regions of the target and the path lengths are in general shorter than those with the geometrical optics-based way of determining the path length. The back-wall signal is always the strongest in areas where the signal is reflected from a concave surface producing a strong lens effect that focusses the wave and leads to a more intensive signal at the receiver. The head and body regions of the target have high concave curvatures from which the signal can be reflected. These are observed as higher peaks in the data in comparison to the back-wall echoes originating from the smoother, less curved parts of the body of the target.

To investigate the multiple scattering and multipath effect suppression on the reconstruction outcome, the full-wave data of the difference of the measured DM and simulated HM for each measurement point was then filtered according to the five cases outlined in Table 3. The normalised difference curves of both the in-phase and quadrature components are shown for the complete signal (Case I), and the filtered signal data are given for the three different thresholds γ (Cases II–IV), and the 80% maximum peak detection (Case V) are shown in Fig. 6. The peaks in the data show the relative size of the difference between the measured DM and the simulated HM signals, suggesting that where the difference is high, there are internal details producing a difference in the received wave. The difference between the distributions of Cases II–IV and Case V adds to the evidence already seen in Fig. 5 that the two different filtering techniques produce different data leading to an expectation that the reconstructions are different. In general, the 80% maximum peak identification technique more often produces longer and shorter paths, although the maximum peaks do coincide roughly in all histograms. Moreover, there are clear differences in the difference data curves of the two different filtering methods. The difference curves tend to peak earlier with the geometrical optics-based thresholding than with the 80% maximum peak detection technique, suggesting that the maximum peak detection may not detect the back wall of the signal well enough for areas where the curvature of the back wall is not sufficient to produce a clear lens effect.

Fig. 4 Visualisation of wave propagation inside the target at two time points showing how the complex shape of the target affects the propagation paths. The wave, originating from outside of the lower left corner of the images, enters the target at two interfaces: first at the left body of the target and later at the head region in which the second wavefront is ahead of the first, which is travelling more slowly inside the target (left image). The scattered wavefronts can also eventually mix (right) due to a reflection back from the interior walls, producing constructive and destructive interference in addition to multipath effects of direct scattering.

Fig. 5 Comparison of the two-way travel times of path lengths from the transmitter to the back wall, assuming the homogeneous interior of the HM based on the geometrical optics (top) and data-driven maximum peak identification (bottom) techniques. The topographies are shown for the 92-point system (left), and the 363-point system (middle). The two longest (red) and shortest (blue) paths for both the filtering methods are shown on the right.

3.3 Signal-to-noise ratio and propagation paths

Earlier research on simulated data suggests that a signal-to-noise ratio (S/N) below 10 dB has a deteriorating effect on the reconstruction (Sorsa et al. 2019). Therefore, we investigated the S/N of the difference between the measured and simulated DM signals considering a constant time interval of 1.8 ns beginning from the first arrival of the wave at 11.2 ns. Figure 7 shows the topographical maps of the S/N for the obtained data over this interval and the topographical presentation of the travel time to the peak S/N for the same data and both signal configurations 92/1 and 363/1. It also includes the temporal signals used for the S/N determination for the points corresponding to the highest and lowest S/N measurement points in the data set, showing a comparison of how the measured and simulated temporal data curves differ in these points.

The topographies show that the agreement between the measurement and simulation is appropriate throughout the model. The local high and low points are found in the regions where the model surface has peculiar curvature, that is, in the neck region in the smaller system case, and in the body region where there is a sharp edge. The highest and lowest S/Ns were observed close to the regions where the topography obtains its maximum and minimum, respectively. Similar to the topographic maps related to the path lengths (Fig. 5), the S/N plots show that the main differences between signal configurations 92/1 and 363/1 relate mostly to a smoother distribution of values in the case of the larger system, and the maximum and minimum areas are found in similar locations providing evidence that the data in both signal configurations are consistent.

Fig. 6 Normalised difference curves of the measured DM and simulated HM data of the 92 measurement points showing the in-phase and quadrature components of the QAM-modulated full-wave data (Case I, top row), the effect of the threshold parameter γ (Cases II–IV, middle rows), and the 80% maximum peak detection filtering (Case V, bottom row). The histograms of the filtered cases show how the distribution of the data changes between the two filtering techniques.

3.4 Reconstructions

In addition to the visual inspection of the reconstructions, the quality of a reconstruction was evaluated quantitatively by volumetric overlap of the reconstructed interior to the exact locations of the details in the analogue. The centre-of-mass difference and the root-mean-squared error (RMSE) and the correlation between the reconstruction and the exact structure were calculated. The results are shown in Table 4.

The overlap was evaluated by two relative overlap values for (1) the exact location of the void in the DM versus reconstruction and (2) the union of the exact location of the void and mantle in the DM versus reconstruction. In cases I and II, the overlap was calculated in relation to the equally large sets, referred to as reconstruction layer 1 and 2, respectively, where the absolute value of the reconstruction is greater than in its complement.

The third measure of the goodness of the reconstruction is the distance between the centre of mass (MC) of the reconstruction set and that of the exact DM. The MC shows that the averaged sparsity-based reconstructions were more regular when signal thresholding based on path length was applied. The RMSE values reflect the large error in localising the void and the surface layer (where the overlap between the reconstructed and the exact details are small, the quantitative error is large).

Figure 8 shows the best two reconstructions obtained by configurations 92/1 and 10/10 as evaluated by these quantitative measures, and it simultaneously shows the effect of marginalising the phase errors by sampling the measurement point set and averaging the reconstruction outcomes to yield the final result. Configuration 92/1, Case IV, yielded the maximum overlap for reconstruction layers 1 and 2 (47.1% and 63.5%, respectively), and the second one, 10/10, Case IV, gave the corresponding values of 6.56% and 61.3%. These values were obtained with geometric optics-based travel times of up to 75% of the two-way trip between the back wall and antenna, and this approach minimised the artefacts in reconstruction layer 1 (void), that is, parts not corresponding to the void. The regularising effect of the averaging process can be clearly observed in the shape of the reconstruction layers.

The effect of selecting the travel-time threshold can be observed in Figs. 9 and 10 (showing reconstruction layers 1 and 2) for each signal configuration and a signal thresholding case as 2D projection and equisurfaces in 3D, respectively. The longer path lengths tend to result in artefacts, especially with an inclusion in the smaller of the lobes, which suggests that the artefact may be caused by some other process than a multiple scattering effect.

Decreasing the value of the threshold γ results in the diminishing of this artefact with configurations 92/1 and 10/10, while it also affects the shape of the inclusion(s) corresponding to the void. This is obvious with the averaged reconstructions, for which the overlap between the void and reconstruction layer 1 decays as γ decreases. The threshold has a smaller effect on the artefacts in the case of the denser point configurations 363/1 and 10/40, which is due to more points where the geometric optics-based travel-time estimate for the path length to the back wall differs at least partly from the realised scattering.

While the reconstructions of the 10/10 Case III appears as the best based on visual inspection, especially regarding the shape localisation of the void, its quantitative measures on the overlap volume 2 and MC difference are inferior to Case IV. Also, the artefact in the smaller lobe of the target decreases the quantitative value for goodness.

For each measurement configuration, the best overlap and MC were obtained when the travel time was thresholded to exclude the multipath effects and the back-wall echo. The preferable threshold value γ was in the range between 0.75 and 0.95, suggesting that it is advantageous to filter the temporal signal well before the expected back-wall reflection.

Although the overlap values given in Table 4 suggest that, with full data, the reconstructions obtained with configurations 92/1 and 363/1 were superior to those obtained via averaging, that is, configurations 10/10 and 10/40, a visual inspection of Figs. 9 and 10 reveals that the smoother reconstructions of the averaged cases have a deteriorating effect on the overlap values. Notably, the averaged reconstructions are able to locate the void in nearly all cases (excluding case V). However, the exact localisation of the void is not optimal, which is reflected well in the quantitative measures.

The overall best reconstruction by both quantitative and qualitative inspection is configuration 92/1 in Case IV. The reconstruction shows the location of the void well and the overlap value (47.1%) is distinctively superior to any other reconstruction. Unfortunately, similar performance and quality of reconstruction cannot be reached when the number of measurements increases, as shown by the same case IV of configuration 363/1, where the void overlap is 0.0% and the overlap of reconstruction layer 2 is also smaller in comparison to the other methods.

Overall, the correspondence between reconstruction layer 2 and the surface layer seems weaker than that between layer 1 and the void. Figures 9 and 10 show that reconstruction layer 2, which optimally covers both the void and the surface layer, is only partly distributed over the surface layer. The non-averaged reconstructions seem to find parts of the surface layer, especially in Case IV of configuration 92/1, which corresponds to the maximum overlap. The averaged reconstructions do not find the surface layer despite their overlap; that is, reconstruction layers 1 and 2 are connected without a gap between them.

Fig. 7 Topographical presentation of the peak S/N and the travel time to peak S/N for the 92/1 (top) and 363/1 (bottom) configurations along with the signal curves for the highest (red plot) and lowest (blue plot) S/N points in the measurement system.

Fig. 8 Demonstration of the effect of averaging over sub-configurations on the reconstruction yielding the best overlap values for the mantle and the void (Table 4). Left: reconstruction obtained in Case IV of signal configuration 92/1. The asteroid mantle is depicted by the green surfaces and the void by the blue ellipsoid. Right: enhanced regularity of reconstruction layer 1 was obtained with the averaged reconstruction corresponding to configuration 10/10 in Case IV, although the overlap between the void and the reconstruction is inferior to that of configuration 92/1.

3.5 Computational cost

The wave simulation was performed in the Puhti cluster of CSC, where each wavefield was modelled on a GPU (Nvidia Volta V100 GPU with 32 GB of HBM2 memory). The inversion computations were performed with a workstation equipped with a 20-core processor (Intel i9-10900X, 3.70GHz) and 128 GB RAM. The inversion process consisted of the BA and back-propagation stages, of which the former required the highest computation cost and longest processing time, while the latter could be performed as given by Eq. (2) in less than one second for a given data vector.

4 Discussion

This study concerned reconstructing the internal structure of a permittivity-controlled plastic asteroid analogue model based on a large spatio-temporal set of experimental quasi-monostatic microwave measurements. The focus was on combining a full-wavefield modelling approach with the multiple scattering and MES technique and marginalisation of the numerical phase errors. Indeed, MES can be considered an inevitable part of the modelling process due to the complex geometry and high contrast of the permittivity distribution, which causes significant signal deviations between different measurement points and intense scattering peaks, in particular the echo from the concave back wall opposite to the wall facing the measurement point. MES was approached via a geometric optics-based thresholding path length to the back wall and determining the back wall based on the maximum peak detection. Then, in the inversion stage, averaging the results for multiple different subsets of measurement data was used to marginalise modelling errors to improve the reconstructions and enable the use of a sparse data set.

The results show the importance of both path-length-based thresholding and marginalisation of modelling errors by averaging, which were found to be necessary to (1) reduce the artefacts resulting from modelling errors and (2) regularise the outcome of the reconstruction process. Obtaining the best possible estimate for the void requires estimating the travel time of the signal and averaging. A geometric optics-based estimate for the travel time was found to be preferable compared to the maximumamplitude thresholding of the full wave. This is most probably due to the difficulty of determining the back-wall echo in a data-driven way in cases where the curvature of the backwall is small and varied, and hence the lens effect, which produces the most notable peaks in the received amplitude, cannot be accurately found. The preferable threshold of the travel time obtained via thresholding the path length to the back wall was found to be from 75 to 95% to exclude the effect of back-wall echo on the reconstruction. Averaging i.e, marginalising the reconstruction over several subsets of data as suggested in Yusuf et al. (2022) resulted in the small centre of mass difference between reconstruction layer 1 and the void, while the non-averaged approach yielded the maximal overlap for both reconstruction layers 1 and 2.

For the dense and sparse measurement point clouds 363/1 and 92/1, the greater effect of thresholding was observed in the case latter, which, due to its lower point density, includes fewer points in which the geometrical optics-based estimate differs significantly from the actual wave propagation. For the present shape model of asteroid 25143 Itokawa, these differences may be assumed to occur, especially, at the measurement points facing the smaller of the two lobes, where on one hand the length of the paths is lower than average, and, on the other hand, some paths extend across the full asteroid body, and, moreover, the angular difference between the long and short paths is small due to the high curvature of the surface. Thus, some paths, when treated as rays, match for longer propagation times, although they might at least partly correspond to significantly shorter propagation due to the non-linearity of the wave equation. The intense anomaly observed in the smaller lobe can be recognised as a result of this effect. For the sparser data set, this anomaly could be significantly reduced by controlling the threshold of the travel time, while for the denser set it was not possible.

Compared to the thresholding based on geometrical optics, the maximum peak thresholding resulted in shorter travel times in general, and especially in the direction of the elongated axis of the Itokawa model. This data-driven method could best perform with symmetrical and near-spherical targets where the curvature of the back wall is fairly constant throughout the object, and hence a maximum peak in the data would most likely be caused by the back wall. This is also why the selected maximum peak threshold did not detect the back-wall echo uniformly in this target model; this can be attributed to the complexity of the geometry due to which the effect of the back wall varies depending on the measurement location.

As shown by our analysis, part of the wavefronts propagating inside the asteroid can be divided into two branches occurring in the two lobes simultaneously. In addition to the path length, branching is a potential source of artefacts, since an obstacle (void) affecting one of the branches can be reflected in another one. Branching was observed to correlate with the direction of the elongated diameter, as the wavefront inside the body propagates more slowly than the one outside, allowing the latter one to reach the second lobe before the wave exits the first one. The artefact observed in the smaller lobe may also partly originate from branching. Nevertheless, since only some of the paths are branched, we deem branching as a nuisance effect to be a minor one compared to the back-wall echo, which can mask smaller amplitude echoes from the interior of the target. Branching and back-wall scattering together explain the lower S/N that was observed to roughly match the elongated direction based on the topographical maps.

The analysis presented in this study is relevant regarding the design of planetary missions targeting small bodies such as the HERA mission, which will carry the Juventas radar (Herique et al. 2019b, 2020). Due to the complex effects observed with the model, advanced wave propagation models will be necessary in such a context. Based on our results, full-wavefield modelling is an advantageous method to optimise the S/N of the wave propagation, which, in the context of planetary missions, has also been proposed and analysed recently in Deng et al. (2021b); Sava & Asphaug (2018a,b); while we present the first results for a highly controlled multipoint asteroid analogue using a time-domain approach.

Radar measurements made by a spacecraft will obviously involve more sources of uncertainty, for example, positioning and scaling errors, than the laboratory measurements processed in this study. While those are not covered in study, our observation with the current experiment setting was that the effect of scaling on the reconstructions is not significant if the scaling error is less than 3%.

Various methods and combinations might be applicable to process the simulated and measured data to filter out the nuisance of the full wavefield. The combination of geometric optics-based temporal data filtering and full-wavefield modelling applied here has been investigated in various studies such as Mocker et al. (2015). For the above justification, an advanced ray-tracing model (Gassot et al. 2020) applied for a comet in Herique et al. (2019a) might perform sufficiently in the inversion stage. A ray-tracing model to be used in connection with rough surfaces has been proposed in Cocheril & Vauzelle (2007). Both full-wavefield and ray-tracing approaches are generally used in advanced tomographic imaging. A recent comparison for an ultrasonic tomography setting can be found in Bancel et al. (2021).

Fig. 9 Cut-view images for the yz, xz, and xy planes (left, middle, and right, respectively) of the target. Reconstruction layers 1 (red) and 2 (grey) are visualised for each reconstruction case ((I)–(V)) and all four measurement configurations. The exact structures of the mantle and the void are depicted by the green and black areas, respectively. Optimally, layers 1 and 2 exactly match these areas.

Fig. 10 Three-dimensional visualisation of reconstruction layers 1 (red) and 2 (grey) as equisurfaces. The exact location of the asteroid mantle and interior are depicted by the two green surfaces and the void by the blue ellipsoid.

5 Conclusion

Imaging of the internal structure of an asteroid requires advanced data analysis methodologies and high-performance computing to account for the complexities of surface and internal scattering in a realistic target object. The significant back-wall echoes caused by the lens effect originating from scattering from concave back wall of the target body can mask the other echoes from internal structures. Therefore, filtering that signal from the data by the presented multiple scattering and multipath effect suppression (MES) method, combined with the simple back-propagation inversion of the data, can yield meaningful information on the interior structure of a complex tomographic target.

This paper, together with its companion paper (Dufaure et al. 2023), investigates methods used to image the internal structure of an asteroid analogue. This second part of the study is concentrated on a special case where the reconstruction is sought from time-domain data using a Born approximation and a linearised inversion model with a sparse data set. As shown in the results, a meaningful reconstruction can be found, and the deep interior void is detected in an asteroid analogue. The localisation of the void, however, is not perfect, as reflected by the relative overlap values and the error quantities, in addition to the qualitative reconstructions. The finer detail of the surface layer could not reliably be reconstructed with the studied method. In this study, we also observed geometry-related artefacts caused by the bi-lobed and elongated structure of the analogue. Further studies are needed to develop a methodology to differentiate such artefacts and actual interior details.

A future work involves performing a similar analysis and comparative analysis with alternative full-wave techniques for this multipoint tomographic measurement and inversion. The first such study in the frequency domain is presented in the joint paper (Dufaure et al. 2023).

Acknowledgement

The authors acknowledge the opportunity provided by the Centre Commun de Ressources en Micro-ondes (CCRM) to use its fully equipped anechoic chamber. Centre de Calcul Intensif d'Aix-Marseille is acknowledged for granting access to its high performance computing resource. L.-I.S. was supported by a Young Researcher's Grant from Emil Aaltonen Foundation. Y.O.Y. was supported by the Magnus Ehrnrooth Foundation through the graduate student scholarship. L.-I.S., Y.O.Y., and S.P. were also supported by the Centre of Excellence in Inverse Modelling and Imaging (Academy of Finland 2018-2025, project number 312341) and the Academy of Finland, project number 336151. L.-I.S. and S.P. acknowledge CSC – IT Center for Science Ltd., for providing computing services on the Puhti supercomputer.

References

All Tables

All Figures

Fig. 1 Analogue based on the shape model of the asteroid Itokawa, whose larger lobe is called the body and the smaller one the head. Top: target detailed model (DM) analogue finite element structure showing the surface layer and the deep interior void. Bottom: wireframe surface structure of the finite element mesh in which the edges have been substituted with prisms to enable 3D printing of the object.

In the text

Fig. 2 Schematic illustrations on the choice of measurement configuration and phase difference. (a): a single transmitter (blue) and receiver (red) pair orbiting a non-convex asteroid body. (b): two signal paths (solid and dashed), whose difference might be significant with respect to the carrier wavelength and thereby its phase. (c): outline of a demodulated measured signal y (red) and numerically simulated (dashed blue) curve $\tilde{y}$ that is in phase with respect to the measurement. (d): simulated signal partly in phase (left) and partly out of phase (right), while otherwise being similar to the signal. Despite their difference with respect to the phase, the signals in (c) and (d) can correspond to nearby measurement points. Consequently, sparse measurement configurations might be advantageous in reconstructing coarse details (oval in (a)), as incorporating data from multiple measurement points to the reconstruction might introduce artefacts in a reconstructed distribution (Yusuf et al. 2022) when found using the Born approximation of the simulated full-wave signal (striped detail in (b)). The density of the transmitter and measurement points in each configuration is depicted by the sectors of the encircling circle.

In the text

Fig. 3 Signal configurations used in the study. (a): Transmitter-receiver pairs of the 92 measurement points which are used in the 92/1 and 10/10 data sets. (b): the 363 transmitter-receiver pairs forming the 363/1 and 10/40 sets. (c): example of a sparse 10-point-pair sub-configuration. The different colours in (a) and (b) indicate the sub-systems into which each measurement point belongs to. The target has been magnified for illustration purposes.

In the text

Fig. 4 Visualisation of wave propagation inside the target at two time points showing how the complex shape of the target affects the propagation paths. The wave, originating from outside of the lower left corner of the images, enters the target at two interfaces: first at the left body of the target and later at the head region in which the second wavefront is ahead of the first, which is travelling more slowly inside the target (left image). The scattered wavefronts can also eventually mix (right) due to a reflection back from the interior walls, producing constructive and destructive interference in addition to multipath effects of direct scattering.

In the text

Fig. 5 Comparison of the two-way travel times of path lengths from the transmitter to the back wall, assuming the homogeneous interior of the HM based on the geometrical optics (top) and data-driven maximum peak identification (bottom) techniques. The topographies are shown for the 92-point system (left), and the 363-point system (middle). The two longest (red) and shortest (blue) paths for both the filtering methods are shown on the right.

In the text

Fig. 6 Normalised difference curves of the measured DM and simulated HM data of the 92 measurement points showing the in-phase and quadrature components of the QAM-modulated full-wave data (Case I, top row), the effect of the threshold parameter γ (Cases II–IV, middle rows), and the 80% maximum peak detection filtering (Case V, bottom row). The histograms of the filtered cases show how the distribution of the data changes between the two filtering techniques.

In the text

	Fig. 7 Topographical presentation of the peak S/N and the travel time to peak S/N for the 92/1 (top) and 363/1 (bottom) configurations along with the signal curves for the highest (red plot) and lowest (blue plot) S/N points in the measurement system.
In the text

Fig. 8 Demonstration of the effect of averaging over sub-configurations on the reconstruction yielding the best overlap values for the mantle and the void (Table 4). Left: reconstruction obtained in Case IV of signal configuration 92/1. The asteroid mantle is depicted by the green surfaces and the void by the blue ellipsoid. Right: enhanced regularity of reconstruction layer 1 was obtained with the averaged reconstruction corresponding to configuration 10/10 in Case IV, although the overlap between the void and the reconstruction is inferior to that of configuration 92/1.

In the text

Fig. 9 Cut-view images for the yz, xz, and xy planes (left, middle, and right, respectively) of the target. Reconstruction layers 1 (red) and 2 (grey) are visualised for each reconstruction case ((I)–(V)) and all four measurement configurations. The exact structures of the mantle and the void are depicted by the green and black areas, respectively. Optimally, layers 1 and 2 exactly match these areas.

In the text

	Fig. 10 Three-dimensional visualisation of reconstruction layers 1 (red) and 2 (grey) as equisurfaces. The exact location of the asteroid mantle and interior are depicted by the two green surfaces and the void by the blue ellipsoid.
In the text