# Frequency-domain modelling and inversion of electromagnetic data for 2D permittivity and conductivity imaging: An application to the Institut Fresnel experimental dataset

François Lavoué, Romain Brossier, Ludovic Métivier, Stéphane Garambois, Jean Virieux

<figure><img src="image_1.png" /></figure>

► To cite this version:

François Lavoué, Romain Brossier, Ludovic Métivier, Stéphane Garambois, Jean Virieux. Frequency-domain modelling and inversion of electromagnetic data for 2D permittivity and conductivity imaging: An application to the Institut Fresnel experimental dataset. Near Surface Geophysics, 2015, 13 (3), pp.227-241. 10.3997/1873-0604.2015004. hal-02009471

HAL Id: hal-02009471

https://hal.science/hal-02009471v1

Submitted on 3 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Frequency-domain modelling and inversion of electromagnetic data for 2D permittivity and conductivity imaging: An application to the Institut Fresnel experimental dataset

François Lavoué<fnref n="1" />, Romain Brossier<fnref n="1" />, Ludovic Métivier<fnref n="1" />,<fnref n="2" />, Stéphane Garambois<fnref n="1" />, and Jean Virieux<fnref n="1" />

Accepted 2014 October 12. Received 2014 August 1; in original form 2013 October 31.

## ABSTRACT

The need for a quantitative imaging of the near subsurface leads to the development of inversion algorithms to infer ground properties from recorded data. The aim of this study is to validate an inversion method recently developed for the simultaneous imaging of dielectric permittivity and electrical conductivity from 2D ground-penetrating radar measurements. The validation is performed using electromagnetic data collected in a well-controlled laboratory environment. In this experiment, the knowledge of the targets enables a quality control of the inversion results. In addition, the free space environment and the measurement of the incident field simplify the choice of a starting model for the inversion, as well as the calibration of the data with respect to the source signature and to the geometrical spreading. To perform accurate and efficient forward simulations, we use a frequency-domain finite-difference scheme whose stencil coefficients can be optimized for each simulated frequency. As the objects of interest are locally concentrated at the center of the acquisition array, it is possible to restrict the computation domain to a small region enclosing the targets using an integral representation of the analytical incident field coming from the sources and of the scattered field that we analytically propagate towards the receivers. An analysis of the numerical errors done on synthetic data shows that this strategy provides an error level that is low enough to not perturb the inversion, while dramatically decreasing the computational cost compared to a full-domain simulation. The monoparameter reconstruction of a purely dielectric target recovers permittivity values in very good agreement with the expected ones, as well as a very satisfying data fit. We also validate our strategy for multiparameter inversion on targets involving both a purely dielectric cylinder and a purely metallic copper tube, although the optimization cannot recover the exact conductivity of copper.

**Key words:** Ground penetrating radar (GPR), full waveform inversion (FWI), microwave imaging, 2D electromagnetic modelling, multiparameter inversion, controlled experiment.

## INTRODUCTION

The quantitative imaging of the near subsurface through non-invasive prospecting techniques appears as a crucial challenge in many fields of application. Among the geophysical methods available to explore the subsurface, ground-penetrating radar (GPR) is particularly interesting for its high resolution properties, despite its depth penetration limits. In geological, hydrological and geotechnical investigations, GPR can provide a qualitative view of the geometry of the sounded medium as well as a quantitative information

on the dielectric permittivity $\varepsilon$ [F/m] and on the electrical conductivity $\sigma$ [S/m] in the subsurface. These properties can in turn be interpreted in terms of the composition of the material (Deeds and Bradford 2002; Ihamouten et al. 2012) or of its water content (Garambois et al. 2002; Huisman et al. 2003; Day-Lewis et al. 2005; Weihermüller et al. 2007).

Recently, Lavoué et al. (2013, 2014) proposed an imaging method based on the full waveform inversion (FWI) of GPR data for the reconstruction of permittivity and conductivity in 2D sections of the subsurface. FWI is a state-of-the-

<footnote n="1">¹ISTerre, Univ. Grenoble Alpes - CNRS, 38041 Grenoble, France. E-mail: francois.lavoue@ujf-grenoble.fr</footnote>

<footnote n="2">²LJK, Univ. Grenoble Alpes - CNRS, 38041 Grenoble, France</footnote>
art quantitative imaging that aims at exploiting the whole information of the radargrams. Recent applications of FWI to GPR data have already been efficient for water content estimation in the first centimeters of agricultural soils (e.g. Lambot et al. 2006; Minet et al. 2010) and for the estimation of permittivity and conductivity in stratified structures such as concrete (Kalogeropoulos et al. 2011; Patriarca et al. 2011) or layered soils (Busch et al. 2012). FWI has also been applied to cross-hole radar data for the imaging of permittivity and conductivity in 2D and pseudo-3D (e.g. Ernst et al. 2007; Meles et al. 2010; Ellefsen et al. 2011; Klotzsche et al. 2013). In their previous contributions, Lavoué et al. (2013, 2014) perform numerical analysis on synthetic data to address the problem of multiparameter imaging of permittivity and conductivity by FWI in surface-to-surface acquisition configuration (on-ground GPR) that leads to a decreased illumination of the targets in the subsurface. To tackle the multiparameter FWI problem, they propose a joint optimization in the permittivity-conductivity parameter space. In this approach, permittivity and conductivity models are updated using a quasi-Newton scheme that considers the effect of an approximated Hessian matrix on the descent direction. The Hessian matrix (second-order derivatives of the misfit function) is expected to be crucial for multiparameter problems as it accounts for parameter dimensionalities and for possible trade-offs between different parameter types (Pratt et al. 1998; Operto et al. 2013). The strategy of Lavoué et al. (2013, 2014) involves parameter scaling and regularization factors to determine optimal solutions for the permittivity and conductivity distributions, based on the analysis of the data fit. In the present study, we would like to validate this strategy for multiparameter imaging by considering the inversion of well-controlled data acquired in a laboratory environment that provides a complete illumination of the targets.

The experimental scattering database collected by the Institut Fresnel (Marseille, France) is an interesting tool for testing and validating inversion algorithms on real physical data (see the dedicated special sections Belkebir and Saillard 2001, 2005; Litman and Crocco 2009). In this work, we focus on a dataset acquired on inhomogeneous targets (Belkebir and Saillard 2005), with a particular interest for those involving both dielectric and metallic objects which enable to address the problem of multiparameter reconstruction. The targets are cylinders located in free space (i.e., surrounded by air) and elongated perpendicularly to the observation plane, so that the problem can be considered to be 2D (see Geffrin et al. 2005, for the description of the experimental setup). The use of such experimental data greatly simplifies some issues that usually constitute major obstacles in FWI applications. First, the free space environment eliminates the crucial question of designing a suitable initial model. Second, the measurement of the incident field (in the absence of the targets) enables to accurately characterize the source signature and to calibrate the recorded data.

The present paper is organized in two parts, related to the two components of our imaging algorithm, namely the forward and the inverse problems. A preliminary but important aspect of the validation task concerns the design of the forward modelling that should accurately and efficiently reproduce the observed data. Our frequency-domain finite-difference (FDFD) modelling uses a mixed-grid stencil with

optimized coefficients (Jo et al. 1996; Hustedt et al. 2004) that can simulate the total field in the entire domain without significant errors due to numerical dispersion and numerical anisotropy, but it implies a significant computational effort. In a free space environment, it is highly desirable to restrict the computational domain to a small region enclosing the targets. Following the work of Wilcox and Velichko (2010), we use an integral representation of the fields to inject an analytical incident field in the reduced computation domain where FDFD computations are performed. A similar integral can be used to propagate analytically the recorded scattered field back to the real receiver locations. However, the discretization of the integral leads to numerical errors whose magnitude must be compared with the dispersion errors of the pure numerical approach to assess the advantage of one strategy over the other. In a second part of the paper, we challenge our imaging algorithm by reconstructing first a purely dielectric target (monoparameter inversion) and secondly a target containing both dielectric and metallic objects (multiparameter inversion).

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 1.** (a) Acquisition setup and computation domain. The cross located at $\theta_s = 45^\circ$ indicates a source and the dots between $\theta_r = 105^\circ$ and $\theta_r = 345^\circ$ the corresponding receiver locations. The 80-cm-by-80-cm black box corresponds to the reduced domain for FDFD computation. (b) Zoom on the computation domain. The 30-cm-by-30-cm green box delimits the zone reconstructed during the inversion. The red circle indicates the loading contour (393 virtual sources) and the blue one the monitoring contour (472 virtual receivers).</figcaption></figure>

PRESENTATION OF THE DATA

