Liisa-Ida Sorsa

Mika Takala

Christelle Eyraud

Sampsa Pursiainen

A Time-Domain Multigrid Solver with Higher-Order Born Approximation for Full-Wave Radar Tomography of a Complex-Shaped Target

Keywords: Radar Tomography, Small Solar System Bodies, Finite Element Time-Domain (FETD), Multigrid, Graphics Computing Unit (GPU)

This paper introduces and evaluates numerically a multigrid solver for non-linear tomographic radar imaging. Our goal is to enable the fast and robust inversion of sparse timedomain data with a mathematical full-wave approach utilizing a higher-order Born approximation (BA). Full-wave inversion is computationally expensive, hence techniques to speed up the numerical procedures are needed. To model the wave propagation effectively, we use the finite element time-domain (FETD) method, which is equipped with a multigrid scheme to enable the rapid evaluation of the higher-order BA. As a potential application, we consider the tomography of small solar system bodies (SSSBs) and asteroid interiors particular, the latter of which can contain internal details observable by radar, e.g., layers, voids and cracks. In the numerical experiments, we investigated monostatic, bistatic and multistatic measurement configurations. The results obtained suggest that, with a relevant noise level, the tomographic reconstruction quality can be improved by applying the higher-order BA in comparison to the first-order case, which we interpret as a linearization of the inverse problem. Our open multigrid-FETD solver for Matlab (The Mathworks Inc.) is available online. It applies Matlabs features for graphics computing unit acceleration to enhance computational performance.

I. INTRODUCTION

This paper concerns tomographic full-wave radar imaging in which the internal permittivity distribution of a given domain is to be reconstructed via transmitting and measuring electromagnetic waves penetrating through a domain [START_REF] Yu | Microwave imaging in layered media: 3-d image reconstruction from experimental data[END_REF], [START_REF] Eyraud | A large 3d target with small inner details: A difficult cocktail for imaging purposes without a priori knowledge on the scatterers geometry[END_REF], [START_REF] Haynes | real-time microwave imaging of differential temperature for thermal therapy monitoring[END_REF], [START_REF] Eyraud | Influence of the description of the scattering matrix on permittivity reconstruction with a quantitative imaging procedure: polarization effects[END_REF]. We consider inverting a wave equation in the time domain for a sparse set of measurements to reconstruct the interior structure of a complex target object, for example, a small solar system body (SSSB), which can be sounded by radar in a planetary space mission [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF], [START_REF] Hérique | Direct observations of asteroid interior and regolith structure: science measurement requirements[END_REF], [START_REF] Sava | 3d radar wavefield tomography of comet interiors[END_REF], [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF]. The first attempt aiming at tomographic reconstruction of the interior of a SSSB was the Comet Nucleus Sounding Experiment by Radio-wave Transmission (CONSERT) [START_REF] Kofman | Comet nucleus sounding experiment by radiowave transmission[END_REF], [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF], [START_REF] Kofman | The comet nucleus sounding experiment by radiowave transmission (CONSERT): A short description of the instrument and of the commissioning stages[END_REF], a part of the European Space Agency's (ESA) Rosetta mission, which visited the comet 67P/Churyumov-Gerasimenko in 2014. The density estimates available today [START_REF] Carry | Density of asteroids[END_REF] suggest that SSSBs can be highly porous and contain a significant amount of void space, which might be detectable via tomographic radar imaging.

In this study, we aim at developing a mathematical and numerical higher-order Born approximation (BA) [START_REF] Verma | Elementary Theory of Scattering[END_REF], [START_REF] Born | Quantenmechanik der stoßvorgänge[END_REF] for the effective modeling and inversion of non-linear fullwave scattering in the time-domain for a complex-shaped or structured target object. Data processing in the time domain has been applied, e.g., in CONSERT [START_REF] Kofman | The comet nucleus sounding experiment by radiowave transmission (CONSERT): A short description of the instrument and of the commissioning stages[END_REF], [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF], as it is beneficial in the presence of complex scattering; the noisy parts of the signal can be filtered out of the data based on their travel time, including, e.g., the effects of anisotropic structures or highly uncertain scattering. As of other studies concentrating on the present theoretical context, analytical and computational approaches to solve scattering problems [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], [START_REF] Tarantola | Inverse problem theory and methods for model parameter estimation[END_REF] utilizing BA have been previously introduced for various different cases including, among other things, special domain structures such as the cylindrical geometry [START_REF] Leone | Inverse scattering under the distorted born approximation for cylindrical geometries[END_REF], [START_REF] Brancaccio | Information content of born scattered fields: results in the circular cylindrical case[END_REF], [START_REF] Pierri | Dielectric profiles reconstruction via the quadratic approach in 2-d geometry from multifrequency and multifrequency/multiview data[END_REF], [START_REF] Pierri | Second-order iterative approach to inverse scattering: numerical results[END_REF]; advanced inversion of spectral information, e.g., via multiple signal classification (MUSIC) or other methods in the time-reversal of signals [START_REF] Devaney | Time reversal imaging of obscured targets from multistatic data[END_REF], [START_REF] Ciuonzo | On time-reversal imaging by statistical testing[END_REF], [START_REF] Ciuonzo | Noncolocated time-reversal music: Highsnr distribution of null spectrum[END_REF], [START_REF] Ciuonzo | Performance analysis of time-reversal music[END_REF]; and non-linear processes, e.g., distorted Born iterative (DBI) techniques [START_REF] Leone | Inverse scattering under the distorted born approximation for cylindrical geometries[END_REF], [START_REF] Haddadin | Adaptive regularization of a distorted born iterative algorithm for diffraction tomography[END_REF], [START_REF] Lavarello | A study on the reconstruction of moderate contrast targets using the distorted born iterative method[END_REF], [START_REF] Lavarello | Tomographic reconstruction of threedimensional volumes using the distorted born iterative method[END_REF], [START_REF] Hesford | Fast inverse scattering solutions using the distorted born iterative method and the multilevel fast multipole algorithm[END_REF], [START_REF] Chew | Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method[END_REF]. Moreover, regularization techniques such as TV constraints can be used in reconstructing structural distributions based on electromagnetic measurements [START_REF] Abubaker | Total variation as a multiplicative constraint for solving inverse problems[END_REF].

Differing from the aforementioned studies, we approach the tomography of a per se complex-e.g., non-convex and irregular-target applying a multigrid version of the finite element time-domain (FETD) method [START_REF] Li | Time-domain finite element methods for Maxwell's equations in metamaterials[END_REF], [START_REF] Jin | The Finite Element Method in Electromagnetics[END_REF]. Multigrid-FETD enables the accurate modeling of an arbitrary target, maintaining its actual shape with various different finite element (FE) mesh resolutions. It also allows one to speed up the inversion process, which can be performed using a coarse mesh resolution determined by the size of the smallest detectable details. We have previously introduced a multigrid-FETD inversion approach [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF] in which the first-order BA of the wave scattering is used. It has been successfully applied in linearized and iterative TV regularization to reconstruct a synthetic permittivity distribution within a real 3D asteroid shape in [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF], [START_REF] Sorsa | Bistatic full-wave radar tomography detects deep interior voids, cracks and boulders in a rubble-pile asteroid model[END_REF]. Here, building on our previous work, we introduce an iterative and fully non-linear inverse solver in which the unknown permittivity is updated via subsequent linearized approximations akin to the DBI approach and the numerical wave-field can be updated through a higher-order BA. This update requires performing a computationally intensive deconvolution routine. The multigrid-FETD provides the The total wave-field at r. p

