Bistatic full-wave radar tomography detects deep interior voids, cracks and boulders in a rubble-pile asteroid model

TL;DR

This study demonstrates that sparse bistatic full-wave radar tomography can effectively detect internal features like voids, cracks, and boulders in a complex rubble-pile asteroid model, enhancing asteroid interior imaging capabilities.

Contribution

It shows that bistatic CRT can identify internal asteroid features with sparse measurements and noise robustness, advancing asteroid interior exploration methods.

Findings

Detects surface layers, voids, and boulders with -10 dB noise level

Bistatic setup improves inversion robustness over monostatic

Reconstructing smooth backgrounds remains challenging, suggesting need for complementary techniques

Abstract

In this paper, we investigate full-wave computed radar tomography (CRT) using a rubble-pile asteroid model in which a realistic shape (Itokawa) is coupled with a synthetic material composition and structure model. The aim is to show that sparse bistatic radar measurements can distinguish details inside a complex-structured rubble-pile asteroid. The results obtained suggest that distinct local permittivity distribution changes such as surface layers, voids, low-permittivity anomalies, high-permittivity boulders, and cracks can be detected with bistatic CRT, when the total noise level in the data is around -10 dB with respect to the signal amplitude. Moreover, the bistatic measurement set-up improves the robustness of the inversion compared to the monostatic case. Reconstructing the smooth Gaussian background distribution was found to be difficult with the present approach, suggesting that complementary techniques, such as gravimetry, might be needed to improve the reliability of the inference in practice.

Tables

Table 1. Table 1: Details in the asteroid models (A)–(E)

				Detail sizeaaSizes are indicated in meters, with respect to the wavelengths of a 20 MHz center frequency signal, and with respect to the wavelengths of a 2 MHz bandwidth signal.
Model	Description	Detail	${\epsilon }_r$	Length	Width	Depth
(A)	Single void	Ellipsoid	1	120 m	120 m	55 m
				8.0$\lambda$	8.0$\lambda$	3.7$\lambda$
				0.8$\lambda$	0.8$\lambda$	0.4$\lambda$
(B)	Highly porous inclusion	Ellipsoid	2	120 m	120 m	55 m
				11.3$\lambda$	11.3$\lambda$	5.2$\lambda$
				1.1$\lambda$	1.1$\lambda$	0.5$\lambda$
(C)	High-permittivity boulder	Ellipsoid	15	120 m	120 m	55 m
				31.8$\lambda$	31.8$\lambda$	14.1$\lambda$
				3.1$\lambda$	3.1$\lambda$	1.4$\lambda$
(D)	Deep crack	Crack	1	185 m	10–55 m	10–50 m
				12.3$\lambda$	0.7–3.7$\lambda$	0.7–3.3$\lambda$
				1.2$\lambda$	0.07–0.4$\lambda$	0.07–0.3$\lambda$
(E)	Shallow crack	Crack	1	185 m	10–55 m	10–50 m
				12.3$\lambda$	0.7–3.7$\lambda$	0.7–3.3$\lambda$
				1.2$\lambda$	0.07–0.4$\lambda$	0.07–0.3$\lambda$

Table 2. Table 2 : Distinguishability of model details

	Surface layer		Interior
Model	DataaaBased on the data curve in Figure 4(a).	Reconst.bbBased on the reconstructions in Figures 5 and 6.	DataaaBased on the data curve in Figure 4(a).	Reconst.bbBased on the reconstructions in Figures 5 and 6.
(A)	Clear	Moderate	Clear	Clear
(B)	Clear	Moderate	Clear	Moderate
(C)	Clear	Moderate	Moderate	Moderate
(D)	Clear	Moderate	Weak	Weak
(E)	Clear	Moderate	Weak	Moderate

Table 3. Table 3: Mean square errors (MSE), mean absolute errors (MAE), and relative mean absolute errors (MAE-R)

Model	Error	Detail	Surface layer	Global
(A)	MSE	8.21	0.99	1.19
	MAE	2.80	0.89	0.93
	MAE-RaaNormalized relative to the difference between the detail permittivity and the global mean of the Gaussian random field (4).	0.93	…	…
(B)	MSE	3.65	0.89	0.98
	MAE	1.90	0.92	0.91
	MAE-RaaNormalized relative to the difference between the detail permittivity and the global mean of the Gaussian random field (4).	0.95	…	…
(C)	MSE	124.67	0.92	3.53
	MAE	11.12	0.87	1.08
	MAE-RaaNormalized relative to the difference between the detail permittivity and the global mean of the Gaussian random field (4).	1.01	…	…