Fig. 1(a) shows the configuration of the experiment in the (xy)-plane (see Geffrin et al. 2005, their fig. 1, for a 3D representation). The targets, located at the center of the acquisition array, consist of cylinders elongated in the z-direction. Several experiments were performed for different
targets (Geffrin et al. 2005). For each experiment, the transmitting and receiving antennas are moved on a circular frame of radius $d_e = 1.67$ m around the targets. Fig. 1(a) exemplifies one source location (black cross) and the corresponding receiver locations (black dots), sweeping successively every 1°. Following Geffrin et al. (2005), we denote $\theta_s$ and $\theta_r$ the angular positions of sources and receivers, respectively, and we call angle between source and receiver the angle $|\theta_s - \theta_r|$ formed by the source, the target and the receiver (see Fig. 1a). During the experiment, the array is rotated to acquire data at different source positions, the angle between source and receiver always being comprised between 60° and 300°. For each source and receiver position, two measurements are performed: one for the incident field $u_{inc}$ in free space (without the target) and one for the total field $u_{tot}$ (with the target), enabling to deduce the field scattered by the targets $u_{sc} = u_{tot} - u_{inc}$. Finally, measurements for the two modes TE and TM are available. In this study, we consider only the measurements performed in TM mode which involves the electric field component $E_z$ polarized perpendicularly to the observation plane ($xy$). Please note that we follow the same terminology as Geffrin et al. (2005) by adopting the convention of the electromagnetic community for the TM mode (it is the TM$_z$ mode defined in Taflove and Hagness 2005, p. 55). It corresponds to the TE mode for the geophysical community (see e.g. Jol 2009, fig. 1.8, p. 13).

Fig. 2 shows synthetic models corresponding to the targets we will consider in the following, labelled as FoamDielExt (Fig. 2a) and FoamMetExt (Fig. 2b,c) by Gefrin et al. (2005). The target FoamDielExt consists of two purely dielectric cylinders made of foam ($\varepsilon_r = 1.45 \pm 0.15$) and of plastic ($\varepsilon_r = 3 \pm 0.3$), of diameters 8 cm and 3.1 cm, respectively. Their conductivity is considered to be zero. Geffrin et al. (2005) perform numerical modelling on these synthetic models using a method of moments and domain integral formulations, so that we can compare our simulation results. Regarding the inversion, we will use the experimental dataset FoamDielExt in order to validate the monoparameter reconstruction of permittivity. The target FoamMetExt consists of the same purely dielectric foam cylinder ($\varepsilon_r = 1.45 \pm 0.15$) and of a purely metallic copper tube ($\varepsilon = 1$, $\sigma \gg 1$ S/m) of diameter 2.85 cm and of thickness 2 mm. We will use this dataset to challenge our strategy for multiparameter inversion. Of course, we cannot hope to reconstruct the exact amplitude of the conductivity anomaly (which is around $10^6$ S/m for pure copper). But we can analyse the trade-off effect between permittivity and conductivity, i.e. the trend to recover erroneous non-zero conductivity values for the foam cylinder, or permittivity values $\varepsilon_r \neq 1$ for the copper.

Table 1 sums up the measurement properties. For the target FoamDielExt, measurements are performed for 9 frequencies between 2 and 10 GHz, and repeated over 8 source positions with an interval of 45°. For the target FoamMetExt, measurements are performed for 17 frequencies between 2 and 18 GHz, and repeated over 18 source positions with an interval of 20°. In the latter case, the extension of the frequency bandwidth and the refinement of the source sampling aim at better constraining the multiparameter inversion for the reconstruction of the complex target FoamMetExt (Geffrin et al. 2005).

In the frequency range we consider, target sizes are of the order of — or even below — the smallest wavelength in free space, while the propagation distance goes beyond 100 wavelengths (see Table 2). Two aspects may require our attention. First, we cannot expect to recover the exact thickness of the copper tube which is of the order of $\lambda/8$ for the highest available frequency, given that the maximal resolution of our imaging algorithm, based on the diffraction principle, is of the order of $\lambda/2$. Secondly, we can expect to face a numerical challenge to simulate the wave propagation over 100 wavelengths between the target and the receivers (even over 200 wavelengths on the path source-target-receiver).

## FORWARD PROBLEM

### Numerical strategy

A first challenge consists in accurately and efficiently performing synthetic simulations to be compared with experimental data in the inversion process. In this study, we use a FDFD scheme based on the optimized mixed-grid stencil of Hustedt et al. (2004), which leads to the linear system

$$A(\omega, \mathbf{m}) \mathbf{u}(\omega) = \mathbf{S}(\omega), \quad (1)$$

where $\omega$ is the simulated angular frequency (in rad/s), $\mathbf{m}$ represents the physical model, $A$ is the impedance matrix resulting from the FDFD scheme, $\mathbf{u}$ is the simulated wavefield (i.e. the component $E_z$ in TM mode) and $\mathbf{S}$ denotes the source term. We solve the linear system (1) for each frequency through a LU factorization using the direct solver MUMPS (MUMPS-team 2011).

In the finite-difference scheme, the optimization of the weighting coefficients of the stencil for a given $\lambda/h$ ratio (with $h$ the grid step) enables to minimize the errors due to numerical dispersion and numerical anisotropy (Jo et al. 1996; Hustedt et al. 2004). It is thus possible to achieve an error of less than 1.5% on the incident field recorded by the receivers. However, the computation in the entire domain of Fig. 1(a) implies a significant computational effort, both in terms of CPU time and of memory requirement, especially because the large number of degrees of freedom requires the use of double precision arithmetic to avoid instabilities in MUMPS. In view of the inversion, the efficiency of the forward problem is crucial. In a free space environment, it is highly desirable to restrict the computational domain to a small region enclosing the target, as shown in Fig. 1(b).

To perform the inversion in this reduced domain, we must be able:

(i) to inject in the domain an arbitrary incoming wavefield emitted from a remote source of the real acquisition array located outside the domain,

(ii) to propagate analytically the field scattered by the target and computed numerically in the reduced domain back to the real receiver locations outside the domain.

To do so, we use an integral representation of the fields, following the work of Wilcox and Velichko (2010) and Velichko and Wilcox (2010). Integral representations are usual in non-destructive testing applications (NDT). For instance, it has been recently used by Zhao et al. (2013) for GPR data re-datuming in a tunnel grouting test context. In our case, it
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Synthetic models corresponding to the experimental targets. (a) FoamDielExt (permittivity model). (b,c) FoamMetExt (b, permittivity; c, conductivity). In (a) and (b), the central foam cylinder has a permittivity $\epsilon_r = 1.45$.</figcaption></figure>

<table><caption>**Table 1.** Measurements properties (after Geffrin et al. 2005, their table 1).</caption><thead><tr><th rowspan="2">Dataset</th><th colspan="2">Frequencies</th><th colspan="2">Sources</th><th colspan="2">Receivers</th></tr><tr><th>Nb</th><th>Min-max</th><th>Interval</th><th>Nb</th><th>Interval</th><th>Nb</th><th>Interval</th></tr></thead><tbody><tr><td>FoamDielExt</td><td>9</td><td>2-10 GHz</td><td>1 GHz</td><td>8</td><td>45°</td><td>241</td><td>1°</td></tr><tr><td>FoamMetExt</td><td>17</td><td>2-18 GHz</td><td>1 GHz</td><td>18</td><td>20°</td><td>241</td><td>1°</td></tr></tbody></table>

is particularly suited because we know the exact analytical solution for the wave propagation in free space between the acquisition array and the target region. The Helmholtz-Kirchhoff integral provides an expression for the scalar wave-field $u$ at any point $\mathbf{r}$ within a closed contour $C$, as a function of an arbitrary field $u_{inc}$ on the boundary:

$$u(\mathbf{r}) = \int_C [G(\mathbf{r}, \mathbf{s}) \nabla u_{inc}(\mathbf{s}) - u_{inc}(\mathbf{s}) \nabla G(\mathbf{r}, \mathbf{s})] \cdot \mathbf{n}dS, \quad (2)$$

where $\mathbf{n}$ is the outward normal-contour vector and $G(\mathbf{r}, \mathbf{s})$ denotes the Green function, i.e. the field value at point $\mathbf{r}$, resulting from a source located at $\mathbf{s}$. The vector $\mathbf{r}$ locates a point inside the contour $C$, while the vector $\mathbf{s}$ points on the contour. In other words, we can deduce the field $u$ inside the contour from the superposition of fields generated by monopole sources ($G$ terms) and dipole sources ($\nabla G$ terms) located on the contour. In the following, we will refer to these sources as *virtual* sources (in opposition with real ones, located on the acquisition array outside the computation domain, see Fig. 1a). In expression (2), the Green function $G(\mathbf{r}, \mathbf{s})$ is computed numerically with the FDFD scheme. The values of the incident field on the contour $u_{inc}(\mathbf{s})$ are known analytically as the solution of 2D wave propagation in free space between the real and the virtual source arrays (see e.g. Taflove and Hagness 2005, §5.3.1, p.172).