The discretized wave-field at -th time point q Scalar and gradient perfectly matched layer matrix A, B (k) Diagonal weight and gradient evaluation matrix, resp.

Q(i)

Restriction matrix with a single non-zero entry

Q (i) i,i = 1 u, g
Electric field and its gradient integrated over time v, w Scalar and vector valued test function, resp. T , T j Finite element mesh and its j-th element, resp.

ϕ i i-th nodal FE basis function χ j
Characteristic function of the j-th element . essential means to tackle the computational cost of this update. Our implementation is available online as an open toolbox for Matlab (The Mathworks Inc.). To enhance the computational performance, Matlab's features for graphics processing unit (GPU) acceleration are applied in both the forward modeling and inversion stage.

The numerical experiments involve a two-dimensional test domain and simulated data. The parameter selection in the numerical experiments was made regarding the radar tomography of an SSSB-a small asteroid in particular-as a potential application. The results suggest that, using the present approach, a candidate solution for the full non-linear inverse problem can be found in a sufficiently short time with improved accuracy when compared to solving a linearized inverse problem. Furthermore, it seems that maximizing the benefit of the potential future bistatic [START_REF] Herique | A Direct Observation of the Asteroids Structure from Deep Interior to Regolith: Two Radars on the AIM Mission[END_REF], [START_REF] Bambach | Discus-the deep interior scanning cubesat mission to a rubble pile near-earth asteroid[END_REF] and multistatic measurements [START_REF] Geffrin | Optimization of a bistatic microwave scattering measurement setup: From high to low scattering targets[END_REF] yielding a sparse set of data might necessitate applying a higher-order BA. This article is organized as follows. Section II describes mathematically the BA, forward simulation and inversion process. Additionally, the underlying wave propagation model is briefly reviewed in Appendix A. Section III includes the numerical results, and Section IV discusses the outcome. Table I lists the essential mathematical symbols of this study.

II. METHODS

A. Sparse full-wave tomography in time-domain

We aim at resolving the relationship between the relative electrical permittivity distribution ε r , and the full timedependence of the wave defined in domain Ω and received at a given set of measurement points. For reconstructing ε r , the ability to model full-wave measurements is essential, since it allows distinguishing signal fluctuations arriving from different parts of the tomography target D, e.g., from the surface or the deep interior. It is necessary to tackle the potential modeling errors arising from the complexity, for example, the effects related to anisotropic structures or scattering effects with high uncertainty. In particular, to optimize the reconstruction quality in a backscattering measurement, it is necessary to filter out the signal peak arising from the wave scattered by the rear wall of D [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF], [START_REF] Sorsa | Bistatic full-wave radar tomography detects deep interior voids, cracks and boulders in a rubble-pile asteroid model[END_REF]. Namely, the amplitude of such a peak can be larger compared to the earlier arriving peaks which are essential for detecting the scatterers in the interior structure. Time-domain full-wave modeling allows this via restricting the investigated signal length as depicted in Fig. 1. Right: The signal (solid grey wave), i.e., the electric field E, corresponding to the measurement configuration on the right has peaks at time (t) points corresponding to the locations (A), (B) and (C). The intense and potentially noisy peak (C) arriving from the rear wall of D needs to be excluded from the final data in order to improve the signal-to-noise ratio of the measurement. The first-order BA takes into account the scattered wave-front which originates at (A) and interacts with the domain (dashed grey wave). The second-order BA takes into account also the two times scattered part of the wave-field, thereby leading to a further correction and an improved approximation for the later arriving peak (dashed black wave) scattered from (B).

B. Green's functions

Assuming that the signal f is transmitted at the point p 1 and d is received at p 2 , the problem of modeling the full wave can be associated with the task of finding the Green's function [START_REF] Duffy | Green's Functions with Applications[END_REF] G p1, p2 predicting the wave received according to the equation d = G p1, p2 * f . The Green's function is a functional of the relative permittivity ε r and a nuisance parameter θ, i.e., G = G[ε r , θ]. That is, the Green's function is not exactly known because of the different uncertainties in the mathematical model. The unknown exact effect of θ on the outcome is related to (numerical) modeling errors, for example, to ambiguities arising from signal attenuation effects which can be related to various factors, e.g., to the conductivity σ and its indirect relationship to ε r [START_REF] Barriot | A two dimensional simulation of the CONSERT experiment (radio tomography of comet wirtanen)[END_REF], propagation losses, reflections, and refractions. For simplicity, it is modeled here with additive random effects by assuming that

G[ε r , θ] = G[ε r ] + E, (1)
where E is an unknown modeling error. In this study, the numerical approximation of the Green's function is first obtained in a forward simulation process after which it is used in the inversion stage. In the latter case it is significant that the dependence of G on ε r is non-linear. Consequently, in order to optimize the modeling accuracy, the approximation for G[ε r ] needs to be updated during an iterative inversion process in which ε r is gradually updated.

C. Higher-order Born approximation of wave scattering

Assume that the Green's function G[ε bg ] of a background permittivity distribution ε bg is given and the task is to approximate G[ε r ], where ε r = ε bg + ρ ε with ρ ε denoting a local perturbation of the relative permittivity at r. That is, a small scattering obstacle at r acts as a point source, whose amplitude is proportional to the local wave-field h. Since the wavefield is altered only at r, Green's function equals to G[ε bg ] elsewhere. The perturbed wave-field is, consequently, of the form (Fig. 2)

d = G[ε bg ] p1, p2 * f + c G[ε bg ] r, p2 * h, (2)
where c is a constant contrast factor, whose numerical dependence on the perturbation 2 is to be determined in the following sections, and

h = G[ε r ] p1, r * f . (3)
The n-th order Born approximation (BA) of h is to assume that h is an n-th degree polynomial with respect to the first term of (2). The first-order BA follows from substituting

G[ε r ] with G[ε bg ], i.e., hBA,1 = G[ε bg ] p1, r * f . (4)
The higher-order (n-th order) approximation can be derived from the following recursive equation:

hBA,n = G[ε bg ] p1, r * f + c G[ε bg ] r, r * hBA,n-1 . (5)
As depicted in Fig. 2, the first-order BA is based on G[ε bg ] and, therefore, it cannot reproduce the non-linear propagation effects in which the path of the altered wave crosses the point r two or more times. Taking such effects into account necessitates a second or higher-order BA which can improve the accuracy of the signal tail, i.e., later arriving peaks, as illustrated in Fig. 1.

D. Evaluation of Green's function via regularized deconvolution

Evaluation of Green's function for any complex-structured ε r necessitates advanced computations which cannot be performed exactly within a feasible time for each scattering point r. Therefore, it is approximated by the following regularized deconvolution process [START_REF] Dimri | Deconvolution and Inverse Theory: Application to Geophysical Problems[END_REF] (Fig. 2):

1) Evaluate the terms h = G[ε r ] p1, r * f and p = G[ε r ] p2, r * f for each scattering point r by propagating a single wave from both p 1 and p 2 .

2) Estimate the Green's function

g = G[ε r ] r, p2 = G[ε r ] p2, r (6)
by using Tikhonov regularized deconvolution with a suitably chosen regularization parameter δ. The reciprocity of the wave propagation ensures that the equation ( 6) holds. The estimate obtained for g can be then applied to estimate the wave d = g * f originating at the scattering point r.

E. Numerical higher-order Born approximation

We assume that an incident wave-field p for time points = 1, 2, . . . , n T has been modeled numerically for a given distribution ε bg and the task is to obtain it for ε r = ε bg + ρ ε , where ρ ε corresponds to a small scattering obstacle. The timeevolution of the wave-field obeys the following so-called leapfrog formulae derived in Appendix A:

q (k) + 1 2 = q (k) -1 2 +∆t a (k) -1 2 p +1 = p +∆tC -1 (f + b + 1 2 ). (7)
Here q

(k) -1 2
is the gradient of p integrated over time and a

(k) -1 2
and b + 1 2 are auxiliary time-evolution vectors (Appendix A). Matrix C is a permittivity-weighted positive definite mass matrix [START_REF] Braess | Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics[END_REF] with entries of the form C i,j = Ω ε r ϕ i ϕ j dΩ, where ϕ i and ϕ j denote FE basis functions (Appendix A). It can be decomposed as

C = C 1 + C 2 , (8)
where C 1 and C 2 correspond to ε bg and ρ ε , respectively. If ρ ε is small enough so that C -1 1 C 2 < 1, the inverse of C can be expanded via the geometric Born series expression for

(I + C -1 1 C 2 ) -1 as given by C -1 = (C 1 + C 2 ) -1 = (I + C -1 1 C 2 ) -1 C -1 1 = (I + C -1 1 C 2 + C -2 1 C 2 2 + • • • )C -1 1 . (9)
The first-degree polynomial approximation of this identity is of the form

C -1 ≈ (I + C -1 1 C 2 )C -1 1 .
When substituted into the leap-frog formulae [START_REF] Sava | 3d radar wavefield tomography of comet interiors[END_REF], this results in the following expression:

q (k) + 1 2 = q (k) -1 2 +∆t a (k) -1 2 p +1 = p +∆tC -1 1 (h + f + b + 1 2 ), (10)
where

h = C 2 C -1 1 (f + b + 1 2
) denotes an auxiliary source vector. We denote

h BA,1 = C 2 C -1 1 (f + b + 1 2 ) (11)
to emphasize that this particular choice for the auxiliary source results in the first-order BA. That is, at each order, the change of the total field

C -1 (f + b + 1 2 ) is approximated with that of the incident field C -1 1 (f + b + 1 2 )
. Formula (10) can be applied recursively n times to obtain the n-th order BA. The auxiliary source corresponding to the result of the recursion is: where P n is a matrix-valued polynomial of the form

h BA,n = C 2 P n (C -1 1 C 2 )C -1 1 (f + b + 1 2 ), (12)
P n (C -1 1 C 2 ) = I + C -1 1 C 2 + • • • + C -n+1 1 C n-1 2 . ( 13
)
The formula ( 12) can be derived inductively as follows. Using the first-order formula two times recursively results in the second order form:

h BA,2 = C 2 C -1 1 (h BA,1 + f + b + 1 2 ) = C 2 C -1 1 (C 2 C -1 1 (f + b + 1 2 ) + f + b + 1 2 ) = C 2 (I + C -1 1 C 2 )C -1 1 (f + b + 1 2 ) = C 2 P 2 (C -1 1 C 2 )C -1 1 (f + b + 1 2 ). (14)
Further, the assumption that the formula

h BA,n = C 2 P n (C -1 1 C 2 )C -1 1 (f + b + 1 2 ) (15)
holds for an arbitrary natural number n, is verified by the induction step:

h BA,n+1 = C 2 C -1 1 (h BA,n + f + b + 1 2 ) = C 2 C -1 1 (C 2 P n (C -1 1 C 2 )C -1 1 (f +b + 1 2 ) +f + b + 1 2 ) = C 2 P n+1 (C -1 1 C 2 )C -1 1 (f + b + 1 2 ). (16)
The n-th order approximation [START_REF] Verma | Elementary Theory of Scattering[END_REF] tends to the exact solution as n → ∞. This can be shown as follows:

C -1 1 (h BA,∞ +f +b + 1 2 ) = (C -1 1 C 2 P ∞ (C -1 1 C 2 ) + I) •C -1 1 (f +b + 1 2 ) = P ∞ (C -1 1 C 2 )C -1 1 (f +b + 1 2 ) = C -1 (f +b + 1 2 ), (17)
where the last identity follows from the fact that

C -1 = P ∞ (C -1 1 C 2 )C -1 1 .

1) Convergence of the higher-order Born approximation:

The condition number κ of an unweighted mass matrix is independent of the mesh size within a quasi-uniform FE mesh, i.e., a mesh in which the ratio between the maximum and minimum element size is uniformly bounded [START_REF] Kamenski | Conditioning of finite element equations with arbitrary anisotropic meshes[END_REF]. That is, for a quasi-uniform mesh, κ = C -1 C = constant, if ε r is constant. It follows that the weighted mass matrices C 1 and C 2 arising from the same FE discretization satisfy the following inequality: showing that the condition C -1 1 C 2 < 1 for the convergence of the Born series ( 9) and, thereby, the higher-order Born approximation, can be satisfied only if the magnitude of the perturbation ρ ε is moderate compared to the background ε bg . The minimum min x ε bg can be considered in the local neighborhood of the perturbation, since any vector multiplied by C -1 1 C 2 is negligible far from the perturbation. Namely, first multiplication by C 2 restricts any vector to the perturbed neighborhood which is, then, spread by some amount, when multiplied by C -1 1 , since the solution of the linear system defined by a mass matrix decays, when moving away from the source.

C -1 1 C 2 ≤ C -1 1 C 2 ≤ constant × max x |ρ ε | min x ε bg , (18)
As the BA is here found through the Tikhonov regularized deconvolution process (Section II-D) and the source h BA,n of BA is linear respect to C 2 , selecting the magnitude of ρ ε can be associated with choosing an appropriate regularization parameter value δ. Namely, an update of the form ρ ε → γρ ε with scaling γ > 0 corresponds to h BA,n → γh BA,n and alternatively to δ → γ -1 δ, if the estimate p obtained for the wave is updated as p → γp . That is, the same effect which follows by decreasing the perturbation can be generated via increasing the regularization parameter. Therefore, in the numerical implementation, it suffices to assume that ρ ε = 1 in the numerical evaluation of the BA, and to select an appropriate regularization parameter with respect to that.

2) Multigrid formulation of the permittivity: We define the permittivity distribution with respect to a coarse and nested dsimplex mesh T (Fig. 3), i.e., ε r = M j=1 c j χ j and assume that the permittivity distribution is piecewise constant as given by ε r = m j=1 c j χ j . The multigrid formulation, i.e., the use of a coarse inversion mesh, is applied to set the resolution of the unknown permittivity distribution on a suitable level. Namely, the maximal theoretical reconstruction accuracy obtainable with full data is generally lower than what is needed for propagating the wave.

In order to obtain a feasible performance with the deconvolution approach presented in II-D, it is assumed that ε r is Forward simulation

Source points

Model the full-wave propagation using each of the following points as the source point p1:

• transmission positions • the nodes of the mesh T .

Wave propagation

Perform leap-frog iteration in T .

Data points

Store the full-wave data using each of the following points as the data point p2:

• receiver positions • the nodes of the mesh T .

Inversion process

Data points

Use the full-wave measurements obtained at the receiver positions as the data for the inversion process.

Total variation regularization

Find a new estimate for εr by performing the iterative TV regularized inversion process for one or more multigrid meshes nested with respect to T can be used as explained in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF]. For each mesh, form a Jacobian matrix with respect to εr as described in Section II-F.

Update wave-field

After the TV regularized iteration, correct the wave-field to correspond εr at an update restricted to T j , we define the following element-wise auxiliary source vector which is used in the deconvolution process:

h (i,j) = C(j) 2 Q (i) C -1 1 (f + b + 1 2 ). (19)
Here, the matrix

Q (i) ∈ R n×n has one nonzero entry Q (i)
i,i = 1 and the vector h (i,j) differs from zero only if the i-th node belongs to T j . The solution of the system (10) with respect to h (i,j) can be obtained as

p ≈ ri∈T j , i≤N d (i,j) = d+1 k=1 d (i k ,j) , (20)
where the number of terms is d + 1, that is, the number of nodes in T j , and d (i,j) is the regularized deconvolution-based solution of the following auxiliary system:

r (i,j,k) + 1 2 = r (i,j,k) -1 2 + ∆ta (i,j,k) -1 2 , d (i,j) +1 = d (i,j) + ∆tC -1 (h (i,j) + b (i,j) + 1 2 ). (21)
This system can be derived from (10) simply by substituting h (i,j) as the source. The definition for the auxiliary vectors

a (i,j,k) -1 2 and b (i,j) + 1 2
follow directly by substituting p and d

(k) -1 2
with d (i,j) and r

(i,j,k) -1 2
, respectively, in Equation ( 28) and (29) (Appendix A). Note that only d

(i,j)
+1 out of the auxiliary variables r (i,j,k) + 1 2 and d (i,j) +1 needs to be evaluated. This is approximated via the regularized deconvolution process. Hence, with the multigrid approach, one can perform the deconvolution process with respect to the coarse system ( 21) by applying the wave data obtained with the dense one. A graphical description of the multigrid solver is shown in the Fig. 4.

F. Jacobian matrix

As shown in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF], the wave equation can be linearized with respect to the permittivity distribution as follows:

∂q (k) + 1 2 ∂c j = ∂q (k) -1 2 ∂c j +∆t ∂a (k) -1 2 ∂c j ∂p +1 ∂c j = ∂p ∂c j +∆tC -1 (h diff + ∂b + 1 2 ∂c j ), (22)
where

h diff = ∂C ∂c j C -1 (f + b + 1 2 ) ( 23
)
is an auxiliary source function implied by the identities (∂C -1 /∂c j )(f +b + 1 2 ) and ∂C -1 /∂c j = C -1 (∂C/∂c j )C -1 the latter one of which can be obtained via differentiating both sides of CC -1 = I which gives ∂C -1 /∂c j = C -1 ∂C/∂c j C -1 . Based on the definition of the permittivity distribution and the matrix C, it holds that (∂C/∂c j ) k, = Tj ϕ k ϕ dΩ. That is, the entry (∂C/∂c j ) k, is non-zero only if for which ϕ k and ϕ are supported on T j .

When interpreted in terms of Section II-E, h diff can be interpreted as a special case of the first-order BA [START_REF] Carry | Density of asteroids[END_REF]. Morever, the multigrid formulation presented in Section II-E2 is valid, since the update (∂C/∂c j ) differs from zero only in the element T j ∈ T .

G. Non-linear inversion process

The data was inverted using the following fully non-linear inversion process (Figure 4):

1) Choose a constant initial guess x (0) for the unknown permittivity distribution, a desired number of iterations N , and set = 1. 2) Find an estimate x ( ) for the permittivity distribution using the iterative total variation (TV) regularization algorithm described in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF] with one or more multigrid meshes nested with respect to T . For each mesh, evaluate a Jacobian matrix with x ( -1) as the linearization point. 3) If < N , update the wave-field via BA with x ( ) as the reference permittivity distribution, set → + 1 and repeat the second step. Perform the update sequentially element-by-element with the element-wise strategy described in Section II-E2. 4) Associate the final estimate with the reconstruction of the permittivity distribution.

This algorithm is a fully non-linear generalization of the reconstruction method presented in the study [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF] which is restricted to the linearized case N = 1. Note that the iterative TV regularization process of the second step can be non-linear also when N = 1. It also allows using an arbitrary number of multigrid hierarchy levels in a straigthforward manner as explained in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF]. Here, for simplicity without losing the generality, the investigation is restricted to two hierarchy levels T and T .

H. Computational complexity 1) Forward approach: The computational complexity of the forward solver follows from the signal's shortest wavelength which determines the mesh size, h. The number of spatial degrees of freedom n in the FE mesh is proportional to h -d , in which d is the dimension of the mesh. That is, the number of non-zeros in the system matrices grow as O(n). Hence, the complexity of each time step is O(n), including the multiplication by C -1 . Namely, the complexity of the solution of the linear system determined by C by preconditioned conjugate gradient method [START_REF] Golub | Matrix Computations[END_REF], applied in inverting the mass matrix C, is O(n √ κ) and the condition number κ of C is independent of n in a quasi-uniform mesh (Section II-E1). The total complexity of the forward solution is, thus, O(n 1+1/d ) since the number of time steps is depends on the degrees of freedom along a single dimension, i.e. it is of the complexity of O(n 1/d ). [START_REF] Annan | Electromagnetic principles of ground penetrating radar[END_REF].

2) Inverse approach: In the present multigrid scheme, the resolution of the mesh T is not bound by that of the forward simulation but rather by that of the desired imaging resolution. Following from the radar range resolution, the maximal imaging resolution is determined by the signal bandwidth [START_REF] Daniels | Ground Penetrating Radar (2nd Edition)[END_REF]. Evaluation of the first-order BA for a given signal transmission requires performing the regularized deconvolution process (Section II-D), i.e., accounting the effect of the signal propagating through the target object D for each of the M elements in T at each of the K measurement points, resulting in the complexity of O(KM ). Since the number of elements is comparable to that of the nodes, this can be expressed also as O(KN ), where N refers to the number of nodes. When evaluating a higher-order transform, the regularized deconvolution is evaluated for the combined set of nodes and measurement points. That is, the complexity is O(KN 2 ) with respect to the spatial degrees of freedom. The computational cost is also directly proportional to the order of the transform. The evaluation of a BA essentially determines the computational complexity of the inversion process (Section II-G).

I. Numerical implementation with the GPU-accelerated Tomographic Radar Reconstruction (GPU-ToRRe) toolbox

Numerical forward and inverse methods were implemented for the MATLAB (The MathWorks Inc.) programming environment which natively allows using GPU accelerated algorithms as a part of the code. The GPU-accelerated Tomographic Radar Reconstruction (GPU-ToRRe) toolbox 1 de-1 https://github.com/sampsapursiainen/GPU-Torre veloped in this study is available online in GitHub. The documentation can be found therein. The numerical results of this study have been computed using this toolbox.

J. Domain

Numerical experiments were performed in a twodimensional origin-centric square Ω including the tomography target D (Figure 5) in the center part and a perfectly matched layer (PML) [START_REF] Schneider | Understanding the FDTD Method[END_REF] near the boundaries. The PML was embedded in the model as shown in Appendix to simulate an open-field wave propagation, i.e., to dampen echoes from the boundaries back to the center. For the generality of the results, we used a unitless set of parameters t, x, ε r , σ and c = ε -1/2 r (velocity) which can be scaled to SI-units via the expressions (µ 0 ε 0 ) 1/2 st, s x, ε 0 ε r , (ε 0 /µ 0 ) 1/2 s -1 σ, and (ε 0 µ 0 ) -1/2 c (c = 1 for ε r = 1), respectively, with s denoting a spatial scaling factor (meters), ε 0 = 8.85 • 10 -12 F/m the electric permittivity of vacuum, and µ 0 = 4π • 10 -7 H/m the magnetic permeability which is assumed to be constant in Ω. By tuning s, the system can be scaled to match with different scales and applications. That is, the actual target and its analog scale model can be modeled through a single system but two different values of s. Table II shows the scaling for s = 1, s = 500 m and s = 0.5 m. The first one of these is the unitless representation, the second one corresponds to a radio frequency measurement performed for a SSSB, and the last one to a laboratory-scale microwave experiment. Conductivity was considered as an unknown latent nuisance parameter predicted by the equation σ = 5ε r , as it difficult to be inverted as shown, e.g., in [START_REF] Barriot | A two dimensional simulation of the CONSERT experiment (radio tomography of comet wirtanen)[END_REF].