(D)	MSE	7.10	0.93	1.05
	MAE	2.49	0.89	0.90
	MAE-RaaNormalized relative to the difference between the detail permittivity and the global mean of the Gaussian random field (4).	0.83	…	…
(E)	MSE	6.86	0.90	1.04
	MAE	2.38	0.87	0.89
	MAE-RaaNormalized relative to the difference between the detail permittivity and the global mean of the Gaussian random field (4).	0.79	…	…

Equations

$$\mathbf{y}=\mathbf{Lx}+{\mathbf{y}}^{(bg)}+\mathbf{n}.$$

$${\mathbf{x}}_{\ell }=({\mathbf{L}}^T\mathbf{L}+\alpha \mathbf{D}{\mathbf{\Gamma }}_{\ell }\mathbf{D}{)}^{-1}{\mathbf{L}}^T\mathbf{y},{\mathbf{\Gamma }}_{\ell }=\mathrm{diag}(|{\mathbf{Dx}}_{\ell }|{)}^{-1},$$

$$D_{i,j}=\beta {\delta }_{i,j}+\frac{{\ell }^{(i,j)}}{{\max }_{i,j}{\ell }^{(i,j)}}(2{\delta }_{i,j}-1),$$

$${\delta }_{i,j}=\{\begin{array}{l}1,\text{ if }j=i,\\ 0\text{ otherwise. }\end{array}$$

Full text

Asteroid — Interior — Radar — Computed Radar Tomography — Rubble-Pile Asteroid

1 Introduction

This paper investigates spaceborne computed radar tomography (CRT) applied to a rubble-pile asteroid model. Spaceborne CRT inherits from the airborne ground penetrating radar (GPR) which was originally developed to serve in cost-effective surveying of the underground in applications which entail the necessity to work with an antenna not in contact with the surveyed structure (Catapano et al., 2012a, b; Soldovieri et al., 2017). Such airborne GPR systems have been validated by airborne measurement data and numerical experiments (Fu et al., 2014). Space mission concepts to perform radio tomography of small asteroids have been proposed by, for example, Asphaug et al. (2001), Safaeinili et al. (2002), Snodgrass et al. (2018), and Bambach et al. (2018).

The first attempt to reconstruct the deep interior of a small solar system body (SSSB) was the Comet Nucleus Sounding Experiment by Radio-wave Transmission (CONSERT), a part of European Space Agency’s (ESA) Rosetta mission to explore the comet 67P/Churyumov-Gerasimenko, in which a radio signal was transmitted between the orbiter Rosetta and the lander Philae (Kofman et al., 1998, 2007, 2015). Many more missions aiming at exploring the structure and composition of SSSBs are currently ongoing or being planned. In August 2018, the Osiris-REx mission by NASA (Berry et al., 2013) begun its rendezvous with the asteroid 101955 Bennu to measure its physical, geological and chemical properties and collect a sample of the asteroid surface regolith (Lauretta et al., 2017). The Hayabusa mission by the Japanese Aerospace Exploration Agency (JAXA) explored the asteroid Itokawa already in 2005 (Japan Aerospace Exploration Agency, 2003; Okada et al., 2006; Nakamura et al., 2011; Tsuchiyama et al., 2011) and retrieved surface regolith for analysis on Earth (Fujiwara et al., 2006), confirming the Itokawa’s Earth-based classification as an S-type asteroid. The images taken by Hayabusa on Itokawa have been used to analyze the size distribution of boulders on its surface. The results indicate that the boulders cannot solely be the product of cratering but that many of them originated from the disruption of a larger body (Saito et al., 2006; Michikami et al., 2008). JAXA’s mission Hayabusa2 arrived at the asteroid 16173 Ryugu in mid-June 2018. It will survey its target for a year and a half, returning back to Earth in December 2020 (Tsuda et al., 2013). These missions by NASA and JAXA concentrate on the surface properties of the target asteroids and do not carry instruments which could be used to explore the deep interior structures.

The next candidate to deploy a CRT system to explore the deep interior of a SSSB is ESA’s asteroid mission Hera which targets the binary near-Earth asteroid system Didymos and is due to launch in 2023 (Michel et al., 2016). The current plan allows deployment of a 6U CubeSat form factor payload (Carnelli, 2018) which has been recently proposed for the radar-carrying spacecraft concept DISCUS (Deep Interior Scanning CUbeSat) (Bambach et al., 2018). Furthermore, ESA’s recent Concurrent Design Facility (CDF) study report (European Space Agency, 2018) suggests that small planetary platforms involving CubeSats can provide future opportunities for CRT.

In this paper, CRT is applied to a rubble-pile asteroid model utilizing the mathematical approach proposed in Pursiainen & Kaasalainen (2016) and Takala et al. (2018b). We investigate a synthetic framework featuring the radar specifications of the DISCUS concept, and a target with a complex shape, structure and material composition. The shape of the asteroid Itokawa is coupled with a Gaussian random field model for the relative dielectric permittivity distribution together with added structural details. We base our model on the recent observations and impact simulations which suggest that the internal porosity of the asteroid body varies, increasing towards the center, and that it may have a detailed structure (Carry, 2012; Deller, 2017; Jutzi & Benz, 2017).

The aim is to validate numerically the mathematical approach of Pursiainen & Kaasalainen (2016) and Takala et al. (2018b) and also to explore, whether sparse bistatic radar measurements can distinguish details inside a complex-structured rubble pile asteroid.

2 Materials and Methods

In this study, we utilize the finite element time-domain (FETD) approach (Pursiainen & Kaasalainen, 2016) equipped with a far-field model proposed in Takala et al. (2018b). In FETD, the radar signal is propagated over a suitably chosen time interval within a volumetric finite element (FE) mesh which can be adapted accurately to a given set of surface features. The permittivity structure can be found via a multigrid approach (Takala et al., 2018a). That is, the signal is inverted using a coarser grid which is a nested strucure with respect to the more refined FE mesh.

2.1 Model of the asteroid Itokawa

The FE mesh was constructed based on the detailed shape model of the asteroid Itokawa (Figure 1, left) which is openly available as a triangular stereolitography (STL) mesh (Hayabusa Project Science Data Archive, 2007). The scaling of the asteroid was assumed to match with the actual size of Itokawa, whose longest diameter is 535 meters (Fujiwara et al., 2006). The unstructured triangulated asteroid surface was imported to Meshlab (Cignoni et al., 2008) and resampled with a Poisson-disk sampling algorithm to obtain a uniform mesh with a desired resolution of 5 762 nodes and 11 520 faces. The volumetric asteroid model was divided into a surface and interior compartment. The boundary between these two parts was obtained by scaling a resampled and smoothed-out version of the outer surface by the factor of $0.9$. An additional ellipsoidal compartment was utilized as a void detail, when simulating the measurement data.

2.2 The structure and the 3D domain structure

The computation domain $\Omega$, depicted in Figure 1 (second from the left) was constructed in the Gmsh software (Geuzaine & Remacle, 2009), in which the structure of the computation domain $\Omega$ can be described within a geometry (GEO) file and the compartments of the asteroid model can be obtained via their STL surface grids. Two coarse FE meshes with 38k nodes and 212k tetrahedra were generated using Gmsh. One of these was employed for simulating the measurements, and the other one, not including any interior details, was used in the inversion stage. The reason to apply two different grids was to avoid the so-called inverse crime, i.e., overly good match between the actual permittivity distribution and its reconstruction. The permittivity was modeled as a piecewise constant function within the coarse grid. A nested two times uniformly refined grid was utilized, when simulating the signal propagation, in order to achieve a sufficient resolution with respect to the wavelength.

2.3 The set of measurement points

The measurement point set was modeled after the DISCUS mission concept (Takala et al., 2018b) in which a master and slave CubeSat both equipped with a half-wavelength dipole antenna orbit the target asteroid at a 5 km distance. The angle between the orbiters is 25 degrees with respect to the asteroid’s center of mass (Figure 1, third from the left). The measurement configuration is bistatic: the master CubeSat both transmits and records the signal, and the slave provides an additional receiver. The angle $\theta$ between the measurement plane normal and the asteroid spin, determining the angular coverage of the measurements, is assumed to be $70$ degrees. Optimally, a full coverage could be obtained with $\theta=90^{\circ}$. A practically obtainable orbiting direction may be expected to have a somewhat but not very much lower value, as the target asteroid will be likely to have a close-to retrograde spin (La Spina et al., 2004). In total, 64 measurement points were included in the point cloud. The resulting bistatic set of measurement points and antenna orientations has a sparse and limited-angle spatial coverage with an aperture around the z-axis (Figure 1, right).

2.4 Radar specifications

Inheriting from the DISCUS concept (Bambach et al., 2018), the radar is assumed to have a transmission power of 10 W, a 2 MHz total signal bandwidth, and a relatively low 20 MHz center frequency. In practice, $\leq 100$ MHz will be necessary in order to achieve appropriate signal penetration and to minimize solar noise (Kofman, 2012). The maximal range (imaging) resolution following from these parameters is about 35-40 m inside an asteroid with relative permittivity approximately 4.

2.5 Computing Gaussian random field inside a given asteroid geometry

A Gaussian random field for the dielectric relative permittivity, $\varepsilon_{r}$, was first generated in a regular 20-by-20-by-20 lattice (Figure 2, left). The mean of the random field was set to $\varepsilon_{r}=4$, assuming a standard deviation equal to one. The correlation value between the adjacent lattice points, was set to $0.2$ based on visual examination. The correlation was chosen to be isotropic. These choices were thought to roughly account for the recent impact simulation results in which large pieces of porous rubble, bound together by gravity, are concentrated in the interior (Jutzi & Benz, 2017), yet allowing some clear randomness of the distribution due to the fact that the exact interior structure is mainly unknown (Carry, 2012).

Implementation of the Gaussian random field algorithm was retrieved from Constantine (2012). The asteroid interior geometry was then fitted inside the cubic mesh and the values were assigned to the finite element mesh nodes via the nearest neighbour interpolation principle. Each finite element mesh tetrahedron was assigned a value which is the average of the electrical permittivities of each node constituting the element. Finally, the random field was restrictred to permittivity values greater or equal to one, which is the dielectric permittivity of vacuum and, therefore, the physical minimum of relative permittivity.

2.6 Structural details

This section describes the details of the five asteroid models (A)–(E) applied in this study. The details of these have been given also in Table 1. In each model, the permittivity distribution in the interior part is determined by a Gaussian random field with values approximately between two and six (Figure 2, right). Furthermore, each one includes a surface layer with a lower permittivity to model the higher porosity near the surface predicted in the impact simulations (Jutzi & Benz, 2017). In each case, a different random realization of the Gaussian random field is used.

2.6.1 Model (A): A single void

In (A), to investigate anomaly detection, a single ellipsoid with a relative permittivity value $\varepsilon_{r}=1$ (vacuum) was placed deep in the interior (Figures 3(a) and 3(b)). This anomaly could be, for example, a void between two large boulders constituting a part of the body of the asteroid.

2.6.2 Model (B): A Highly porous inclusion

In (B), we study a situation in which the anomaly is not as distinct from the surroundings as a void space (Figure 3(c)). That is, the interior may include areas which are occupied by very porous or otherwise low permittivity material. Therefore, we investigate also an ellipsoid, whose relative dielectric permittivity is adjusted to $\varepsilon_{r}=2$ instead that of of the vacuum ($\varepsilon_{r}=1$).

2.6.3 Model (C): A High-permittivity boulder

An asteroid may also be formed around a single large boulder with dielectric permittivity significantly higher than that of its surroundings. Such system is investigated with model (C) in which the dielectric permittivity of the enclosed ellipsoid is assigned with the value $\varepsilon_{r}=15$ which is significantly higher than the surrounding other material (Figure 3(d)). This value corresponds, e.g., to solid basalt (Hansen et al., 1973), and was included in the study to investigate whether such permittivity structures can be reconstructed and to provide a more complete view of the capabilities of radar tomography.

2.6.4 Model (D): A Deep crack

In (D), to model a random crack inside the asteroid, a simple nearest-neighbour walking algorithm was applied within the set of FE mesh nodes. Each node included in the path created by the random walk was assigned with the permittivity value of one (vacuum), and the permittivity of each tetrahedron in the FE mesh was determined by the average permittivity of its nodes.

The resulting crack is irregularly shaped and smaller in comparison to the void (Figures 3(e) and 3(f)). The width of the crack is in the order of 10-50 meters, depending on the location within the asteroid.

2.6.5 Model (E): A Shallow crack

In (E), the random walk approach was utilized to create a shallow crack. This crack was also allowed to bifurcate to form a two-branched structure (Figures 3(g) and 3(h)).

2.7 Conductivity distribution

The conductivity distribution causing signal attenuation was assumed to be a latent parameter, i.e., an uninteresting nuisance parameter, depending on the permittivity according to $\sigma=5\xi\varepsilon_{r}$, where $\xi=(\mu_{0}/\varepsilon_{0})^{-1/2}$ s*-1* with spatial scaling parameter $s=2100$ m (Takala et al., 2018b). If $\varepsilon_{r}=4$, this results in a conductivity of around $2.5\times 10^{-5}$ S m*-1*, which matches roughly with an attenuation rate of 25 dB km*-1* approximating a low-frequency signal decay within a porous body (Kofman, 2012).

2.8 Forward model

The inversion process requires a reference permittivity $\varepsilon^{(bg)}_{r}$, a background model. Such a model is constructed similar to the detailed, structural models with the difference that the interior permittivity is constant, here $\varepsilon^{(bg)}_{r}=4$. We utilize a linearized forward model given by Pursiainen & Kaasalainen (2016):

[TABLE]

Here, the vectors $\mathbf{y}$ and $\mathbf{y}^{(bg)}$ contain the measured and simulated data for the actual permittivity distribution $\varepsilon_{r}$ and its constant-valued approximation $\varepsilon_{r}^{(bg)}$, respectively. The vector $\mathbf{x}$ determines the discretizised permittivity distribution, $\mathbf{L}$ denotes the Jacobian matrix resulting from the linearization at $\varepsilon_{r}^{(bg)}$, and $\mathbf{n}$ contains both the measurement and forward modeling errors.

The noise vector $\mathbf{n}$ is here assumed to be a zero-mean Gaussian random variable due to the various potential and unknown noise formation processes, e.g., signal attenuation due to latent factors. Following from the rough estimates for the cosmic background noise level with the present radar and mission specifications (Takala et al., 2018b), the standard deviation was chosen to be -15 dB with respect to the maximal entry of the difference $|\mathbf{y}-\mathbf{y}^{(bg)}|$ between the measured and simulated signal.

2.9 Inversion approach

We use a simple total variation (TV) based iterative regularization approach presented in Pursiainen & Kaasalainen (2016). An estimate of $\mathbf{x}$ in equation 1 can be produced via the iterative regularization procedure

[TABLE]

in which $\bm{\Gamma}_{0}=\mathbf{I}$ and $\mathbf{D}$ is of the form:

[TABLE]

The first term is a weighed identity operator limiting the total magnitude of $\mathbf{x}$, whereas the second term penalizes the jumps of $\mathbf{x}$ over the edges of the mesh, the $\ell^{(i,j)}$ being the edge length.

The inversion process in the equation 2 minimizes the function $F(\mathbf{x})=\|\mathbf{Lx}-\mathbf{y}_{bg}-\mathbf{y}\|_{2}^{2}+2\sqrt{\alpha}\|\mathbf{Dx}\|_{1},$ in which the second term equals the total variation of $\mathbf{x}$ if $\beta=0$. The coefficients $\alpha$ and $\beta$ are the reqularization parameters.

The total variation penalty function is evaluated with respect to the tetrahedral inversion mesh as shown in Takala et al. (2018b). The inversion computation was run with regularization parameters values $\alpha=0.01$ and $\beta=\{0.005,0.0075\}$ depending on the model. The parameters were adjusted based on preliminary experiments with the goal to maximize the distinguishability of the interior details and to obtain an appropriate range of permittivity values with respect to the actual distribution.

3 Results

The full-wave computed radar tomography simulation was performed for five different model structures (A)–(E) covering the following features: (A) a void space, (B) a highly porous inclusion, (C) a high-permittivity boulder, (D) a deep crack, and (E) a shallow crack inside the asteroid. The results of the simulations can be found in the Figures 4–7 and the Table 2. The results suggest that deep interior permittivity anomalies can be detected inside a rubble-pile asteroid with limited-angle set of bistatic full-wave data. Nevertheless, detecting smooth structures was found to be difficult in comparison to how the void, boulder and crack details were distinguished. The model-specific results for (A)–(E) are described in the below sections.

3.1 Simulated measurement data

Time series of the normalized data at a single measurement point is shown in the Figure 4. The shown time interval 1.3–4.3 $\mu$s shows signal propagation within the asteroid. Extending the duration beyond 4.3 $\mu$s was found to be unfeasible due to the noise effects caused by the inward-directed reflection peak from the outer surface facing away from the measurement position. The first echo at time point 2 $\mu$s originates from the surface of the asteroid. In the background model (red line), the permittivity of the asteroid interior is constant $\varepsilon_{r}=4$ across the whole asteroid interior and hence a more distinct echo is recorded in comparison to each the detailed model (blue line) which includes a surface layer with a lower permittivity value ($\varepsilon_{r}=3$). Hence, also the amplitude of the reflected signal is lower in the latter case.

The figures 4(a) and 4(b) show how ellipsoids and cracks, respectively, can be distinguished in the raw signal at a single measurement point. The differences in the curves are due to the differences in how the signal propagates through the asteroid. The non-linearity of the radar signal propagation is revealed by a mutual comparison between the amplitudes of the signals for models (A), (B), and (C). The geometrical models of these are identical, but the signal curves, however, differ also regarding the shape of the measured signal. The differences between the signals for the deep and shallow crack (D) and (E) (Figure 4(b)) are minor but discernible. The shallow crack produces a more distinct echo which can be observed at the time point of approximately 2.7 $\mu$s.

3.2 Model (A): Void detection

A comparison between the reconstructed permittivity and the exact 3D distribution for model (A) is shown in the Figures 5(a) and 5(b), respectively. It is evident that the (vacuum) void space of (A) is detectable. The quality of the radial (depth) accuracy with respect to the asteroid’s body was observed to be better than in the tangential one. This is obviously due to the sparse spatial distribution of the measurement points compared to the relatively dense time-resolution of the data which results in reliable localization capability.

3.3 Model (B): Detection of a highly porous inclusion

The task of detecting highly porous materials ($\varepsilon_{r}=2$), case (B), was also found to be feasible. However, the reconstruction quality obtained for (B) was weaker compared to that of (A) (Figure 5(e)). The exact shape of the inclusion found for (B) is not as obvious as in the case of (A), but the location of the anomaly is evident as shown by the Figures 5(e) and 5(h).

3.4 Model (C): Boulder detection

The Figures 5(g) and 5(i) show that a high-permittivity boulder can be detected. Similar to the void detection, also here the best detection accuracy was obtained in the transversal direction. The surface layer in the reconstruction (Figure 5(i)) is less obvious than for the cases of (A) and (B).

3.5 Model (D): Deep crack detection

The reconstruction for model (D) shown on the left in the Figure 6 suggests that the the deep crack is practically invisible for the radar. Obviously, its location in the deep interior, and the limited-angle data, i.e., the absence of the measurement points around the z-axis, are potential factors limiting the detectability of the crack.

3.6 Model (E): Shallow crack detection

In contrast to the deep crack, a moderately clear echo was detected for the shallow crack in the case (E) as shown on the right in the Figure 6. Based on the reconstruction (Figure 6(d)), the complex shape can also be discerned. However, the exact shape and, especially, the branches of the crack cannot be accurately distinguished.

3.7 Goodness of reconstructions

The qualitative assessment of the goodness of the reconstructions presented in the Figures 5 and 6 is shown in the Table 2. The quantitative goodness of a reconstruction is measured by the mean square error (MSE) and the mean absolute error (MAE) computed separately for the anomaly detail and the surface layer, in addition to the global reconstruction containing the detail, the surface layer and the remaining asteroid interior of the asteroid model. The relative mean absolute error (MAE-R) is computed for the details only to compare the goodness of the detail reconstruction between the models (A)–(E).

The computed mean and absolute error values (Table 3) reflect the qualitative reconstructions. The errors are greatest in the detail areas and smallest in the surface layer. As is evident in the reconstruction images in the Figures 5 and 6, the majority of the volumes in the ellipsoids and the cracks is not captured by the reconstruction and hence the permittivity values of these areas differ from the exact model very clearly. This is especially the case with the boulder model (C), in which the difference between the permittivity values in the exact model and the reconstruction varies greatly.

3.8 Comparison of bistatic to monostatic measurement

Model (A) was used to compare reconstructions obtained with monostatic (single-satellite) and bistatic measurement data. The results show (Figures 7(a) and 7(b)) that both satellite configurations enable the detection of the void. The bistatic measurement appears to show slightly more prominent shape of the anomaly, but this may not be significant in practical applications. Overall, bistatic measurement approach appears to provide robustness to the reconstruction process. The same finding was previously reported and quantified in Takala et al. (2018b) and can be accounted for by the fact that, in addition to the second receiver, the bistatic measurement set-up includes also the monostatic transmitter-receiver data collection. Hence, the bistatic measurement includes more information on the object.

3.9 Effect of noise

Finally, noise levels between -25 and 0 dB were investigated to determine the effect of noise on the reconstruction accuracy. Based on the results (Figures 7(c)-7(h)), the artifacts appear at the -8 dB level after which they increase rapidly along with the noise.

4 Discussion

The aim of this paper was to numerically validate the full-wave computed radar tomography (CRT) model developed by Pursiainen & Kaasalainen (2016) and Takala et al. (2018b) for complex-structured asteroid models and especially to estimate the reconstruction capability of the CRT for rubble-pile asteroids. The surface model of the asteroid Itokawa was used because it is one of the few asteroids which the surface and bulk properties are known (Abe et al., 2006; Fujiwara et al., 2006; Okada et al., 2006; Barnouin-Jha et al., 2008; Nakamura et al., 2011; Tsuchiyama et al., 2011), and for which a high-resolution surface model is available (Hayabusa Project Science Data Archive, 2007). The structural models (A)–(E) of this study were generated by combining the data and knowledge available on the physical properties of small solar system bodies and on how electromagnetic waves behave in such media (Kofman, 2012; Kofman et al., 2015; Herique et al., 2018).

To create a complex model for the deep interior part, roughly matching with the mass-concentration estimates obtained from the asteroid impact evolution studies (Jutzi & Benz, 2017), the dielectric properties of the asteroid were modelled by a Gaussian random field, which is a well-established approach to modelling random spatial structures in geostatistics and earth sciences (Lantuéjoul, 2002; Christakos, 2005; Chilès & P, 2012). Gaussian random fields are also used to generate random 3D porous structures (Roberts & Garboczi, 2002) and modelling concrete (Bićanić et al., 2010) in civil engineering applications. Models (A)–(E) were obtained adding a surface layer together with deep interior details to a Gaussian random field. Special interest was paid on modelling the inhomogeneous deep interior permittivity and void space accounting for bulk macroporosity of the body, which can be expected to be, for example, as high as 41 percent in the rubble-pile Itokawa (Fujiwara et al., 2006).

The results suggest that permittivity anomaly details can be detected using CRT within a complex-structured asteroid, excluding the deep areas around the center of mass which are challenging due to signal attenuation and the limited-angle measurement point distribution with an aperture around the z-axis. However, it also appears that the smooth Gaussian random field structure is difficult to reconstruct accurately. Compared to the earlier work (Pursiainen & Kaasalainen, 2016; Takala et al., 2018b), the present inversion scheme is more realistic, showing that the practically achievable reconstruction quality will largely depend on the shape and the dielectric properties of the asteroid body. Moreover, it is obvious that the structures closer to the surface are detected more clearly than those lying deep inside the body.

Reconstructing the interior details was shown to be feasible up to approximately -10 dB of noise. For a carrier frequency of 20 MHz, the galactic noise can be estimated to be around $5\times 10^{-20}$ W m*-2Hz-1*. At a distance of 1 AU, the dominating measurement noise is likely to originate from the Sun which radiates with the magnitude of $5\times 10^{-19}$ W m*-2Hz-1* and $4\times 10^{-23}$ W m*-2Hz-1* for its active and inactive (quiet) phase of sunspot activity, that is, for surface temperatures $10^{6}$ and $10^{10}$ K, respectively (Kraus, 1967; Barron et al., 1985). During the active phase there are radio emissions in time scales varying from seconds to hours. This noise can be diminished, by some amount, via orienting the antenna in an optimal way. When moving away from the Sun, the radiation, in general, decreases proportionally to the inverse square of the radius. This means that at 3 and 30 AU, the radiation levels due to the Sun will be -19 and -59 dB, respectively, with respect to those experienced at 1 AU. The latter one of these values is a rough estimate for the Kuiper belt.

The imaging resolution was found to be higher in the radial direction than in the tangential one. The details recovered also seem to be elongated in the tangential direction, exhibiting similar wave-front properties as conventional radargrams (Daniels & of Electrical Engineers, 2004). The elongation can be interpreted to be due to the sparsity of the spatial measurement points. It is also present in other full-wave applications of CRT (Gueting et al., 2015, 2017), suggesting that enhancing the results significantly without changing the radar specifications might be difficult. Consequently, it seems that exact shapes cannot be accurately reconstructed, for example, structures oriented towards the center of the asteroid. This is also reflected in the mean square errors and mean absolute errors computed for the reconstructions. However, based on the crack detection results it seems also reasonable to assume that even complex structures near the surface can be detected. The quality of reconstruction and the exact permittivity values may also be assumed to depend on the applied inversion technique and regularization parameters. In the present study, some of the details, most prominently the surface layer, are more obvious in the data (Figure 4) than in the final reconstructions (Figures 5 and 6), suggesting that the inversion quality might be improved for example via depth weighting or noise-robust techniques such as multigrid, Markov chain Monte Carlo sampling or expectation maximization techniques (Liu, 2008; Tilley et al., 2017; Takala et al., 2018a).

Other future work on this topic will be to investigate methods to utilize complementary data in the inversion process. One way is to use multiple radar frequencies as suggested by Herique et al. (2018). For example, a frequency range of 25, 50 and 100 MHz can be considered to provide optimization of radar penetration depth. However, because the radar signal attenuation is directly proportional to signal frequency (Kofman, 2012), the scattering and attenuation effects are stronger in the higher end of such a frequency range. In the present Itokawa model, extra signal attenuation in the deep part at the ellipsoidal void detail is 10 to 30 dB in the center frequency range of 50 to 100 MHz (Kofman, 2012), respectively. Therefore, due to the lower penetration capability, it is possible that 50 and 100 MHz frequencies are suitable to exploring only the surface layer of asteroids and the deep interior details can only be recovered by low-frequency waves. The analysis in this work involved only inversion of the modulated signal across the whole bandwidth. By combining the data from the modulated signal and the envelope, the accuracy of the method might be further improved, for example in improving crack detection. Moreover, gravimetry combined with radar sensing might help to improve the visibility of smooth structures or layers (Fregoso & Gallardo, 2009; Park et al., 2014). Resolving the effect of surface scattering from anomalies such as large boulders, steep hills or crevasses on the quality of the inversion, and enhancing the noise model to distinguish different error sources, will also be important future work direction.

5 Conclusion

The results of this paper show that bistatic computed radar tomography can detect deep interior voids, cracks and high permittivity boulders inside a complex rubble pile asteroid model when the total noise level in the data is around -10 dB with respect to the signal amplitude. The results suggest that permittivity anomaly details can be detected within a complex-structured asteroid by using CRT, excluding the deep areas around the center of mass which are challenging due to signal attenuation and the limited-angle measurement point distribution with an aperture around the z-axis.

Furthermore, the bistatic measurement set-up improves the robustness of the inversion compared to the monostatic case. However, recovering the smooth Gaussian background distribution was found to be difficult with the present approach, suggesting that complementary techniques, such as gravimetry, might be needed to improve the reliability of the inference in practice. The inversion quality might be improved via method development and/or analyzing complementary data, such as a higher-frequency signal, carrier wave envelope or gravity field measurements.

6 Acknowledgments

L.-I.S., M.T. and S.P. were supported by the Academy of Finland Centre of Excellence in Inverse Modelling and Imaging.

Appendix A Computational framework and performance

For inverting the data, we used the computational full-wave radar tomography approach developed in Pursiainen & Kaasalainen (2016), Takala et al. (2018a), and Takala et al. (2018b). The computationally intensive part of the simulation was run in the GPU partition of Tampere Center for Scientific Computing (TCSC) Narvi cluster which consists of 8 GPU nodes with 20 CPU cores, 4 GPU each, totalling 32 NVIDIA Tesla P100 16 GB GPUs. The inversion reconstruction procedures were run on a Lenovo P910 workstation equipped with two Intel Xeon E5 2697A $2.6$ GHz 16-core processors and 128 GB RAM.

The creation of the system from node and tetrahedra data, including all parts of the asteroid geometry, orbiter points and interpolation of the finite element mesh required approximately 11 GB of memory space and 6 700 seconds (approximately 1 hrs 52 minutes). Out of this, the actual system creation took 652 seconds (approximately 10 minutes) and creating the spatial interpolation matrix accounted for the rest of the computational time.

The minimum computation time of one transmitted signal was 11 400 seconds (approximately 3 hrs 10 minutes). The system size was approximately $4.8$ GB. Because the forward simulation was parallelized into 16 processes, the total computation time of the entire simulation, 64 transmitter points, took a minimum of 45 400 seconds (approximately 12 hrs 37 minutes). A typical expected computation time for the entire forward simulation was approximately $13.5$ hours.

The present study shows that full-wave (full-bandwidth) data can be computed and inverted for a rubble-pile asteroid with a realistic size and shape using a state-of-the-art computing cluster. GPU acceleration was observed to be necessary in the forward simulation stage in order to achieve feasible computation times and a speed-up by at least a factor of 10. In radar applications, similar experiences have recently been reported, for example in (Cordua et al., 2012). While the GPU-based forward simulation approach performs well regarding the speed, it is at the moment restricted with respect to the system size. In our case, doubling the asteroid diameter would not have been possible due to the limited GPU RAM which in our cluster is currently 12 GB. Since GPUs are a rapidly progressing field of technology, the available memory capacity may be expected to be significantly larger within a few years. Further development of the present computation framework may involve a pipeline for simulating realistic asteroid interiors based on a finite element mesh structural model.