A similar integral representation can be used to describe the scattered field $u_{sc}(\mathbf{r})$ outside the contour $C$ if the scattered field on the contour $u_{sc}(\mathbf{s})$ is known (Velichko and Wilcox 2010):

$$u_{sc}(\mathbf{r}) = - \int_C [G(\mathbf{r}, \mathbf{s}) \nabla u_{sc}(\mathbf{s}) - u_{sc}(\mathbf{s}) \nabla G(\mathbf{r}, \mathbf{s})] \cdot \mathbf{n}dS, \quad (3)$$

where the vector $\mathbf{r}$ locates now a point outside the contour. In eq. (3), the scattered field $u_{sc}(\mathbf{s})$ on the contour is computed numerically with the FDFD scheme (it corresponds in fact to the field $u(\mathbf{r})$ in eq. 2), while the propagator $G(\mathbf{r}, \mathbf{s})$ and its gradient across the contour $\nabla G(\mathbf{r}, \mathbf{s})$ are computed

analytically. The minus sign before the integral is due to the orientation convention of the outward-pointing vector $\mathbf{n}$.

We use the windowed Sinc interpolation proposed by Hicks (2002) to locate accurately the virtual sources on the integration contour, that does not coincide with the finite-difference grid. This interpolation consists in spreading the source excitation on several grid nodes around the true source location to mimic a bandlimited version of a monopole or dipole point source. Because the implementation of point sources creates singularities in the numerical solution, it is not possible to record the scattered field $u_{sc}(\mathbf{s})$ of eq. (3) in the neighbourhood of nodes where the incident field $u_{inc}(\mathbf{s})$ is injected (eq. 2). Following Velichko and Wilcox (2010), we therefore implement two distinct contours: a *loading* contour $C_l$ to inject the incident field and a *monitoring* contour $C_m$, located outside the loading contour, to record the scattered field and apply eq. (3) to propagate it towards the receiver in the far-field region (see Fig. 1).

The discretization of the loading contour $C_l$ in $N_l$ segments of length $\delta C$ then leads eq. (2) to become

$$u(\mathbf{r}) = \delta C \sum_{i=1}^{N_l} \left[ \cos(\theta_i) \frac{\partial u_{inc}}{\partial x}(\mathbf{s}_i) + \sin(\theta_i) \frac{\partial u_{inc}}{\partial z}(\mathbf{s}_i) \right] G(\mathbf{r}, \mathbf{s}_i) \\ - u_{inc}(\mathbf{s}_i) \left[ \cos(\theta_i) D_x + \sin(\theta_i) D_z \right] G(\mathbf{r}, \mathbf{s}_i), \quad (4)$$

where the operators $D_x$ and $D_z$ correspond to x- and z-oriented dipoles using the windowed Sinc interpolation and $\theta_i$ is the angular position of the $i^{th}$ virtual source. The partial derivatives of the incident field $u_{inc}$ are computed analytically.

Conversely, when implementing eq. (3), the windowed Sinc interpolation is used to extract values of the scattered field at virtual receivers on the monitoring contour $C_m$ and the derivatives of the propagator $G$ are computed analytically.
<table><caption>**Table 2.** Target sizes and propagation distances in terms of wavelengths λ propagated in free space.</caption><thead><tr><th rowspan="2">Objects/lengths</th><th rowspan="2">(cm)</th><th>2 GHz</th><th>10 GHz</th><th>18 GHz</th></tr><tr><th>λ = 15 cm</th><th>λ = 3 cm</th><th>λ = 1.7 cm</th></tr></thead><tbody><tr><td>Grid step h</td><td>0.1</td><td>λ/150</td><td>λ/30</td><td>λ/16</td></tr><tr><td>Copper thickness</td><td>0.2</td><td>λ/75</td><td>λ/15</td><td>λ/8</td></tr><tr><td>Copper diameter</td><td>2.85</td><td>~ λ/5</td><td>~ λ</td><td>~ 1.8λ</td></tr><tr><td>Plastic diameter</td><td>3.1</td><td></td><td></td><td></td></tr><tr><td>Foam diameter</td><td>8.0</td><td>λ/2</td><td>2.7λ</td><td>5λ</td></tr><tr><td>Distance target-receiver</td><td>167</td><td>11λ</td><td>56λ</td><td>100λ</td></tr><tr><td>Distance source-target-receiver</td><td>334</td><td>22λ</td><td>112λ</td><td>200λ</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** (a,c) Wavefields emitted by the virtual source array (loading contour, in red). The blue circle indicates the monitoring contour. (b,d) Relative difference with the corresponding analytical wavefields. (a,b) 2 GHz, (c,d) 18 GHz. Excepted at nodes where sources are injected, the errors are below 0.47% at 2 GHz (b) and below 0.29% at 18 GHz (d) (to be compared with the errors done when performing the FDFD computation in the entire domain, which are below 0.19% and 0.14%, respectively). Note that the error ϵ is not computed in the same way inside and outside the loading contour: inside the contour, $ϵ = |u_{numerical} - u_{analytical}|/|u_{analytical}|$, and outside, $ϵ = |u_{numerical}|/|u_{analytical}|$ (because $u_{numerical}$ should be zero).</figcaption></figure>

cally:

$$u_{sc}(\mathbf{r}) = -\delta C \sum_{i=1}^{N_m} \left[ \cos(\theta_i) D_x + \sin(\theta_i) D_z \right] u_{inc}(\mathbf{s}_i) G(\mathbf{r}, \mathbf{s}_i) - u_{inc}(\mathbf{s}_i) \left[ \cos(\theta_i) \frac{\partial G}{\partial x}(\mathbf{r}, \mathbf{s}_i) + \sin(\theta_i) \frac{\partial G}{\partial z}(\mathbf{r}, \mathbf{s}_i) \right]. (5)$$

Figs 3(a) and (c) show the incident fields resulting from eq. (4) at frequencies 2 and 18 GHz respectively, using a contour element $\delta C = 4$ mm (= $\lambda/4$ at 18 GHz). Perfectly Matched Layers (PMLs, Bérenger 1994) are used to absorb the waves on the edges of the computation domain (not shown). In Figs 3(b) and (d), the relative errors against analytical solution are shown. Using the integral representation to inject the incident field in the reduced domain, we achieve an error level which is slightly higher than when computing the field in the entire domain ($\simeq 0.5\%$ vs. $\simeq 0.2\%$ respectively), but it still remains of the same order of magnitude and we consider it as satisfactory.

An important property of the integral representation, eq. (2), is that the incident field synthesized by the virtual array exists only inside the contour C, where eq. (2) is valid (see Fig. 3). However, if a scatterer is placed inside the contour, then the scattered field generated by the incident field of eq. (2) exists both inside and outside the contour (because the scatterer acts as an internal source). This property is illustrated in Fig. 4, where we inject the incident field of Fig. 3(a) in a domain that contains the target FoamDielExt of Fig. 2(a). The scattered field is then analytically propagated using eq. (3) from the monitoring contour (blue circle in Fig. 4) up to the receiver positions on the acquisition ar-

<figure><img src="image_7.png" /><figcaption>**Figure 4.** Example of wavefield scattered by the synthetic target FoamDielExt of Fig. 2(a), illuminated by the incident field of Fig. 3(a), at 2 GHz. Inside the loading contour (red circle), the total field is computed, whereas only the diffracted field propagates outside and is recorded by the monitoring contour (blue circle).</figcaption></figure>

ray (see Fig. 1a). Note that if the monitoring contour were placed inside the loading contour, recording the total field, we would still propagate only the scattered field outwards the contour using eq. (3) (this is a property of the Helmholtz-Kirchhoff integral). We place the monitoring contour outside the loading contour to avoid the cumulative errors of both integrals on the propagation of the incident field.

In order to evaluate the numerical error done on the scattered field computed in the reduced domain and propagated towards the receivers, we perform the following simple test: We put a point source at the center of the acquisition ar-
<table><caption>**Table 3.** Computation cost and accuracy of full-domain vs. reduced-domain FDFD simulations. The LU factorization and the resolution of the linear system are performed by MUMPS in parallel on 16 MPI processes. We report elapsed times using one node of 16 cores of the Froggy cluster from the CIMENT platform (Univ. Grenoble). Note that the full-domain computation is performed with double-precision arithmetic while the reduced-domain computation uses single precision to achieve the level of accuracy indicated in the table.</caption><thead><tr><th></th><th>Full domain</th><th>Reduced domain</th></tr></thead><tbody><tr><td>Number of grid points (including PMLs)</td><td>4051 × 4051</td><td>851 × 851</td></tr><tr><td>Elapsed time for LU factorization</td><td>71 s</td><td>0.67 s</td></tr><tr><td>Time for FoamMetExt simulation (18 src, 17 freq.)</td><td>2100 s</td><td>98 s</td></tr><tr><td>MUMPS memory requirement (LU)</td><td>59 Go</td><td>1.5 Go</td></tr><tr><td>Error on the incident field in the reduced domain</td><td>&lt; 0.2%</td><td>&lt; 0.5%</td></tr><tr><td>Error on the scattered field at receiver positions</td><td>&lt; 1.5%</td><td>&lt; 0.3% (0.8%)</td></tr></tbody></table>