In the unitless coordinates, the diameter of D was approximately 0.28, or in terms of wavelength, 8.7λ. The relative permittivity distribution of D included a surface layer with the thickness around 0.02 (0.6λ) and three inclusions (voids) with maximum diameter of 0.01-0.09 (0.3-2.8λ) (Figure 5. The values of ε r for the interior part, the surface layer and the voids were chosen to be 4, 3 and 1, respectively. Outside D, ε r was set to be one, i.e., that of air or vacuum. Inside D, conductivity causing signal energy loss was assumed to be a nuisance parameter of the form σ = 5ε r , and vanish elsewhere, i.e., σ = 0.

1) Numerical discretization: To avoid the inverse crime, i.e., the overly good data fit in the inversion process, the exact and background wave data were simulated using two different triangular (shape-regular and unstructured) FE meshes consisting of 93 475 and 40 715 nodes together with 186 544 and 81 040 triangles. The time increment ∆t in the leap-frog iteration was set to be 6.25E-5 and 2.5E-4, respectively.

Each triangular mesh applied in the wave propagation process (forward simulation) was obtained by refining a coarse one uniformly two times. The permittivity distribution was reconstructed for an original coarse mesh (996 triangular elements, 552 nodes) which was nested with respect to the one corresponding to the background data.

2) Signal pulse and measurement: The Blackman-Harris window [START_REF] Irving | Numerical modeling of ground-penetrating radar in 2-D using MATLAB[END_REF], [START_REF] Harris | On the use of windows for harmonic analysis with the discrete fourier transform[END_REF], [START_REF] Nuttall | Some windows with very good sidelobe behavior[END_REF] with the duration T 0 was used as the shape of the source function f (t) for t ≤ T 0 (0.67 ns), and f (t) = 0, otherwise. The duration of the pulse was chosen to be T 0 = 0.1 and the duration of the measurement T = 1.1. The time interval between each data sampling point was set to be 0.005 corresponding to a 2.5 oversampling rate relative to the Nyquist criterion. The signal properties and classification can be found in Table III.

The signal was transmitted and received at 0.32 diameter circular path centered at the origin and enclosing D. We investigated the following four different spatial measurement configurations depicted in Figure 6: 1) The monostatic configuration is constituted by a single spacecraft.

2) The bistatic I configuration includes two spacecraft with a constant 22.5 degree angle in between them. 3) In the bistatic II configuration, the separating angle is 90 degrees. 4) In the multistatic configuration, a 90 degrees angle is covered by altogether five equally spaced spacecraft, each two separated by a 22.5 angle. In each one of these, the red one both transmits and receives the signal while the other spacecraft are used as additional receivers. The dataset included a total number of 16 transmission points equally distributed around the target D. This sparse distribution not satisfying the Nyquist criterion is used in order to take into account the in-situ restrictions of a radar measurement performed during a planetary space mission. Achieving a full spatial measurement coverage would necessitate using a point density which would oversample the Nyquist criterion by factor two in the vicinity of D. The corresponding number of points can be obtained dividing the wavelength of the highest frequency signal component in vacuum (here 0.03) by four times the circumference of the circumcircle containing D (here 0.9), resulting here to about 120 equally spaced points [START_REF] Bucci | Electromagnetic inverse scattering: Retrievable information and measurement strategies[END_REF].

3) Noise: Zero-mean Gaussian white noise with a fixed standard deviation (STD) was added to the measurements. To investigate the effect of the noise on the signal, the peak-topeak signal-to-noise ratio (PPSNR) was evaluated for each measurement configuration. It was defined as the decibel value of the relative noise peak level with respect to the amplitude of the initial data vector, i.e., the difference between the exact signal and the background data. The noise peak level was defined to be the 95 % quantile of the Gaussian noise distribution. When selecting the noise STD, the targeted level of PPSNR was between 10 and 20 dB (with 0 dB referring to the level of the signal peak) which is known to allow finding a tomographic reconstruction and also seems a potential range for measurements performed for an SSSB based on the CONSERT data [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF], [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF]. The motivation to use Gaussian noise is the significance of the modeling errors, whose net effect approach a Gaussian random variable under the additive uncertainty model, assuming that the errors are, additionally, independent and identically distributed.

4) Inverse estimates: The final reconstruction of the permittivity distribution was obtained via three steps of the nonlinear inversion algorithm in which, on each step, a regularized permittivity estimate was found through a single step of the iterative TV regularization scheme presented in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF]. In this algorithm, the regularizing function is a sum of TV and L2-norm penalty term. The TV term corresponds to the Euclidean norm of the permittivity gradient integrated over D. It regularizes the jumps between adjacent elements, while the L2-norm penalizes the total magnitude of the distribution. These terms were scaled by the parameter values α = 2E-1 and β = 1E-3, respectively. The predominating parameter α affects the reconstruction quality and the role of β is mainly to ensure the numerical stability of the inversion process, i.e., the boundedness of the reconstruction. The value of the Tikhonov regularization parameter for the deconvolution process was set to be δ = 1E-4. Each parameter value was selected based on preliminary tests. The reconstructions were produced and analyzed separately for the first, second and third-order Born approximation.

5) Inverse error measures: The quality of the inverse estimates was measured by evaluating the structural similarity (SSIM) [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF] between the exact and the reconstructed permittivity distribution denoted here by ε r and ε * r . To analyze the value accuracy and localization of the permittivity details, we evaluated also the mean squared error (MSE) and the relative overlap error (ROE). The latter one of these is defined as the relative error in overlap between ε r and ε * r for the surface layer and voids. We define ROE as the percentage

ROE i = 100 1 - Area(A) ∩ Area(S i ) Area(S i ) for i = 1, 2, (24)
where S 1 denotes the voids and S 2 the surface layer, A = (S 1 ∪ S 2 ) ∩ R is the overlap between the target set S 1 ∪ S 2 and the set R in which a given reconstruction is smaller than a limit such that Area(R) = Area(S 1 ∪ S 2 ).

III. RESULTS

The results of the numerical experiments can be found in Tables IV and V and Figure 7. Both the signal configuration and the order of the Born approximation were observed to have a significant effect on the reconstruction quality obtained in the non-linear inversion process.

A. Peak-to-peak signal-to-noise ratio

The Bistatic I and multistatic signal configuration, i.e., those including a measurement point within a 22.5-degree angle from the transmission location, produced a higher measurement PPSNR than the monostatic one. The PPSNR obtained with Bistatic II, in which the additional measurement point is at 90 degrees with respect to the transmitter, was equal to that of the monostatic case.

B. SSIM

A visual inspection of the reconstructions suggests the two and multi-point approaches provide an increased sharpness and distinguishability of the permittivity details. This observation is supported by the SSIM results, showing that the structural similarity between the reconstruction and the exact permittivity distribution increases along with the measurement point count and also with the order of the BA. Based on the SSIM, the Bistatic configurations I and II were found to improve the reconstruction quality by a somewhat similar marginal as compared to the monostatic case, while the superior results were obtained with the multistatic configuration.

