Issue 
A&A
Volume 674, June 2023



Article Number  A73  
Number of page(s)  12  
Section  Planets and planetary systems  
DOI  https://doi.org/10.1051/00046361/202244778  
Published online  02 June 2023 
Imaging of the internal structure of an asteroid analogue from quasimonostatic microwave measurement data
II. The time domain approach
^{1}
Computing Sciences, Tampere University (TAU),
PO Box 692,
33101
Tampere, Finland
email: liisaida.sorsa@tuni.fi
^{2}
AixMarseille Univ, CNRS, Centrale Marseille, Institut Fresnel,
Marseille, France
Received:
19
August
2022
Accepted:
13
January
2023
Context. The internal structures of small solar system bodies (SSSBs) are still poorly understood. In this paper, we find an experimental tomographic reconstruction of coarse highcontrast details inside a complexstructured target object using multipoint fullwave radar data.
Aims. Our aim is to advance the development of inversion techniques to be used in potential planetary scientific radar investigations targeting SSSBs, which have complex shapes and whose internal structure is largely unknown. Finding out the structure is an important scientific objective of Solar System research in order to understand its history and evolution.
Methods. This is the second part (Paper II) of a joint study considering the methods to analyse and invert quasimonostatic microwave measurement data of an asteroid analogue. We focused on incorporating advanced, fullwave, forward simulation in time domain with experimental data obtained from multiple measurement points. In particular, this study investigates multiple scattering and multipath effect suppression (MES) to reduce artefacts in the reconstructions. MES is necessary since the highcontrast and complexshaped target and, especially, its back wall in high curvature regions cause intense reflections that deteriorate the reconstruction quality if not treated correctly. We considered the following two approaches to obtain MES: (i) geometrical opticsbased pathlength thresholding and (ii) a peak detection method to investigate whether a datadriven approach could be used. At the inversion stage, we investigated marginalisation of random effects due to modelling by splitting a larger point set into several sparse sets of measurements.
Results. Based on the results, MES is crucial to localise a void inside the complex analogue target. A reconstruction can be found when the maximum signal propagation time approximately matches that of the first backwall echo for each measurement point. The marginalisation approach allows us to find a reconstruction that is comparable in quality to the case of full data, while reducing the computation effort per subsystem, which is advantageous when inverting a large data set.
Key words: scattering / methods: data analysis / minor planets / asteroids: general / planets and satellites: interiors / techniques: image processing / waves
© The Authors 2023
Open Access article, published by EDP Sciences, under the terms of the Creative Commons Attribution License (https://creativecommons.org/licenses/by/4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
This article is published in open access under the Subscribe to Open model. Subscribe to A&A to support open access publication.
1 Introduction
Radar tomography constitutes an illposed and illconditioned 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 nonunique, 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 macroscale 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) Hayabusa2 mission and the National Aeronautics and Space Administration's (NASA) OSIRISREx, as well as radar missions such as ESA's future confirmed HERA mission, whose Juventas Radar (JuRa; Herique et al. 2019b) will perform lowfrequency 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 highcontrast 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 quasimonostatic backscattering measurements incorporating advanced fullwave 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 opticsbased pathlength thresholding, which allows the estimation of the twoway 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 datadriven 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 backpropagation as proposed in Sorsa et al. (2021b). The results suggest that a reconstruction from multiple points can be successfully combined and that highcontrast details inside a complex body can be localised, while demonstrating the advantage of the approaches (i) and (ii) in MES, especially in filtering the backwall 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 lowfrequency signals (Sorsa et al. 2019) and in an analysis of signal propagation and backpropagation in the case of a singlepoint 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 subcatastrophic impacts of elongated bodies can lead to bilobed structures where the porosity of the body decreases towards the centre of the body. As our model is a rubblepile 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 rubblepile structure.
To obtain a 3Dprintable structure, the edges of the conforming the finite element mesh constituting the numerical asteroid model were inflated to yield a permittivitycontrolled 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 lowloss Preperm ABS450 filament (Premix Oy, Finland) whose relative electric permittivity at 2.4 GHz frequency is 4.50 + j0.019. The analogue is a lowloss 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 Microondes (CCRM) in Marseille, France (Geffrin et al. 2009). The experiment was carried out in an iglootype spherical configuration of 2372 quasimonostatic 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 frequencyper 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 lowfrequency 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).
Realsize asteroid measurement parameters corresponding to the laboratory parameters.
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 fullwave approach using the finiteelement timedomain (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 BlackmanHarris 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 highcontrast 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 backpropagation analysis relying on the firstorder 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 firstorder BA (Sorsa et al. 2021b) provides a linearised mapping: (1)
There, 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 denotes a discretised estimate of the perturbation between the actual permittivity distribution and its background value, n is a measurement noise vector, and is a vector containing the fullwave measurements. A backpropagated estimate for a given data can be obtained as (2)
2.4 Marginalising phase errors by sampling and averaging
The measured data sets y_{1}, y_{2},…, y_{n}, where n is the total number of independent data sets, can be interpreted to correspond to data set k specific data vectors according to the formula (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 nonlinearity 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 (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 illposedness 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 approximates a normally distributed random variable with mean x^{*} and a standard deviation proportional to (Chung & Zhong 2001; Liu 2008).
In the present study, incorporating fullwave 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 positionbased 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).
Spatial configurations of transmissionmeasurement point pairs used in this study.
2.5 Sparse measurement configuration point sets
The large laboratory measurement point set was downsampled to obtain a computationally convenient number of points for the purpose of present fullwave 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 subconfigurations 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 subconfiguration and the number of subconfigurations in the total data set, resulting in the data sets 92/1, 363/1, 10/10, and 40/10. In each configuration and subconfiguration, 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 subconfigurations 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 is calculated based on the real part of the permittivity value of the HM. The velocity of the wave inside the target is obtained as , where c_{0} is the wave velocity in free space. The time interval of the measurement data utilised in the inversion process is denoted as (6)
where T_{path} is the twoway traveltime estimate obtained with the geometric opticsbased 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 nonconvex 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 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 fullwave 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): Transmitterreceiver pairs of the 92 measurement points which are used in the 92/1 and 10/10 data sets. (b): the 363 transmitterreceiver pairs forming the 363/1 and 10/40 sets. (c): example of a sparse 10pointpair subconfiguration. The different colours in (a) and (b) indicate the subsystems into which each measurement point belongs to. The target has been magnified for illustration purposes. 
Five different timedomain filtering techniques applied to the fullwave signal.
2.7 Investigated reconstruction cases IV
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 backwall scattering. The first of these, (I), uses full data including all multipath effects carried by the full wave. In cases II–IV, the fullwave 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 twoway 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 datadriven way of determining the backwall 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 timedomain 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) datadriven 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 363point 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 pathlength distribution of the larger point system is smoother. In the geometrical opticsbased pathlength 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 datadriven 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 opticsbased way of determining the path length. The backwall 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 backwall 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 fullwave 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 inphase 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 opticsbased 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 twoway 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 datadriven maximum peak identification (bottom) techniques. The topographies are shown for the 92point system (left), and the 363point system (middle). The two longest (red) and shortest (blue) paths for both the filtering methods are shown on the right. 
3.3 Signaltonoise ratio and propagation paths
Earlier research on simulated data suggests that a signaltonoise 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 inphase and quadrature components of the QAMmodulated fullwave 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 centreofmass difference and the rootmeansquared 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 sparsitybased 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 opticsbased travel times of up to 75% of the twoway 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 traveltime 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 opticsbased traveltime 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 backwall 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 backwall 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 nonaveraged 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 subconfigurations 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. 
Goodness of the reconstruction of the interior details.
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 20core processor (Intel i910900X, 3.70GHz) and 128 GB RAM. The inversion process consisted of the BA and backpropagation 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 permittivitycontrolled plastic asteroid analogue model based on a large spatiotemporal set of experimental quasimonostatic microwave measurements. The focus was on combining a fullwavefield 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 opticsbased 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 pathlengthbased 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 opticsbased 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 backwall echo in a datadriven 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 backwall 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 nonaveraged 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 opticsbased 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 nonlinearity 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 datadriven method could best perform with symmetrical and nearspherical 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 backwall 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 backwall echo, which can mask smaller amplitude echoes from the interior of the target. Branching and backwall 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, fullwavefield 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 timedomain 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 opticsbased temporal data filtering and fullwavefield modelling applied here has been investigated in various studies such as Mocker et al. (2015). For the above justification, an advanced raytracing model (Gassot et al. 2020) applied for a comet in Herique et al. (2019a) might perform sufficiently in the inversion stage. A raytracing model to be used in connection with rough surfaces has been proposed in Cocheril & Vauzelle (2007). Both fullwavefield and raytracing 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 Cutview 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 Threedimensional 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 highperformance computing to account for the complexities of surface and internal scattering in a realistic target object. The significant backwall 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 backpropagation 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 timedomain 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 geometryrelated artefacts caused by the bilobed 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 fullwave 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 Microondes (CCRM) to use its fully equipped anechoic chamber. Centre de Calcul Intensif d'AixMarseille 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 20182025, 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
 Asphaug, E., Ryan, E. V., & Zuber, M. T. 2002, Asteroids III, 1, 463 [Google Scholar]
 Asphaug, E., Scheeres, D., & Safaeinili, A. 2008, 37th COSPAR Scientific Assembly, 1320 July 2008, Montréal, Canada, 37, 139 [Google Scholar]
 Bancel, T., Houdouin, A., Annic, P., et al. 2021, IEEE Transactions on Ultrasonics, Ferroelectrics, and Frequency Control, 68, 2554 [CrossRef] [Google Scholar]
 Britt, D. T., & Lebofsky, L. A. 1997, Asteroid: Compositional Structure and Taxonomy (Dordrecht: Springer Netherlands), 33 [Google Scholar]
 Bucci, O. M., & Isernia, T. 1997, Radio Sci., 32, 2123 [CrossRef] [Google Scholar]
 Carry, B. 2012, Planet. Space Sci., 73, 98 [CrossRef] [Google Scholar]
 Chung, K. L., & Zhong, K. 2001, A Course in Probability Theory (Academic press) [Google Scholar]
 Cocheril, Y., & Vauzelle, R. 2007, Prog. Electromagnet. Res., 75, 357 [CrossRef] [Google Scholar]
 Deng, J., Kofman, W., Zhu, P., et al. 2021a, 15th Europlanet Science Congress 2021, held virtually, 1324 September 2021, online at https://www.epsc2021.eu/, EPSC202150 [Google Scholar]
 Deng, J., Rogez, Y., Zhu, P., et al. 2021b, Comput. Phys. Commun., 265, 108002 [NASA ADS] [CrossRef] [Google Scholar]
 Dufaure, A., Eyraud, C., Sorsa, L.I., et al. 2023, A&A, 674, A72 (Paper I) [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Eyraud, C., Sorsa, L.I., Geffrin, J.M., et al. 2020, A&A, 643, A68 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Fujiwara, A., Kawaguchi, J., Yeomans, D., et al. 2006, Science, 312, 1330 [NASA ADS] [CrossRef] [Google Scholar]
 Gaskell, R., Saito, J., Ishiguro, M., et al. 2008, NASA Planetary Data System HAYAAMICA5ITOKAWASHAPEV1.0. [Google Scholar]
 Gassot, O., Hérique, A., Rogez, Y., et al. 2020, in 2020 IEEE Radar Conference (RadarConf 20), 1 [Google Scholar]
 Geffrin, J.M., Eyraud, C., Litman, A., & Sabouroux, P. 2009, Radio Sci., 44, RS003837 [Google Scholar]
 Herique, A., Agnus, B., Asphaug, E., et al. 2018, Adv. Space Res., 62, 2141 [CrossRef] [Google Scholar]
 Herique, A., Kofman, W., Zine, S., et al. 2019a, A&A, 630, A6 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Herique, A., Plettemeier, D., Kofman, W., et al. 2019b, in EPSC Abstracts, 13, EPSCDPS 2019807 [Google Scholar]
 Herique, A., Plettemeier, D., Kofman, W., Rogez, Y., & Goldberg, H. 2020, in EGU General Assembly Conference Abstracts, 19957 [Google Scholar]
 Jutzi, M., & Benz, W. 2017, A&A, 597, A62 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Kofman, W. 2012, in 2012 19th International Conference on Microwaves, Radar Wireless Communications, 2, 409 [Google Scholar]
 Kofman, W., Herique, A., Goutail, J.P., et al. 2007, Space Sci. Rev., 128, 413 [CrossRef] [Google Scholar]
 Kofman, W., Herique, A., Barbin, Y., et al. 2015, Science, 349, aab0639 [Google Scholar]
 Liu, J. S. 2008, Monte Carlo strategies in Scientific Computing (Springer Science & Business Media) [Google Scholar]
 Mocker, M. S., Schiller, M., Brem, R., et al. 2015, in 9th European Conference on Antennas and Propagation (EuCAP), IEEE, 1 [Google Scholar]
 Pursiainen, S. 2008, J. Inv. IllPosed Probl., 16, 873 [Google Scholar]
 Rivkin, A. S., Chabot, N. L., Stickle, A. M., et al. 2021, Planet. Sci. J., 2, 173 [NASA ADS] [CrossRef] [Google Scholar]
 Sava, P., & Asphaug, E. 2018a, Adv. Space Res., 62, 1146 [NASA ADS] [CrossRef] [Google Scholar]
 Sava, P., & Asphaug, E. 2018b, Adv. Space Res., 61, 2198 [NASA ADS] [CrossRef] [Google Scholar]
 Sorsa, L.I., Takala, M., Bambach, P., et al. 2019, ApJ, 872, 44 [NASA ADS] [CrossRef] [Google Scholar]
 Sorsa, L.I., Takala, M., Eyraud, C., & Pursiainen, S. 2020, IEEE Trans. Comput. Imaging, 6, 579 [CrossRef] [Google Scholar]
 Sorsa, L.I., Eyraud, C., Hérique, A., et al. 2021a, Mater. Des., 198, 109364 [CrossRef] [Google Scholar]
 Sorsa, L.I., Pursiainen, S., & Eyraud, C. 2021b, A&A, 645, A73 [NASA ADS] [CrossRef] [EDP Sciences] [Google Scholar]
 Yusuf, Y. O., Dufaure, A., Sorsa, L.I., Eyraud, C., & Pursiainen, S. 2022, Icarus, 387, 115173 [NASA ADS] [CrossRef] [Google Scholar]
All Tables
Realsize asteroid measurement parameters corresponding to the laboratory parameters.
Spatial configurations of transmissionmeasurement point pairs used in this study.
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 nonconvex 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 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 fullwave 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): Transmitterreceiver pairs of the 92 measurement points which are used in the 92/1 and 10/10 data sets. (b): the 363 transmitterreceiver pairs forming the 363/1 and 10/40 sets. (c): example of a sparse 10pointpair subconfiguration. The different colours in (a) and (b) indicate the subsystems 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 twoway 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 datadriven maximum peak identification (bottom) techniques. The topographies are shown for the 92point system (left), and the 363point 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 inphase and quadrature components of the QAMmodulated fullwave 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 subconfigurations 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 Cutview 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 Threedimensional 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 
Current usage metrics show cumulative count of Article Views (fulltext article views including HTML views, PDF and ePub downloads, according to the available data) and Abstracts Views on Vision4Press platform.
Data correspond to usage on the plateform after 2015. The current usage metrics is available 4896 hours after online publication and is updated daily on week days.
Initial download of the metrics may take a while.