frame (instead of the target), we record the resulting field on the monitoring contour, and we propagate it analytically towards the receivers on the acquisition array. We then compare the values obtained at receivers with the analytical solution and we find an error below 0.3%. Note that, when an incident field is injected on the loading contour and impinges a target, the error on the scattered field resulting from eq. (3) should also include the error done on the injected incident field (≈ 0.5% after Fig. 3). So we can estimate the global error done on the scattered field to approximately 0.8%. To compare this error with the full-domain case, we have to consider that, when computing the scattered field in the entire domain, the solution suffers from errors due to numerical dispersion and numerical anisotropy on the entire travel path source-target-receiver (of length ≈ 2d<sub>e</sub>). As a result, in the full-domain case, we consider that the error on the incident data recorded in transmission configuration (that we can compute against analytical solution) is also representative for the error done on the scattered data.

These results are reported in Table 3, together with the associated computational requirements. Regarding the errors done on the incident and scattered fields, the restriction of the computation in the reduced domain using the integral representation appears as a satisfying alternative to the full-domain computation, decreasing the computational time by a factor of more than 20, while dividing the memory requirement of the LU factorization by a factor of 40.

## Data pre-processing

To be able to compare our synthetic data with the observed ones, we have to consider that observed and calculated data differ by three major aspects:

(i) Observed data are 3D whereas our modelling is 2D.

(ii) The antennas used in the experiment have a given frequency signature, whereas the simulated source is a Dirac in time (unity source in the frequency domain).

(iii) The experimental setup use ridged horn antennas with a given radiation pattern whereas the modelling assumes an elementary dipole oriented in the z-direction which radiates isotropically in the (xy)-plane.

These characteristics have to be accounted for to convert the observed data in a value that we can compare with the synthetic one.

Usually, the estimation of the source complex spectrum is performed as a part of the frequency-domain inversion process by linear estimation (Pratt 1999). In our case, however,

it is highly beneficial to apply this estimation to the incident field measured in free space, which enables us both to characterize the antennas and to perform a simple 3D-to-2D conversion. To do so, let us write the observed data $d_{inc}^{obs}$ resulting from the measurements of the incident field, for each angular frequency $\omega$, each source located at angle $\theta_s$, and each receiver at $\theta_r$, as

$$d_{inc}^{obs}(\omega, \theta_s, \theta_r) = S_{obs}(\omega) R_{obs}(\omega, \theta_s - \theta_r) G_{3D}(\omega, d_{sr}), \quad (6)$$

where $S_{obs}(\omega)$ denotes the source complex spectrum, $R_{obs}(\omega, \theta_s - \theta_r)$ accounts for the spatial radiation pattern of the antennas, and $G_{3D}$ is the Green function of electromagnetic wave propagation in 3D, that varies with the distance $d_{sr}$ travelled by the wave from source to receiver.

If we adopt the convention that, for each angular frequency $\omega$, the radiation pattern coefficients are 1 in pure transmission regime (i.e., $R_{obs}(\omega, \theta_o) = 1$ for an angle between source and receiver of $\theta_o = |\theta_s - \theta_r| = 180^{\circ}$), then we can deduce the source spectrum of the antennas from the incident field recorded in transmission configuration:

$$S_{obs}(\omega) = \frac{d_{inc}^{obs}(\omega, \theta_s, \theta_r = \theta_s - \theta_o)}{G_{3D}(\omega, d_{rs} = 2d_e)}. \quad (7)$$

Fig. 5 shows the obtained source spectrum, both in amplitude and phase. Computing the spectrum for all source positions, we find that its variations with respect to source location $\theta_s$ are negligible (of the order of 0.01%), showing the quality of these highly reproducible data.

Once we know the spectrum $S_{obs}(\omega)$ from the measurements in pure transmission configuration, we could deduce the radiation pattern of the antennas from the incident field recorded at the other incidence angles. However, we do not need it for our present study because we only consider the scattered field for the inversion. As the scattered field results from an incident field impinging the target at zero incidence (from the source point of view), and as it is itself recorded at zero incidence by the receiver, we can reasonably neglect the effect of the radiation pattern on the scattered data.

A last step to make the observed and calculated fields comparable is to correct for the 3D geometrical expansion, in order to build 2D observed data. By analogy with expression (6), we assume that the observed and synthetic scattered data can be written as

$$d_{sc}^{obs}(\omega, \theta_s, \theta_r) = S_{obs}(\omega) T(\omega, \theta_s, \theta_r) G_{3D}(\omega, d_{str}), \quad (8)$$

$$d_{sc}^{cal}(\omega, \theta_s, \theta_r) = S_{cal}(\omega) T(\omega, \theta_s, \theta_r) G_{2D}(\omega, d_{str}), \quad (9)$$

where $S_{cal}(\omega)$ is a known factor depending on our FDFD implementation, $G_{2D}$ is the 2D Green function which is known
<figure><img src="image_8.png" /></figure>

<figure><img src="image_9.png" /><figcaption>Figure 5. (a) Amplitude and (b) phase of the antenna complex spectrum $S_{obs}(\omega)$.</figcaption></figure>

analytically, $d_{str}$ is the distance travelled by the scattered wave from source to target and from target to receiver, which can be approximated by $d_{str} \simeq 2d_e$, and T is the response of the targets, which can be considered identical in 2D and 3D since the targets have a 2D geometry. Thus, the transformation to be applied to the observed data to compare them to the synthetic ones is

$$d_{sc}^{\text{obs, 2D}}(\omega, \theta_s, \theta_r) = \frac{S_{\text{cal}}(\omega)}{S_{\text{obs}}(\omega)} \frac{G_{2D}(\omega, 2d_e)}{G_{3D}(\omega, 2d_e)} d_{sc}^{\text{obs}}(\omega, \theta_s, \theta_r), \quad (10)$$

where the 3D-to-2D conversion factor $k_{3D \to 2D} = G_{2D}/G_{3D}$ can be expressed by using the far-field approximations of the Green functions (see Taflove and Hagness 2005, §8.2.2, p.332),

$$k_{3D \to 2D} = 2 \sqrt{\frac{\pi d_e c_o}{\omega}} e^{i\pi/4}, \quad (11)$$

with $c_o$ the velocity of light in free space and $i$ the imaginary unit.

## Simulation of synthetic data

Simulation results performed in the synthetic models of Fig. 2 are shown in Figs 6 and 7 for the targets FoamDielExt and FoamMetExt, respectively. We choose the same frequencies and source positions as shown by Geffrin et al. (2005, their figs. 8 and 11) such that we can compare our results. In view of the inversion, we also perform numerical simulations in modified versions of the synthetic models of Fig. 2,

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Figure 6. (a) Amplitude and (b) phase of synthetic vs. observed data for the target FoamDielExt, for frequency 8 GHz and the source at $\theta_s = 270^{\circ}$. The so-called corrected field (blue dashed line) corresponds to the synthetic field computed in a modified version of the model shown in Fig. 2(a) where we translated the targets by $dx = 2$ mm and $dy = 1$ mm and we replaced the permittivity value of the plastic cylinder by 3.3 instead of 3. The match between the observed diffracted field and the synthetic one is similar as in Geffrin et al. (2005, their fig. 8).</figcaption></figure>

in order to investigate the sensitivity of the data to model variations.

For the target FoamDielExt, our simulation results are very similar to those of Geffrin et al. (2005). As suggested by these authors, we also perform a simulation in a modified version of the model shown in Fig. 2(a), where we translated the targets by $dx = 2$ mm and $dy = 1$ mm, and we replaced the permittivity value of the plastic cylinder by 3.3 instead of 3 to obtain a better match (blue dashed line in Fig. 6).

For the target FoamMetExt (Figs 2b,c), simulation results at 18 GHz are shown in Fig. 7 (black line). In this case, our simulation results differ slightly from those presented by Geffrin et al. (2005, their fig. 11): Contrary to these authors, who calibrate the amplitude of their numerical solution such as to match the central diffraction peak (at 180°), we do not recover the same amplitude for the observed and for the synthetic data at 180°. Nonetheless, we obtain a similar match as Geffrin et al. (2005) for the other angles. These differences in amplitude with respect to angle are probably due to our 3D-to-2D conversion which is based on the incident field in free space (without the target): It may be thus poorly valid to convert the signal recorded at 180°, which is transmitted through the target. Again, we perform simulations in a modified version of the FoamMetExt models of Figs 2(b,c),
<figure><img src="image_12.png" /><figcaption>Figure 7. (a) Amplitude and (b) phase of synthetic vs. observed data for the target FoamMetExt, for frequency 18 GHz and the source at $\theta_s = 240^{\circ}$. The match obtained with $\sigma = 5$ S/m for copper (black line) is similar as in Geffrin et al. (2005, their fig. 11) for extreme angles but we do not recover the same amplitude as these authors at 180$^{\circ}$. A better match is obtained with $\sigma = 100$ S/m (blue dashed line).</figcaption></figure>