C. MSE and ROE

The first-order BA was found to be advantageous with respect to the overlap (ROE), especially, in the cases of the Bistatic I and multistatic configuration, which nevertheless had a lower value-accuracy with regard to MSE as compared to the other two configurations. Increasing the order of the BA led to an enhanced global value-accuracy for each configuration, and based on a visual examination, it also led to a sharper contrast between the image details and the background. The advantage of using a higher-order BA was obvious, especially, with respect to the value of the surface layer, which generally is overly low in the first-order case. Notably, with the second and third-order BA, the Bistatic configuration I performed slightly better compared to II regarding both MSE and ROE.

D. Computing time

The performance of the numerical solver was evaluated using an NVIDIA Quadro P6000 GPU; using it computing the Jacobian took about 45 s, updating the full-wave data 75, 150 and 225 s with the first-, second-, and third-order BA, resulting in a total reconstruction process duration of about 285, 435 and 585 s, respectively. Propagating the wave for a single transmission point took about 45 and 135 s for the background and exact permittivity distribution, respectively. The wave was propagated for each node included in the regularized deconvolution process. The total number of transmission points was 16 outside the target D and 552 inside, corresponding to each node of the coarse inversion mesh. Performing these processes in a standard CPU was observed to generally take 20-100 times the time required by a GPU.

IV. DISCUSSION

In this study, we have shown that the linearized full-wave inversion approach for the tomographic reconstruction of small solar system bodies (SSSBs) presented in [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF] can be derived from the first-order Born approximation (BA) by associating it with the Jacobian matrix of the numerical wave-field. Here, this approach was generalized as a non-linear iteration in which the wave-field-i.e., the forward solution-can be updated via a BA of an arbitrary order. The structural similarity (SSIM) criterion shows an improvement in the reconstructions with the increase in the BA order. We have also prepared a graphics processing unit-accelerated toolbox, GPU-ToRRe for Matlab (MathWorks Inc.), to achieve a sufficiently short computation time for experimentation and further method development.

The structural similarity (SSIM) and overlap (ROE) results suggest that the quality of the reconstruction can be increased via bistatic or multistatic measurement configurations, i.e., dual-or multi-point measurement schemes. Increasing the order of the BA was found to improve the SSIM further, and its role with regard to the value-accuracy (MSE) of the reconstructions was found to be crucial. In the Bistatic I configuration, an angular separation of 22.5 degrees between the transmitter and the additional receiver was observed to result in an enhanced measurement PSNR level compared to the monostatic case, which was recently proposed in [START_REF] Sava | 3d radar wavefield tomography of comet interiors[END_REF], whereas the larger 90-degree angle of the Bistatic II configuration did not. It seems that, whereas the Bistatic I reconstruction is superior to the Bistatic II on the surface part, there is not such a difference in the deep interior, which we interpret to follow from the longer and more "multi-way" signal propagation in comparison to the background prediction in the latter case.

The findings of this study suggest that, in addition to the noise-robustness [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF], [START_REF] Sorsa | Bistatic full-wave radar tomography detects deep interior voids, cracks and boulders in a rubble-pile asteroid model[END_REF], the reconstruction quality can be improved via a bistatic measurement with regard to both structural similarity and value accuracy, especially if a higherorder BA is used. A bistatic configuration had also been chosen for CONSERT to provide tomographic travel-time data [START_REF] Hérique | Direct observations of asteroid interior and regolith structure: science measurement requirements[END_REF], [START_REF] Kofman | Properties of the 67p/churyumov-gerasimenko interior revealed by consert radar[END_REF].

Regarding the practical aspects of a space mission design, placing the transmitter in the mothership might be convenient to guarantee its power supply. Small spacecrafts could carry the additional receivers for the chosen measurement configuration [START_REF] Herique | A Direct Observation of the Asteroids Structure from Deep Interior to Regolith: Two Radars on the AIM Mission[END_REF], [START_REF] Bambach | Discus-the deep interior scanning cubesat mission to a rubble pile near-earth asteroid[END_REF]. Therefore, in addition to the monostatic and bistatic configurations, the multistatic one can also be considered as a potential scenario for a future space mission. In the light of the present numerical results, such a configuration should improve the reconstruction quality as shown by the SSIM, global MSE and ROE criteria.

Referring to the convergence properties, the numerical stability of the BA with the present Tikhonov regularization approach seems appropriate. It converges rapidly as a function of the approximation order, matching with our justification given in Section II-E1 for a well-chosen regularization parameter. Based on our inversion results, it seems that the thirdorder approximation might be applicable for wave propagation investigations, as the effect of increasing the order of the BA is already minor when moving from the second to the third order. The practical limit for increasing the order is set by the noise, the effect of which becomes visible gradually as the order of the BA increases. Here the third-order case was selected as the limit, due to the observation of minor noise effects. Even though the highest SSIM was obtained for the thirdorder reconstructions, they can also be argued to be affected by the noise based on a comparison of MSE values. We deem them to be due to the inversion process, which applies the BA multiple (three) times to correct the wave-field-i.e., the Green's function. The sensitivity of such updates to noise is a well-known property of the DBI techniques [START_REF] Leone | Inverse scattering under the distorted born approximation for cylindrical geometries[END_REF], [START_REF] Haddadin | Adaptive regularization of a distorted born iterative algorithm for diffraction tomography[END_REF], [START_REF] Lavarello | A study on the reconstruction of moderate contrast targets using the distorted born iterative method[END_REF], [START_REF] Lavarello | Tomographic reconstruction of threedimensional volumes using the distorted born iterative method[END_REF], [START_REF] Hesford | Fast inverse scattering solutions using the distorted born iterative method and the multilevel fast multipole algorithm[END_REF], [START_REF] Chew | Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method[END_REF]. Our method is similar to the DBI methods, but it performs the computations in the time domain and allows for complex-shaped target domains. The noise effects could be resolved by strengthening the regularization which, however, would diminish the effect of the higher-order updates, since they would alter and decrease the terms of the Born series. In addition to the noise, another obvious factor affecting the performance of the inversion algorithm is the contrast of the permittivity details-i.e., as shown in Section II-E1, the ratio between the magnitude of the permittivity perturbation and the background permittivity distribution is the key factor determining the level of the Tikhonov regularization necessary for the convergence of the BA. Following from this, in can be challenging to reconstruct high-contrast details [START_REF] Haddadin | Adaptive regularization of a distorted born iterative algorithm for diffraction tomography[END_REF], [START_REF] Lavarello | A study on the reconstruction of moderate contrast targets using the distorted born iterative method[END_REF].

The Born iterative method with a steady Green's function has been shown to be advantageous with respect to the noise robustness of the inversion process compared to a DBI updating method to data in the frequency domain [START_REF] Chew | Reconstruction of two-dimensional permittivity distribution using the distorted born iterative method[END_REF]. This observation is also in agreement with our present timedomain modeling results, suggesting that the update routine necessitates finding a balance between accuracy and noise suppression. The analysis of the noise tolerance of our method with the first-order BA without updating the Green's function (corresponding to the Born iterative method) can be found in our previous work [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF], [START_REF] Takala | Far-field inversion for the deep interior scanning cubesat[END_REF] where a noise level above 8 to 10 dB was found to be necessary for obtaining an appropriate reconstruction quality. This compares well with the findings for microwave tomography of the breast [START_REF] Zeng | Effects of noise on tomographic breast imaging[END_REF], for example. For its superior magnitude, the first-order BA determines the inver-First-order Born approximation sion performance in high-noise cases, where the minor higherorder corrections have a less significant effect in cases where the level of the Tikhonov regularization is sufficient. The firstorder BA is also less demanding regarding computational complexity. Therefore, it might be preferable for dense inversion meshes. Optimizing the performance of the multigrid solver with respect to noise and regularization is a complex topic, as it depends on altogether three regularization parameters: the TV, L2-norm and Tikhonov regularization parameters. It is obvious that the optimal conditions, including the most suitable order of the BA, depend on the application-related parameters, such as the permittivity distribution, wavelength, targeted imaging resolution, and the distribution of the scattering obstacles.