where we replace the conductivity value of the copper tube by 100 S/m instead of 5 S/m (blue dashed line). Doing so, we obtain a better match at extreme angles. It shows that the data are well sensitive to conductivity variations in the copper tube and thus suggests that a quantitative reconstruction can be attempted, even if we cannot expect to recover the true value of conductivity in copper ($\simeq 10^6$ S/m). As another sensitivity test, we also compared the synthetic data of Fig 7 with data computed in a model where we replaced the empty copper tube of Fig. 2(c) by a full copper cylinder of same diameter and properties. An important result is that the full cylinder provides exactly the same scattered field as the empty tube when a high conductivity value is used for copper ($\sigma = 100$ S/m), suggesting that any information inside the real tube will be out of reach. Assuming a lower (non-physical) conductivity value for copper ($\sigma = 5$ S/m), data are slightly sensitive to the filling of the tube, which will have consequences in the reconstructions.

In Fig. 8, we compare in the case FoamDielExt the magnitude of the data residuals, i.e. the difference between observed and synthetic data, with the numerical errors, computed as the difference between the analytical solution of 2D wave propagation in free space and the synthetic incident field computed in the entire domain. According to the numerical analysis of the previous section, this error is an

upper bound for the error done on the synthetic scattered field that will be involved in the inversion (see Table 3). For information, we also indicate the magnitude of the total and scattered fields: As underlined by Geffrin et al. (2005), the high ratio between the magnitude of the scattered and total fields makes the need for an accurate modelling critical. In Fig. 8, the error level is significantly lower than the data residuals for most of the source-receiver pairs (it is also true for other sources and frequencies, not shown in Fig. 8), which means that most of the residuals are due to differences between the synthetic model and the reality. This result suggests that our modelling tool fills the accuracy requirement: The inversion should not be affected by the numerical errors. The very good agreement between observed and synthetic data in Fig. 6, as well as the low error level compared to the residuals (Fig. 8), allow us to envisage the inversion of these data with some confidence.

## DATA INVERSION

### Inverse problem formulation

Here we recall the basic ingredients of the method proposed by Lavoué et al. (2014). At each iteration $k$ of the optimization process, we aim at minimizing a misfit function that we define as the sum over frequencies of a data part $\mathcal{C}_D$ and of a model part $\mathcal{C}_M$:

$$ \mathcal{C}(\mathbf{m}^k) = \mathcal{C}_D(\mathbf{m}^k) + \lambda \mathcal{C}_M(\mathbf{m}^k), \quad (12) $$

$$ = \frac{1}{2} \sum_{\omega} \Delta\mathbf{d}(\omega, \mathbf{m}^k)^{\dagger} \Delta\mathbf{d}(\omega, \mathbf{m}^k) + \frac{\lambda}{2} \mathbf{m}^k T \mathcal{D} \mathbf{m}^k. \quad (13) $$

The first term $\mathcal{C}_D$ accounts for the misfit between the observed and the synthetic data through the $l_2$-norm of the residuals $\Delta\mathbf{d}(\omega, \mathbf{m}^k) = \mathbf{d}^{\text{obs}}(\omega) - \mathbf{d}^{\text{cal}}(\omega, \mathbf{m}^k)$ that measure the difference between the observed data $\mathbf{d}^{\text{obs}}$ and the synthetic data $\mathbf{d}^{\text{cal}}$ computed in the current model $\mathbf{m}^k$. Usually, the residuals involve the total simulated wavefield $u_{\text{tot}}$. In our case, however, it is much more convenient to use the scattered field $u_{\text{sc}} = u_{\text{tot}} - u_{\text{inc}}$ since

(i) it can be computed using the measurements of the incident field $u_{\text{inc}}$,

(ii) it is insensitive to antenna radiation pattern,

(iii) it is the field that is analytically propagated toward the receivers by the integral representation (see eq. 3 and Fig. 4).

Assuming that the incident field can be accurately simulated (which has been demonstrated by our numerical analysis, see Figs 3 and 8), we thus have

$$ \Delta\mathbf{d}(\omega, \mathbf{m}^k) = \mathbf{d}_{\text{tot}}^{\text{obs}}(\omega) - \mathcal{R} \mathbf{u}_{\text{tot}}^{\text{cal}}(\omega, \mathbf{m}^k), \quad (14) $$

$$ = \mathbf{d}_{\text{sc}}^{\text{obs}}(\omega) - \mathcal{R} \mathbf{u}_{\text{sc}}^{\text{cal}}(\omega, \mathbf{m}^k), \quad (15) $$

where $\mathcal{R}$ is the projection operator of the wavefield on the receiver location, which takes the integral representation and the windowed Sinc interpolation into account.

The second term $\mathcal{C}_M$ of the misfit function (13) introduces a Tikhonov regularization (Tikhonov and Arsenin 1977). The operator $\mathcal{D}$ corresponds to the Laplacian so that its minimization tends to provide smooth solutions. The hyperparameter $\lambda$ is a regularization weight that balances the importance given to model smoothness relatively to the data
<figure><img src="image_13.png" /><figcaption>**Figure 8.** Comparison of the magnitude of the total and scattered fields vs. residuals and numerical errors, for frequency 8 GHz and the source at $\theta = 270^{\circ}$. The synthetic scattered field is computed in the modified version of Fig. 2(a). The numerical error is obtained by difference between the analytical solution and the total incident field computed in the entire domain.</figcaption></figure>

misfit. In eq. (13), the symbol $^T$ denotes transposition and the notation $^\dagger$ corresponds to the transpose$^T$-conjugate$^*$ operator.

To minimize the misfit function, we use a quasi-Newton optimization method (L-BFGS-B algorithm, Byrd et al. 1995) which performs a local descent based on the gradient of the misfit function $G(\mathbf{m}^k) = \nabla_{\mathbf{m}^k} \mathcal{C}(\mathbf{m}^k)$. The gradient of the model term is straightforwardly computed by finite differences on the model $\mathbf{m}^k$ whereas the gradient of the data term $G_D(\mathbf{m}^k)$ is computed via the adjoint state method (Plessix 2006). Differentiating eqs (13) and (14) with respect to the model parameters, we have

$$ 
\begin{aligned}
G_{D_i}(\mathbf{m}^k) &= \frac{\partial \mathcal{C}_D(\mathbf{m}^k)}{\partial m_i}, \\
&= \sum_{\omega} \mathcal{R}e \left\{ \frac{\partial \Delta \mathbf{d}(\omega, \mathbf{m}^k)^{\dagger}}{\partial m_i} \Delta \mathbf{d}(\omega, \mathbf{m}^k) \right\}, \\
&= -\sum_{\omega} \mathcal{R}e \left\{ \frac{\partial \mathbf{u}_{tot}^{cal}(\omega, \mathbf{m}^k)^{\dagger}}{\partial m_i} \mathcal{R}^{\dagger} \Delta \mathbf{d}(\omega, \mathbf{m}^k) \right\}, \quad (16)
\end{aligned}
$$

where $\mathbf{u}_{tot}^{cal}$ verifies the equation of the forward problem (1), which gives by differentiation

$$ \frac{\partial \mathbf{u}_{tot}^{cal}(\omega, \mathbf{m}^k)}{\partial m_i} = -\mathcal{A}(\omega, \mathbf{m}^k)^{-1} \frac{\partial \mathcal{A}(\omega, \mathbf{m}^k)}{\partial m_i} \mathbf{u}_{tot}^{cal}(\omega, \mathbf{m}^k). \quad (17) $$

Injecting eq. (17) in eq. (16) yields the expression of the gradient

$$ G_{D_i}(\mathbf{m}^k) = \sum_{\omega} \mathcal{R}e \left\{ \mathbf{u}_{tot}^{cal}(\omega, \mathbf{m}^k)^T \frac{\partial \mathcal{A}(\omega, \mathbf{m}^k)^T}{\partial m_i} \mathbf{v}^*(\omega, \mathbf{m}^k) \right\}, \quad (18) $$

where $\mathbf{v}$ is the adjoint wavefield, that verifies the linear system $\mathcal{A}^\dagger \mathbf{v} = \mathcal{R}^\dagger \Delta \mathbf{d}$ and corresponds to the back-propagation of the residuals in the current model (Plessix 2006).

In expression (18), the diffraction matrix $\partial_{m_i} \mathcal{A}$ characterizes the sensitivity of the data to the parameter $m_i$, that refers either to the permittivity $\varepsilon_i$ or to the conductivity $\sigma_i$ at grid point $i$. Lavoué et al. (2014) promote to consider dimensionless parameters that can be gathered in the same model vector $\mathbf{m}$, thus enabling to perform the optimization in the joint permittivity-conductivity parameter space (in opposition to cascaded or decoupled approaches as used for instance by Ernst et al. 2007; Meles et al. 2010; Ellefsen et al.

2011). Thus, we consider a relative permittivity $\varepsilon_r = \varepsilon/\varepsilon_o$ (with $\varepsilon_o \simeq 8.85 \times 10^{-12}$ F/m the dielectric constant) and a relative conductivity $\sigma_r = \sigma/\sigma_o$. The reference conductivity $\sigma_o$ can be defined as $\sigma_o = \varepsilon_o \omega_o$, by analogy with the contrast function used in the inverse scattering community (see e.g. Abubakar et al. 2005), with $\omega_o$ a reference angular frequency that we take as the mean frequency of the measurements (i.e., we consider $\omega_o = 2\pi \times 6$ GHz for the dataset FoamDielExt and $\omega_o = 2\pi \times 10$ GHz for FoamMetExt). However, these definitions are arbitrary (they are only conventions) and we pointed out in Lavoué et al. (2014) that the inversion is very sensitive to the priority given to the permittivity or to the conductivity update via the parameter scaling. Therefore, we introduce an additional scaling factor $\beta$ and we consider the couple of parameters $(\varepsilon_r, \sigma_r/\beta)$ for the optimization. As detailed in Lavoué et al. (2014), a small value for the factor $\beta$ will give priority to the reconstruction of permittivity, whereas a large scaling factor will enhance the conductivity update.

On synthetic experiments, Lavoué et al. (2014) propose a workflow to choose an adequate value for the scaling factor $\beta$, in conjunction with an appropriate regularization weight $\lambda$, based on the data misfit analysis. In the following, we aim at validating this workflow on the experimental data provided by the Institut Fresnel. In all the following inversion tests, we stop the optimization process when the relative misfit decrease between two iterations becomes smaller than $10^4$ times the machine precision (default stopping criterion in the L-BFGS-B algorithm, Zhu et al. 1997).

## Monoparameter inversion of the dataset FoamDielExt

First of all, we wish to validate the quantitative imaging of permittivity through monoparameter inversion, that does not involve any parameter scaling. The target FoamDielExt is particularly suited for this purpose as the objects are supposed to be purely dielectric and have well constrained permittivity values. Fig. 9 shows the inversion result when inverting the 9 frequencies between 2 and 10 GHz simultaneously, without any regularization ($\lambda = 0$). During the inversion, the misfit function has been decreased by 98.8% within 32 iterations. We recover permittivity values of $\varepsilon_r =$
<figure><img src="image_14.png" /><figcaption>Figure 9. Permittivity reconstruction of the target FoamDielExt by monoparameter inversion. (a) 2D model. The red dashed circles indicate the expected contours of the objects, after the translation by $dx = 2$ mm, $dy = 1$ mm suggested by Geffrin et al. (2005). (b) Logs extracted from the 2D models along the horizontal line $y = 0$. The reconstructed values (blue line) are compared to the synthetic model of Fig. 2(a) (red line) and to low-pass-filtered versions of this model, considering the full frequency bandwidth (black line) or only the lowest frequency (green line).</figcaption></figure>

1.42±0.05 for the foam cylinder and $ε_r = 3.10±0.42$ for the plastic one (mean and variance are computed within the red dashed circles shown in Fig. 9a). These recovered values (in blue on the log in Fig. 9b) are in very good agreement with the expected ones (in red in Fig. 9b). On the other hand, the data fit is also very satisfying, as shown in Fig. 10 which compares the observed data with synthetic data computed in the reconstructed model of Fig. 9 (data are shown for the same frequency 8 GHz and source position $θ_s = 270°$ as in Fig. 6).

It is interesting to note in Fig. 9(a) that the image presents a particular geometrical pattern: oscillations with a radial symmetry that manifest themselves as fluctuations on the blue curve in Fig. 9(b). These fluctuations are absent in the background due to the bound constraint $ε_r ≥ 1$ considered in the L-BFGS-B algorithm. In the target, the oscillating pattern is due to the limited wavenumber coverage of our imaging technique, that depends both on the finite-frequency content of the data and on the discrete spatial sampling of the measurements by the acquisition array (Sirgue and Pratt 2004). For comparison, we indicate in Fig. 9(b) low-pass-filtered versions of the expected (red) log, considering the full frequency bandwidth from 2 to 10 GHz (in black), or only the lowest frequency of 2 GHz (in green). For filtering, we use the maximal reconstructed wavenumber given by Sirgue and Pratt (2004) as a function of frequency

<figure><img src="image_15.png" /><figcaption>Figure 10. (a) Amplitude and (b) phase of observed vs. synthetic data computed in the final reconstructed model of Fig. 9, for frequency 8 GHz and the source at $θ_s = 270°$.</figcaption></figure>

and of the minimal illumination angle:

$$k_{max} = \frac{2\omega}{c_o} \cos \left( \frac{|\theta_s - \theta_r|_{min} = 60^{\circ}}{2} \right). \quad (19)$$

As this filter assumes a continuous coverage of the wavenumbers up to $k_{max}$ (i.e. an infinitesimal spatial sampling), the full-bandwidth low-pass-filtered log over-estimates the accuracy of the reconstruction and the actual reconstructed values exhibit a lower wavenumber content: Although very dense, the acquisition sampling still has an imprint in the final image. In particular, the source intervals are responsible for the radial symmetry observed in Fig. 9(a).

The image of Fig. 9(a) can thus be seen as the convolution of the real model with the resolution operator of the imaging technique. Knowing the frequency content of the source and the illumination of the target by the acquisition array, we could deconvolve the reconstructed image of Fig. 9(a) in an impulsive image that would better reflect the reality, as proposed by Ribodetti et al. (2000). This finite-frequency effect also explains why, when the initial model is good enough, inverting all available frequencies simultaneously yields smoother images than inverting them sequentially as proposed by Pratt and Worthington (1990). Using a sequential strategy, the final reconstructed model results from the inversion of the highest frequency and is subject to mono-frequency fluctuations, whereas the simultaneous frequency strategy benefits from a broadband frequency content (Lavoué et al. 2014).
Multiparameter inversion of the dataset
FoamMetExt

To invert the dataset FoamMetExt, we apply the methodology proposed by Lavoué et al. (2014): As the choice of an adequate scaling factor $\beta$ is not straightforward, we first perform several multiparameter inversions independently, using various scaling values, without regularization. Fig. 11 shows the final misfits obtained with the different scaling factors. Based on the data misfit, a scaling factor $\beta = 10$ seems to be the most adequate value but we can observe that misfits of the same order of magnitude can be obtained for lower values, down to $\beta = 0.5$.

As already observed by Lavoué et al. (2014) on synthetic data, the reconstructed models corresponding to roughly equivalent misfits can be quite different. For instance, Figs 12 and 13 show the reconstructed models obtained using scaling factors $\beta = 10$ and $\beta = 0.5$, respectively. The permittivity solutions appear to be particularly sensitive to the choice of the parameter scaling, because the imprint of the foam cylinder on the data is less important than the one of the strongly diffracting copper tube.

The main difference between the two solutions is that, using a scaling factor $\beta = 0.5$, the optimization artificially creates an erroneous permittivity structure inside the copper tube. Based on our forward simulations, we know that this structure is not reliable because the measured data are not sensitive to the inner filling of the highly-conductive copper tube ($\sigma \gg 100$ S/m). Synthetic data, however, may be sensitive to this structure since the reconstructed conductivity value is lower than the actual one ($\sigma \simeq 5$ S/m), in particular in the early iterations of the inversion. The erroneous permittivity reconstruction inside the tube is thus a nice example of the trade-off between permittivity and conductivity: comparing Figs 12 and 13, it can be observed that both permittivity and conductivity reconstructions vary inside the tube, indicating that variations of one parameter compensate the variations of the other regarding the data misfit. These trade-off effects justify the need for considering the Hessian matrix in multiparameter optimization schemes (Operto et al. 2013; Lavoué et al. 2014).

Apart from the artifact inside the copper tube, the two solutions also differ by the amplitude of the fluctuations of permittivity values in the foam cylinder (again, the absence of fluctuations in the background is imposed by the bound constraints $\epsilon_r \ge 1$ and $\sigma \ge 0$ given to the L-BFGS-B algorithm). Finally, the solutions present a thin circular artifact around the copper tube, a region to which the data should be sensitive. These effects can be mitigated by introducing a Tikhonov regularization in the optimization, in order to obtain smoother permittivity models.