R ∩ S1 R ∩ S2 Monostatic R ∩ S1 R ∩ S2 Bistatic I R ∩ S1 R ∩ S2 Bistatic II R ∩ S1 R ∩ S2 Multistatic Second-order Born approximation R ∩ S1 R ∩ S2 Monostatic R ∩ S1 R ∩ S2 Bistatic I R ∩ S1 R ∩ S2 Bistatic II R ∩ S1 R ∩ S2 Multistatic Third-order Born approximation R ∩ S1 R ∩ S2 Monostatic R ∩ S1 R ∩ S2 Bistatic I R ∩ S1 R ∩ S2 Bistatic II R ∩ S1 R ∩ S2 Multistatic
Regarding the goal of reconstructing the interior structure of an asteroid, it will be important to implement and evaluate the present higher-order BA for a 3D computation geometry. This could be done by extending the linearized multigrid solver introduced and used in our recent studies [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF], [START_REF] Sorsa | Bistatic full-wave radar tomography detects deep interior voids, cracks and boulders in a rubble-pile asteroid model[END_REF]. Applying the multigrid approach is crucial in discretizing the unknown permittivity distribution, as the sparsity of the inversion grid enables the fast performance of the wave-field updating routine. The results of this study together with the adaptability of the multigrid mesh suggest that our approach will be sufficiently fast also in a three-dimensional case where the coarse mesh might consist of a few thousand elements. A three-dimensional inverse solver could be evaluated via microwave-range measurements on an analogue target. Recently, a multistatic inversion approach has been validated for frequency-domain data in such a way [START_REF] Eyraud | Imaging the interior of a comet from bistatic microwave measurements: case of a scale comet model[END_REF]. Thus, a natural next step would be to perform a validation for the time domain. Other interesting research directions would be to investigate carrier wave effects in the full-wave inversion of tomographic radar data, and sampling techniques to improve the noiserobustness of the reconstructions [START_REF] Eyraud | Microwave imaging from experimental data within a bayesian framework with realistic random noise[END_REF].

APPENDIX

A. Wave propagation model

Following [START_REF] Takala | Multigrid-based inversion for volumetric radar imaging with asteroid interior reconstruction as a potential application[END_REF], the wave equation for a transverse electric field E = u e z in the spatio-temporal set [0, T ] × Ω can be formulated as a first-order system ε r ∂u ∂t + σu -∇ • g = f and ∂ g ∂t

-∇u = 0, (25)
where ε r and σ denote the relative permittivity and conductivity distribution, respectively, g = t 0 ∇u(τ, x) dτ , and the boundary conditions are set by g| t=0 = ∇u 0 and u| t=0 = u 1 . The right-hand side is the signal source, i.e., the antenna current density given by f (t, x) = δ p ( x) f (t) transmitted at the point p. Here δ p ( x) is Dirac's delta function with respect to p. Integrating

Here, it is assumed that the domain and the parameters are regular enough, so that the weak form has a unique solution u : [0, T ] → H 1 (Ω) [START_REF] Evans | Partial Differential Equations[END_REF].

To obtain a numerical solution, we assume that the electric and gradient field are finite sums of the form u = n j=1 p j ϕ j and g = d k=1 g (k) e k with g (k) = m i=1 q (k) i χ i , respectively. Here, ϕ 1 , ϕ 2 , . . . , ϕ n are real-valued linear (nodal) basis functions and χ 1 , χ 2 , . . . , χ n are piecewise (element-wise) constant element-indicator functions. Defining test functions v : [0, T ] → V ⊂ H 1 (Ω) and w : [0, T ] → W ⊂ L 2 (Ω) with V = span{ϕ 1 , ϕ 2 , . . . , ϕ n } and W = span{χ 1 , χ 2 , . . . , χ m } the weak form can be written in the Ritz-Galerkin discretized form [START_REF] Braess | Finite Elements: Theory, Fast Solvers, and Applications in Solid Mechanics[END_REF], that is, ∂ ∂t Aq (k) -B (k) p + T (k) q (k) = 0,

∂ ∂t Cp + Rp + Sp + d k=1 B (k) T q (k) = f , (27)
with p = (p 1 , p 2 , . . . , p n ) and q (k) = (q

(k) 1 , q (k) 2 , . . . , q (k)
m ), and f ∈ R n with f i = Ω f ϕ i dΩ denoting the coordinate vectors for u, g and the source function, respectively.

The matrices of the system are given by A ∈ R m×m , B ∈ R m×n , C ∈ R n×n , S ∈ R n×n , T ∈ R m×m . Matrices A and T (k) = ζ (k) A are diagonal with non-zero entries determined by A i,i = Ti dΩ. The matrix B (k) is a projection matrix of the form B (k) i,j = Ti e k • ∇ϕ j dΩ, and C, R and S are mass matrices weighted by ε r , σ and ξ, respectively, as given by C i,j = Ω ε r ϕ i ϕ j dΩ, R i,j = Ω σ ϕ i ϕ j dΩ and S i,j = Ω ξ ϕ i ϕ j dΩ. The matrices S and T (k) correspond to a split-field perfectly matched layer (PML), i.e., the set To discretize the time interval [0, T ], we utilize ∆t spaced regular grid of n T time points and the standard difference approximations for the time derivative. These substituted into [START_REF] Hesford | Fast inverse scattering solutions using the distorted born iterative method and the multilevel fast multipole algorithm[END_REF] lead to the leap-frog formulae [START_REF] Sava | 3d radar wavefield tomography of comet interiors[END_REF] [START_REF] Schneider | Understanding the FDTD Method[END_REF], [START_REF] Bossavit | Yee-like schemes on a tetrahedral mesh, with diagonal lumping[END_REF], [START_REF] Yee | Numerical solution of initial boundary value problems involving maxwell's equations in isotropic media[END_REF], where the auxiliary time-evolution vectors a

{ x ∈ Ω | 1 ≤ max k |x k | ≤ 2 }
a (k) -1 2 = A -1 B (k) p -A -1 T (k) q (k) -1 2 , ( 28
) b + 1 2 = -Rp -Sp - d k=1 B (k) T q (k) + 1 2 . (29)

  component of integrated gradient at -th time point a of p at -th time point h Discretized total wave-field at r at -th time point hBA,n n-order Born approximation of h h (i,j) hBA,1 w.r.t. j-th element and i-th node in a coarse mesh d (i,j) p w.r.t. j-th element and i-th node in a coarse mesh h diff hBA,1 obtained as a differential C, C 1 , C 2 Mass matrices weighted w.r.t. εr, ε bg and ρε, resp. R Mass matrix weighted w.r.t. σ S, T (k)