Fig. 14 shows the final misfits that we obtain when performing the multiparameter inversion with various scaling factors $\beta$ and regularization weights $\lambda$. As artifacts arise mainly in the permittivity image, whereas conductivity is better constrained by the data, we only apply regularization on the permittivity model (a finer implementation might involve regularization weights adapted to each parameter type, depending on their imprint on the data, but this fine tuning is not critical for our purpose). Contrary to the results obtained by Lavoué et al. (2014), the final data misfits are not more sensitive to the parameter scaling with regula-

rization than without, so the regularization does not enable to distinguish between the solutions obtained with different scaling factors. However, looking at the final reconstructed models, we can observe that the permittivity solutions obtained with regularization are now quite similar (see Figs 15 and 16). The use of regularization erases the permittivity artifacts around the copper tube, yielding satisfactory results. Adequate values for the parameter scaling and for the regularization weight $\lambda$ can be selected based on the data misfit (Fig. 14), given that applying an adequate regularization, we obtain similar results for all scaling factors that display a good data misfit. Note that the Tikhonov regularization does not induce a dramatic smoothing of the solutions, although a multiplicative regularization may better preserve the contrasts in the present case of well-delimited, piecewise-constant targets (Abubakar et al. 2005).

Finally, Fig. 17 compares the observed data with the synthetic data computed in the final reconstructed models of Figs 15 and 16. The fit to the data is naturally dominated by the main peak amplitude of the signal at angles around $180^\circ$, while the data at extreme angles are not well fitted. Applying an angle-dependent weighting to the data through a data covariance matrix in the misfit function could enable to better fit the data presenting minor amplitudes.

CONCLUSIONS

In this study, we have presented a comprehensive view of the inversion of electromagnetic data collected during a well-controlled laboratory experiment. For an accurate and efficient resolution of the forward problem, we use a FDFD scheme where the stencil coefficients are optimized to each simulated frequency. An integral representation of the fields enables to reduce the computation domain to a small zone enclosing the targets and has been shown to be an accurate, efficient, and elegant alternative to brute force calculations in the entire domain. An important ingredient of the modeling consists in accurately positioning virtual sources and receivers on the integration contours using a windowed Sinc interpolation.

Scattered data produced by the purely dielectric target FoamDielExt allowed us to validate our algorithm regarding monoparameter inversion. On the reconstructed permittivity image, the targets are well-delimited and the recovered values are very close to the expected ones. Observed and synthetic data are in very good agreement. We also confirm our strategy for multiparameter inversion on the dataset FoamMetExt. As already observed on synthetic data, the reconstructed models are sensitive to the scaling applied between different parameters types, especially for parameters that are less constrained by the data. These model variations can be barely visible in the data misfit, leading to an ambiguity between the different solutions. An adequate regularization weight enables to mitigate the artifacts, so that satisfactory models can be obtained and, more importantly, identified on the basis of the data misfit analysis. Besides, forward simulations on synthetic models are of great help to estimate the sensitivity of the data to model variations, and hence to evaluate the reliability of the inversion results. Synthetic data are almost insensitive to some parts of the model where our data-driven technique is not able to deduce any valuable
<figure><img src="image_16.png" /><figcaption>Figure 11. Final misfit decrease as a function of the scaling factor $\beta$ used in the multiparameter inversion of the dataset FoamMetExt.</figcaption></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /></figure>

<figure><img src="image_21.png" /></figure>

<figure><img src="image_22.png" /></figure>

<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>Figure 12. (a,c) Permittivity and (b,d) conductivity models obtained by multiparameter inversion of the dataset FoamMetExt, without regularization and using a scaling factor $\beta = 10$. The misfit function has been decreased by 98.2% in 120 iterations. The logs of panels (c) and (d) are extracted along the line $y = 0$ in panels (a) and (b), respectively.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 14. Final data misfits as a function of the scaling factor $\beta$ and of the regularization weight $\lambda$ used in the multiparameter inversion of the dataset FoamMetExt.</figcaption></figure>
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

**Figure 15.** (a,c) Permittivity and (b,d) conductivity models obtained by multiparameter inversion of the dataset FoamMetExt, using a regularization weight $\lambda = 10^{-11}$ and a scaling factor $\beta = 10$. The misfit function has been decreased by 97.1% in 59 iterations. The logs of panels (c) and (d) are extracted along the line $y = 0$ in panels (a) and (b), respectively.

**Figure 16.** (a,c) Permittivity and (b,d) conductivity models obtained by multiparameter inversion of the dataset FoamMetExt, using a regularization weight $\lambda = 10^{-11}$ and a scaling factor $\beta = 0.5$. The misfit function has been decreased by 96.5% in 101 iterations. The logs of panels (c) and (d) are extracted along the line $y = 0$ in panels (a) and (b), respectively.

information without introducing some prior information in the misfit function.

The success of inverting these experimental data lets us envisage the inversion of on-ground GPR field data in a near future. The present work made us free from the inverse crime approach and required to accurately simulate the observed data. When dealing with GPR field data, we expect to face other major obstacles. In particular, we may encounter difficulties for estimating the GPR source signature, and for designing a suitable initial model for starting the full waveform inversion process. The acquisition configuration is also very different in on-ground GPR applications where data are acquired only from the surface, and it has a strong effect on the ability of the imaging technique to recover the subsurface targets.

In particular, in the frame of multiparameter inversion, a partial illumination tends to enhance the trade-off between parameters (Hak and Mulder 2010). Our previous study on synthetic data (Lavoué et al. 2014) suggests that multiparameter FWI can be performed from surface data if the information contained in the Hessian of the misfit function is taken into account through quasi-Newton methods. However, an illustration on real field data still has to be performed. In a preliminary stage, the Institut Fresnel database may constitute an interesting tool to test the sensitivity of the inversion with respect to the acquisition configuration, as well as to investigate optimization methods that better take the Hessian information into account (e.g. the truncated Newton method, Métivier et al. 2013, 2014).

## ACKNOWLEDGMENTS