Fig. 1 .

 1 Fig. 1. Left: A schematic picture of a tomography target D, e.g. an asteroid, in which the electromagnetic wave transmitted by an antenna (e.g. an orbiting spacecraft) is scattered at different locations denoted by (A), (B) and (C).Right: The signal (solid grey wave), i.e., the electric field E, corresponding to the measurement configuration on the right has peaks at time (t) points corresponding to the locations (A), (B) and (C). The intense and potentially noisy peak (C) arriving from the rear wall of D needs to be excluded from the final data in order to improve the signal-to-noise ratio of the measurement. The first-order BA takes into account the scattered wave-front which originates at (A) and interacts with the domain (dashed grey wave). The second-order BA takes into account also the two times scattered part of the wave-field, thereby leading to a further correction and an improved approximation for the later arriving peak (dashed black wave) scattered from (B).

Fig. 2 .

 2 Fig. 2. (a) A schematic picture depicting the principle of the Tikhonov-regularized deconvolution process which is applied in the Born approximation (BA) of the wave scattering and in forming the Jacobian matrix (as first-order BA). (b) The first-order BA takes into account the scattering wavefronts (black dashed arrows) originating from a single point and its interaction r w.r.t. the existing details in the computation geometry, including the internal part and the surface of the target D. (c)The second-order BA, takes into account, additionally, the wavefronts which have scattered two-times at r (solid gray arrows). Generally, the order of BA determines the maximal number of times the wave can be scattered at r in the model.

Fig. 3 .

 3 Fig. 3. A schematic illustration of the nested mesh structure. The fine mesh T (left) is used in the forward simulation and the coarse mesh, T (right) in the inversion.

•

  the receiver positions • nodes of T for each signal source point used in the forward simulation.

Fig. 4 .

 4 Fig. 4. Graphical description of the forward and inversion stages of the present multigrid solver. The unknown permittivity εr is updated sequentially via a linearized approximation akin to the DBI method [28]. The left column visualizes the point sets involved and the right one includes the descriptions.

2

 2

Fig. 5 .

 5 Fig. 5. The exact tomography target D used in the numerical experiments. The white, gray and black area have the relative permittivity value of 4 (solid layer), 3 (porous layer) and 1 (vacuum), respectively.

Fig. 6 .

 6 Fig. 6. The signal configurations of the numerical experiments. In each one of these, the red one both transmits and receives the signal while the other spacecraft are used as additional receivers. The dataset included a total number of 16 transmission points equally distributed around the target D. 1) The Monostatic configuration is constituted by a single spacecraft. 2) The Bistatic I configuration includes two spacecraft with a constant 22.5 degree angle in between them. 3) In the Bistatic II configuration, the separating angle is 90 degrees. 4) In the Multistatic configuration, a 90 degrees angle is covered by altogether five equally spaced spacecraft, each two separated by a 22.5 angle.

Fig. 7 .

 7 Fig. 7. The reconstructions of the two-dimensional permittivity distribution together with the overlap sets R ∩ S 1 and R ∩ S 2 . The color scale is from 1 (black) to 5 (white).

  [START_REF] Lavarello | A study on the reconstruction of moderate contrast targets using the distorted born iterative method[END_REF] multiplied by v : [0, T ] → H 1 (Ω) and w : [0, T ] → L 2 (Ω) and applying the rule of partial integration, one can obtain the following weak form:∂ ∂t Ω g • w dΩ-Ω w•∇u dΩ = 0, ∂ ∂t Ω ε r uv dΩ+ Ω σ uv dΩ+ Ω g•∇v dΩ = f , if x= p, 0, else.

  which eliminates reflections from the boundary ∂Ω back to the inner part of Ω[START_REF] Schneider | Understanding the FDTD Method[END_REF]. For the PML parameters,ξ( x) = ς, if 1 ≤ max k |x k | ≤ 2 , and ζ (k) ( x) = ς, if 1 ≤ |x k | ≤ 2 , and ξ( x) = ζ (k) ( x) = 0, otherwise.

TABLE I LIST

 I OF ESSENTIAL MATHEMATICAL SYMBOLS. ] p 1 , p 2 Green's function between p 1 and p 2 w.r.t. εr f Signal transmitted at p 1 d Signal received at p 2 . h

	Symbol	Description
	εr	Relative electric permittivity
	ε bg	Background permittivity
	ρε	Local permittivity perturbation
	σ	Conductivity distribution
	G[ r

TABLE II THE

 II PROPERTIES OF THE TEST DOMAIN Ω CORRESPONDING TO THE FOLLOWING VALUES OF THE SCALING PARAMETER: s = 1, s = 500 M AND s = 0.5 M. THE FIRST SI-UNIT VALUE MATCHES WITH A RADIO FREQUENCY MEASUREMENT PERFORMED FOR A SMALL SOLAR SYSTEM BODY[START_REF] Kofman | Radar techniques to study subsurfaces and interiors of the solar system objects[END_REF], AND THE LAST ONE TO A MICROWAVE-RANGE LABORATORY MEASUREMENT (TABLEIII).

	Scaling s	D	Surface layer	Voids
	1	Diameter	0.28		0.02 0.01-0.09 m
			58λ		4λ	2-18λ
		Value of εr	4		3	1
		Value of σ	20		15	5
	500 m	Diameter	140 m		10 m	5 -45 m
		Value of εr	4		3	1
		Value of σ	11.0	8.0E-5 S/m	2.8E-5 S/m
	0.5 m	Diameter	14 cm		1 cm	0.5 -4.5 cm
		εr	4		3	1
		σ	0.11 S/m	0.80 S/m	0.028 S/m
			TABLE III
	SIGNAL PROPERTIES AND CLASSIFICATION FOR THE SCALING PARAMETER
		VALUES s = 1, s = 500 M AND s = 0.5 M.
	s	Center freq. Bandwidth	T	λ in D	Class.
	1	15.5	15.5	0.1	0.03	Unitless
	500 m	10 MHz	10 MHz	2.2 µs	15 m Radio wave
	0.5 m	10 GHz	10 GHz 2.2 ns	15 mm	Microwave
	.

TABLE IV THE

 IV RESULTS OF THE NUMERICAL EXPERIMENTS. THE BEST VALUE OF SSIM IS 1 AND THE BEST VALUE OF MSE AND ROE IS 0.

	Configuration BA order PPSNR (dB) SSIM Global MSE Void MSE Surface MSE Void ROE Surface ROE
	Monostatic	1	13.9	0.896	0.344	0.147	0.102	37.1	22.6
		2	13.9	0.914	0.283	0.171	0.0485	37.3	26.0
		3	13.9	0.918	0.288	0.189	0.0495	34.2	32.0
	Bistatic I	1	15.2	0.897	0.419	0.164	0.124	33.6	17.9
		2	15.2	0.916	0.302	0.173	0.0425	39.4	20.6
		3	15.2	0.922	0.285	0.185	0.0387	40.0	23.8
	Bistatic II	1	13.9	0.907	0.339	0.140	0.0647	34.2	26.4
		2	13.9	0.916	0.305	0.181	0.0480	38.0	32.6
		3	13.9	0.920	0.303	0.200	0.0493	39.8	36.0
	Multistatic	1	-15.2	0.908	0.489	0.232	0.0636	37.7	14.5
		2	-15.2	0.924	0.331	0.208	0.0334	36.5	18.0
		3	-15.2	0.925	0.305	0.209	0.0392	35.9	22.0
	TABLE V
	COMPUTING TIMES (S) FOR THE GPU-ACCELERATED ALGORITHMS.
		Process type
		BA order
		1st 2nd	3rd Other
	Jacobian matrix	45
	Wave-field update	75 150 225
	Reconstruction process 285 435 585
	Exact wave data			135
	Background wave data			45