This study was partially funded by the SEISCOPE consortium (http://seiscope2.osug.fr), sponsored by BP, CGG, Chevron, Exxon-Mobil, JGI, Petrobras, Saudi Aramco,

Schlumberger, Shell, Sinopec, Statoil, Total and Woodside. This work was granted access to the high-performance computing facilities of the CIMENT infrastructure (https://ciment.ujf-grenoble.fr), which is supported by the Rhône-Alpes region (grant CPER07-13 CIRA), the OSUG@2020 labex (reference ANR-10 LABX56) and the Equip@Meso project (reference ANR-10 EQPX-29-01) of the program Investissements d'Avenir supervised by the Agence Nationale pour la Recherche, as well as to the HPC resources of CINES/IDRIS under the allocation 046091 made by GENCI (Grand Equipement National de Calcul Intensif). The authors would like to thank the Institut Fresnel for freely providing their high-quality data (http://iopscience.iop.org/0266-5611/21/6/S09/media). F. Lavoué is particularly grateful to C. Eyraud (Institut Fresnel) for useful informations about the database, to L. Moreau (ISTerre) for his advices on fields integral representation, and to S. Operto (Géoazur, Nice-Sophia Antipolis) for sharing his experience on FDFD schemes and for his careful reading of the manuscript. We also thank two anonymous reviewers for their corrections and X. Dérobert (IFSTTAR, Nantes) for his invitation to publish in this special issue.

## REFERENCES

Abubakar, A., van den Berg, P. M., and Habashy, T. M., 2005. Application of the multiplicative regularized contrast source inversion method on TM- and TE-polarized experimental Fresnel data, *Inverse Problems*, **21**, S5–S13.

Belkebir, K. and Saillard, M., 2001. Special section: Testing inversion algorithms against experimental data, *Inverse Problems*, **17**, 1565–1571.

Belkebir, K. and Saillard, M., 2005. Guest Editors' introduction: Testing inversion algorithms against experimental data: inhomogeneous targets, *Inverse Problems*, **21**, 1–3.
<figure><img src="image_34.png" /><figcaption>**Figure 17.** (a) Amplitude and (b) phase of observed vs. synthetic data computed in the final reconstructed model of Figs 15 and 16, for frequency 18 GHz and the source at $\theta_s = 240°$.</figcaption></figure>

Bérenger, J.-P., 1994. A perfectly matched layer for absorption of electromagnetic waves, *Journal of Computational Physics*, **114**, 185-200.

Busch, S., van der Kruk, J., Bikowski, J., and Vereecken, H., 2012. Quantitative conductivity and permittivity estimation using full-waveform inversion of on-ground GPR data, *Geophysics*, **77**(6), H79-H91.

Byrd, R. H., Lu, P., and Nocedal, J., 1995. A limited memory algorithm for bound constrained optimization, *SIAM Journal on Scientific and Statistical Computing*, **16**, 1190-1208.

Day-Lewis, F. D., Singha, K., and Binley, A. M., 2005. Applying petrophysical models to radar travel time and electrical resistivity tomograms: Resolution-dependent limitations, *Journal of Geophysical Research*, **110**(B8).

Deeds, J. and Bradford, J., 2002. Characterization of an aquitard and direct detection of LNAPL at Hill Air Force base using GPR AVO and migration velocity analyses, in *9th International Conference on Ground Penetrating Radar (GPR 2002)*, Santa Barbara, California (USA), vol. 4758 of SPIE proceedings series, pp. 323-329.

Ellefsen, K. J., Mazzella, A. T., Horton, R. J., and McKenna, J. R., 2011. Phase and amplitude inversion of crosswell radar data, *Geophysics*, **76**(3), J1-J12.

Ernst, J. R., Green, A. G., Maurer, H., and Holliger, K., 2007. Application of a new 2D time-domain full-waveform inversion scheme to crosshole radar data, *Geophysics*, **72**(5), J53-J64.

Garambois, S., Sénéchal, P., and Perroud, H., 2002. On the use of combined geophysical methods to assess water content and water conductivity of near-surface formations, *Journal of Hydrology*, **259**, 32-48.

Geffrin, J.-M., Sabouroux, P., and Eyraud, C., 2005. Free

space experimental scattering database continuation: experimental set-up and measurement precision, *Inverse Problems*, **21**, S117-S130.

Hak, B. and Mulder, W. A., 2010. Migration for velocity and attenuation perturbations, *Geophysical Prospecting*, **58**, 939-951.

Hicks, G. J., 2002. Arbitrary source and receiver positioning in finite-difference schemes using Kaiser windowed sinc functions, *Geophysics*, **67**, 156-166.

Huisman, J. A., Hubbard, S. S., Redman, J. D., and Annan, A. P., 2003. Measuring soil water content with ground penetrating radar: A review, *Vadose Zone Journal*, **2**, 476-491.

Hustedt, B., Operto, S., and Virieux, J., 2004. Mixed-grid and staggered-grid finite difference methods for frequency domain acoustic wave modelling, *Geophysical Journal International*, **157**, 1269-1296.

Ihamouten, A., Villain, G., and Dérobert, X., 2012. Complex permittivity frequency variations from multio{t} offset GPR data: Hydraulic concrete characterization, *IEEE Transactions on Geoscience and Remote Sensing*, **61**(6), 1636-1648.

Jo, C. H., Shin, C., and Suh, J. H., 1996. An optimal 9-point, finite-difference, frequency-space 2D scalar extrapolator, *Geophysics*, **61**, 529-537.

Jol, Harry M., 2009. *Ground Penetrating Radar: Theory and Applications*, Elsevier.

Kalogeropoulos, A., van der Kruk, J., Hugenschmidt, J., Busch, S., and Merz, K., 2011. Chlorides and moisture assessment in concrete by GPR full waveform inversion, *Near Surface Geophysics*, **9**(3), 277-285.

Klotzsche, Anja, van der Kruk, Jan, Linde, Niklas, Doetsch, Joseph, and Vereecken, Harry, 2013. 3-D characterization of high-permeability zones in a gravel aquifer using 2-D crosshole GPR full-waveform inversion and waveguide detection, *Geophysical Journal International*, **195**(2), 932-944.

Lambot, S., Weihermüller, L., Huisman, J. A., Vereecken, H., Vanblooster, M., and Slob, E. C., 2006. Analysis of air-launched ground-penetrating radar techniques to measure the soil surface water content, *Water Resources Research*, **42**, W11403.

Lavoué, F., Brossier, R., Métivier, L., Garambois, S., and Virieux, J., 2013. 2D full waveform inversion of GPR surface data: permittivity and conductivity imaging, in *7th International Workshop on Advanced Ground-Penetrating Radar (IWAGPR 2013)*, Nantes (France).

Lavoué, F., Brossier, R., Métivier, L., Garambois, S., and Virieux, J., 2014. Two-dimensional permittivity and conductivity imaging by full waveform inversion of multio{t} offset GPR data: a frequency-domain quasi-Newton approach, *Geophysical Journal International*, **197**(1), 248-268.

Litman, A. and Crocco, L., 2009. Testing inversion algorithms against experimental data: 3D targets, *Inverse Problems*, **25**(2).

Meles, G. A., van der Kruk, J., Greenhalgh, S. A., Ernst, J. R., Maurer, H., and Green, A. G., 2010. A new vector waveform inversion algorithm for simultaneous updating of conductivity and permittivity parameters from combination crosshole/borehole-to-surface GPR data, *IEEE Transactions on Geoscience and Remote Sensing*, **48**, 3391-3407.

Métivier, L., Brossier, R., Virieux, J., and Operto, S., 2013. Full Waveform Inversion and the truncated Newton method, *SIAM Journal On Scientific Computing*, **35**(2), B401-B437.

Métivier, L., Bretaudieu, F., Brossier, R., Operto, S., and Virieux, J., 2014. Full waveform inversion and the truncated Newton method: quantitative imaging of complex subsurface structures, *Geophysical Prospecting*, **62**, 1353-1375.

Minet, J., Lambot, S., Slob, E. C., and Vanblooster, M., 2010. Soil surface water content estimation by full-waveform GPR signal inversion in the presence of thin layers, *IEEE Transactions on Geoscience and Remote Sensing*, **48**(3), 1138-1150.
MUMPS-team, 2011. *MUMPS - Multifrontal Massively Parallel Solver users' guide - version 4.10.0 (May 10, 2011)*, ENSEEIHT-ENS Lyon, http://mumps.enseeiht.fr or http://graal.ens-lyon.fr/MUMPS, Accessed 15 January 2014.

Operto, S., Brossier, R., Gholami, Y., Métivier, L., Prieux, V., Ribodetti, A., and Virieux, J., 2013. A guided tour of multiparameter full waveform inversion for multicomponent data: from theory to practice, *The Leading Edge, Special section Full Waveform Inversion* (September), 1040–1054.

Patriarca, C., Lambot, S., Mahmoudzadeh, M. R., Minet, J., and Slob, E., 2011. Reconstruction of sub-wavelength fractures and physical properties of masonry media using full-waveform inversion of proximal penetrating radar, *Journal of Applied Geophysics*, **74**, 26–37.

Plessix, R. E., 2006. A review of the adjoint-state method for computing the gradient of a functional with geophysical applications, *Geophysical Journal International*, **167**(2), 495–503.

Pratt, R. G., 1999. Seismic waveform inversion in the frequency domain, part I: theory and verification in a physic scale model, *Geophysics*, **64**, 888–901.

Pratt, R. G. and Worthington, M. H., 1990. Inverse theory applied to multi-source cross-hole tomography. Part I: acoustic wave-equation method, *Geophysical Prospecting*, **38**, 287–310.

Pratt, R. G., Shin, C., and Hicks, G. J., 1998. Gauss-Newton and full Newton methods in frequency-space seismic waveform inversion, *Geophysical Journal International*, **133**, 341–362.

Ribodetti, A., Operto, S., Virieux, J., Lambaré, G., Valéro, H.-P., and Gibert, D., 2000. Asymptotic viscoacoustic diffraction tomography of ultrasonic laboratory data: a tool for rock properties analysis, *Geophysical Journal International*, **140**, 324–340.

Sirgue, L. and Pratt, R. G., 2004. Efficient waveform inversion and imaging: a strategy for selecting temporal frequencies, *Geophysics*, **69**(1), 231–248.

Taflove, A. and Hagness, S. C., 2005. *Computational Electrodynamics: The Finite-Difference Time-Domain Method*, Artech House, 3rd edn.

Tikhonov, A. and Arsenin, V., 1977. *Solution of ill-posed problems*, Winston, Washington, DC.

Velichko, A. and Wilcox, P. D., 2010. A generalized approach for efficient finite element modeling of elastodynamic scattering in two and three dimensions, *Journal of Acoustical Society of America*, **128**(3), 1004–1014.

Weihermüller, L., Huisman, J. A., Lambot, S., Herbst, M., and Vereecken, H., 2007. Mapping the spatial variation of soil water content at the field scale with different ground penetrating radar techniques, *Journal of Hydrology*, **340**, 205–216.

Wilcox, P. D. and Velichko, A., 2010. Efficient frequency-domain finite element modeling of two-dimensional elastodynamic scattering, *Journal of Acoustical Society of America*, **127**(1), 155–165.

Zhao, Y., Xie, X., Wu, J., Chen, J., and Ge, S., 2013. Maxwell curl equation datuming for GPR based on the Kirchhoff integral solution and application in a tunnel grouting test, *Near Surface Geophysics*, **11**, 211–219.

Zhu, C., Byrd, R. H., and Nocedal, J., 1997. Algorithm 778: L-BFGS-B, FORTRAN routines for large scale bound constrained optimization, *ACM Transactions on Mathematical Software*, **23**(4), 550–560.