# Electromagnetic bench based on ETSA antennas for civil engineering materials tomography

Théau Cousin, Antoine Tonnoir, Cyrille Fauchard, Pierre Charbonnier,
Christian Gout

► To cite this version:

Théau Cousin, Antoine Tonnoir, Cyrille Fauchard, Pierre Charbonnier, Christian Gout. Electromagnetic bench based on ETSA antennas for civil engineering materials tomography. The IEEE International Conference on Antenna Measurements and Applications (IEEE CAMA), IEEE Antennas and Propagation Society (IEEE AP-S), Nov 2023, Gènes, Italy. pp.568-573, 10.1109/CAMA57522.2023.10352675. hal-04343862

HAL Id: hal-04343862

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

Submitted on 14 Dec 2023

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_1.png" /></figure>

Distributed under a Creative Commons CC0 1.0 - Universal - International License
# Electromagnetic bench based on ETSA antennas for civil engineering materials tomography

Théau Cousin

LMI - ENDSUM

INSA - Cerema - Routes de France

76000 Rouen, France

theau.cousin@insa-rouen.fr

Antoine Tonnoir

LMI

INSA Rouen

76000 Rouen, France

Cyrille Fauchard

ENDSUM

Cerema

76000 Rouen, France

Pierre Charbonnier

ENDSUM

Cerema

67000 Strasbourg, France

antoine.tonnoir@insa-rouen.fr

cyrille.fauchard@cerema.fr

pierre.charbonnier@cerema.fr

Christian Gout

LMI

INSA Rouen

76000 Rouen, France

christian.gout@insa-rouen.fr

**Abstract**—In this work, we are interested in density reconstruction of road samples, aiming to guarantee road reliability. To do so, a joint work between Routes de France and Cerema has led to an electromagnetic (EM) bench, based on Exponentially Tapered Slot Antennas (ETSA). A mathematical inversion model is used to find the permittivity from antennas parameters, using simulation of electromagnetic waves scattering. From this permittivity, one can recover the density via the Complex Refractive Index Model (CRIM) [1].

**Index Terms**—ETSA antennas, S-parameters, Electromagnetic scattering, Inversion, Permittivity.

## INTRODUCTION

Tomography has been widely used in several contexts. The emergence of new mathematical and computational methods has put electromagnetic waves on show, leading to a lot of applications in medicine [2] or geosciences [3]. Yet, to the best of our knowledge, an electromagnetic bench using low cost equipments to find the permittivity of samples with complex structure (from civil engineering or Earth sciences) combined with recent mathematical approach has been rarely considered.

## I. EQUIPMENTS

The EM bench consists of three motors moving the two antennas and the sample's base, allowing x-y-z axis measurements, see fig. 1. A Vector Network Analyzer (VNA) is used for EM waves emission and reception, connected to two ETSA antennas. They have been developed during I. Boughanmi's thesis [4] according to E. Guillanton's work [5].

## II. MATHEMATICAL BACKGROUND

### A. Forward model

We start from Maxwell's second order equation in time-harmonic domain :

$$ \nabla \times \nabla \times \underline{E} - k_0^2 \varepsilon \underline{E} = 0, \quad (1) $$

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

<figure><img src="image_3.png" /><figcaption>Fig. 1. ETSA Antennas used and EM bench.</figcaption></figure>

where $k_0 = \omega\sqrt{\mu_0\varepsilon_0}$ the propagation factor, $\omega = 2\pi f$ the angular velocity and $\varepsilon$ the relative permittivity, the unknown looked for. Solving the scattering problem amounts to find the scattered field $\underline{E}^s$ such that $\underline{E} = \underline{E}^i + \underline{E}^s$ with $\underline{E}^i$ the incident field, leading to :

$$ \nabla \times \nabla \times \underline{E}^s - k_0^2 \varepsilon \underline{E}^s = k_0^2 (\varepsilon - 1) \underline{E}^i. \quad (2) $$

We also consider Perfectly Matched Layers [6] (PML), to emulate an unbounded domain. In time-harmonic regime, the PML formulation consists in applying a complex stretching $x_i \rightarrow L_{x_i}(x_i)x_i$. This change of variables amounts in replacing in equation (1) the derivative with respect to $x_i$ by $L_{x_i}^{-1}\partial_{x_i}$.

$$ \underline{\nabla} \times \underline{\nabla}_L \times \underline{E} - k_0^2 \varepsilon \underline{L} \underline{E} = k_0^2 (\varepsilon - 1) \underline{L} \underline{E}^i, \quad (3) $$

where

$$ \underline{L} = \begin{pmatrix} L_y L_z & 0 & 0 \\ 0 & L_x L_z & 0 \\ 0 & 0 & L_y L_x \end{pmatrix}, $$
and

$$ \underline{\nabla}_L \times = \begin{pmatrix} 0 & -\frac{L_x}{L_z} \partial_z & \frac{L_x}{L_y} \partial_y \\ \frac{L_y}{L_z} \partial_z & 0 & -\frac{L_y}{L_x} \partial_x \\ -\frac{\tilde{L}_z}{L_y} \partial_y & \frac{L_z}{L_x} \partial_x & 0 \end{pmatrix}. $$

The Finite Element Methods (FEM) discretisation is done by Nedelec elements family [6] associated with the spaces :

$$ \begin{cases} V_h = \{\underline{v} \in H(\mathrm{curl}, \Omega) : v|_K \in N_1(K) \forall \in \Omega\} \\ N_1(K) = [P_0(K)]^3 \oplus S_1(K) \\ S_1(K) = \{\underline{s} \in [P_1(K)]^3 : \underline{s}(\underline{x}) \cdot \underline{x} = 0\}, \end{cases} $$

giving the following discretized forward problem :

$$ [\mathbb{K} - k_0^2 \mathbb{M}(\varepsilon)] \underline{c}(\varepsilon) = \underline{S}(\varepsilon). \quad (4) $$

We suppose here that the forward problem is solved for each value of frequency f from the interval chosen by the user.

## B. S-parameters

The VNA used gives S-parameters, that can be represented as a ratio for each frequencies between In and Out fields (X(f), Y(f)) according to the interested antenna, explained in fig. 2.

<figure><img src="image_4.png" /><figcaption>Fig. 2. S-parameters scheme.</figcaption></figure>

In our case, we consider the 1st antenna as the emitter and the 2nd antenna as the receiver, i.e. the parameter $S_{21}(f)$. For better clarity, we will denote $S_{21}(f)$ by $S$. From the electrical field simulated at the receiving antenna, the associated S-parameter is expressed as follow:

$$ S = \langle \underline{E}, \underline{P} \rangle, $$

where $\underline{P}$ is a vector simulating the rotational aspect of the EM bench.

## C. Inverse problem

The inverse problem consists in finding the permittivity $\varepsilon$ of the sample, given $\underline{S}^{\text{meas}}$. Classically, this leads to the following minimisation problem : Find $\varepsilon$ that minimize

$$ \mathcal{M}(\underline{\varepsilon}^*) = \min_{\underline{\varepsilon}} \frac{1}{2} \|\underline{S}^{\text{meas}} - \underline{\mathcal{F}}(\underline{\varepsilon})\|^2 + \frac{1}{2} \alpha \|\underline{\nabla}\underline{\varepsilon}\|^2, $$

where $\underline{\mathcal{F}}(\underline{\varepsilon})$ corresponds to $\underline{S}^{\text{sim}}$ obtained from the forward problem solved given $\varepsilon$, i.e. :

$$ \underline{\mathcal{F}} : \underline{\varepsilon} \in \mathbb{C}^{n_m} \rightarrow \underline{S}^{\text{sim}} \in \mathbb{C}^{n_r \times n_f}. $$

Let $\underline{\mathcal{L}}(\underline{\varepsilon}) = \underline{S}^{\text{meas}} - \underline{\mathcal{F}}(\underline{\varepsilon})$, the resolution of the minimisation problem is done by a Newton type approach.

* Updating the permittivity $\varepsilon$ at each iteration :

$$ \underline{\varepsilon}^{(k+1)} = \underline{\varepsilon}^{(k)} + \delta \underline{\varepsilon}, $$

* where the variation $\delta \underline{\varepsilon}$ given by :

$$ \underline{\mathcal{J}}^t \underline{\mathcal{J}} \underline{\delta \underline{\varepsilon}} + \underline{\mathcal{J}}^t \underline{\mathcal{r}} = 0. $$

The jacobian $\underline{\mathcal{J}}$ (also called sensitivity matrix) is obtained by the adjoint state method [7] :

$$ \underline{\mathcal{J}}(\underline{\varepsilon}) = \langle \underline{\tilde{c}}(\underline{\varepsilon}), \underline{S}'(\underline{\varepsilon}) - k_0^2 \mathbb{M}'(\underline{\varepsilon}) \underline{c}(\underline{\varepsilon}) \rangle, \quad (5) $$

where $\underline{\tilde{c}}(\underline{\varepsilon})$ is solution to the following adjoint problem :

$$ (\underline{\mathbb{K}} - k_0^2 \underline{\mathbb{M}}(\underline{\varepsilon}))^t \underline{\tilde{c}}(\underline{\varepsilon}) = \underline{r}'_1, \quad (6) $$

where $\underline{r}'_1 = \frac{\partial \underline{r}}{\partial \underline{\varepsilon}}$.

A regularization term $\alpha \|\Delta\underline{\varepsilon}\|$ is used to get a well conditioned problem. For a finite element type mesh (non uniform tetrahedron), the Laplacian can be understood as a finite differences scheme between neighbors tetrahedra, i.e. for a tetrahedron with a permittivity $\varepsilon$ and adjacents tetrahedron with permittivities $\varepsilon_j, \varepsilon_k, \varepsilon_l, \varepsilon_m$ :

$$ \Delta\underline{\varepsilon} = \underline{\varepsilon} - b_j\underline{\varepsilon}_j - b_k\underline{\varepsilon}_k - b_l\underline{\varepsilon}_l - b_m\underline{\varepsilon}_m, \quad (7) $$

where $b_i$ can be either equal to 1 (classic regularization), or $b_i = f(\underline{\varepsilon} - \underline{\varepsilon}_i)$ for $f(x) = c_1(c_2 + x^2)^{-1}$. This technic is called semi-quadratic regularization [8], and allows a non-smooth solution. It is applied after a first convergence of the inversion algorithm, working as a post-processing of the data obtained. The coefficient $\alpha$ is updated every time the norm $\|\underline{S}^{\text{meas}} - \underline{S}^{\text{sim}}\|_{L_2}^2$ reaches a plateau and is still less than a tolerance chosen.

# III. CALIBRATION

## A. Transfer functions

Following S. Lambot idea [9], antennas can be represented as source points emitting and receiving two signals $X_i$ and $Y_i$, associated with transfer functions $H$ (fig. 3). These allow to take into accounts several aspects : incident and transmitted waves from the antenna, reflected waves from the scatterer. It also avoids the antennas modeling. We assume here that every functions depend on the frequency $f$.

<figure><img src="image_5.png" /><figcaption>Fig. 3. Antennas as point sources with their associated transfer functions.</figcaption></figure>
The waves that bounces between the scatterer and the antenna is also taken into account, and can be defined inside the scattering problem (⋆) as shown in fig. 4. In this figure, M depends on the theoretical incident field used. Here, $\underline{E}^i$ is of the form $\nabla \times \Phi$, for $\Phi$ solution of (1) with $\varepsilon_r = 1$ (i.e. $\Phi = H_0(k_0r)$) representing spherical waves modelisation in 2D (for geometrical attenuation aspect). The use of basic transfer functions theory leads to:

$$S = \frac{MH_{r2}H_{t1}}{1 - M(H_{f1} + H_{f2})} \qquad (8)$$

<figure><img src="image_6.png" /><figcaption>Fig. 4. Scattering problem associated to the transfer functions.</figcaption></figure>

## B. Procedure

Transfer functions are estimated from the procedure explained in [9], readapted to our needs. Experimentation has been done to generate a data set of different S value for several distances $d \in \{0, ..., D\}$ and 801 frequencies $f \in [1, 8.5]$ GHz:

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

Using these measurements, and reorganizing equation (8) we obtain the following least square problem:

$$S^d = S^d M^d (H_{f1} + H_{f2}) + M^d H_{r2} H_{t1},$$

where we wish to compute the transfer functions $H_{f1}$, $H_{f2}$, $H_{r2}$ and $H_{t1}$. The problem is of the form $\underline{A} \underline{x} = \underline{b}$ with

$$\begin{array}{ll} * \underline{A} \in \mathbb{C}^{D \times 2}, & A_{d,0} = S^d M^d \\ & A_{d,1} = M^d \end{array}$$

$$\begin{array}{ll} * \underline{x} \in \mathbb{C}^2, & x_0 = H_{f1} + H_{f2} \\ & x_1 = H_{r2} H_{t1} \end{array}$$

$$\begin{array}{ll} * \underline{b} \in \mathbb{C}^D, & b_d = S_{21}^d. \end{array}$$

Solving it for each frequency $f$ eventually gives the needed function transfer vectors, used for taking into accounts the antennas and the bench into our measurement.

## C. Filtering

In equation (8), M is considered as the "perfect" signal. In practice, it is considered as the filtered, calibrated signal of the measured data, $S^c$. Singling M out in (8) leads to:

$$S^c = \frac{S}{S(H_{f1} + H_{f2}) + H_{t1} H_{r2}}$$

The graph of fig. 5 shows the comparison between the S-parameter used for the forward problem and the S-parameter calibrated from the machine in free space conditions, in order to validate the transfer functions.

<figure><img src="image_8.png" /><figcaption>Fig. 5. Real and imaginary part of S-parameter used for simulation and measured after calibration.</figcaption></figure>

# IV. EM BENCH APPLICATION

## A. 2D numerical validation

The data set used is the S-parameter measured at the receiving antenna according to the rotation of the base sample $r$ and frequency $f$, giving a vector of measurement $S^{\text{meas}}$ of size $n_r \times n_f$, indexed by $l = f + rn_f$.

$$\begin{cases} r \in [0, 2\pi] & \text{for } n_r = 8, \\ f \in [2, 4.6] & \text{GHz for } n_f = 20. \end{cases}$$

The inversion process is such that we have a first convergence ($||S^{\text{meas}} - S^{\text{sim}}||^2 < 10^{-3}$) with classic $L_2$-norm for standard regularization, then a second with a semi-quadratic regularization.

The measurements are done from the permittivity shown in fig. 6. The permittivity reconstruction over the inversion process is plotted in fig. 7, with the associated $L_2$-norm error convergence between $S^{\text{meas}}$ and $S^{\text{sim}}$ in fig. 8.

<figure><img src="image_9.png" /><figcaption>Fig. 6. Permittivity map to recover.</figcaption></figure>

Fig. 9 shows the permittivity reconstruction with the semi-quadratic (SQ) regularization, applied after the first convergence, with the coefficient $b_i$ related to the matrix regularization. Fig. 10 shows the $L_2$-norm error convergence associated.

## B. Bench application

Two tests have been conducted to validate the approach. The first one with Fontainebleau sand, whose permittivity is known and allows to validate the whole model, and the second
<figure><img src="image_10.png" /><figcaption>Fig. 7. Evolution of map permittivity over the inversion process.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Fig. 8. $L_2$-norm error convergence between $S^{\text{meas}}$ and $S^{\text{sim}}$ during inversion process. At iteration 15, the algorithm considers it has reached a plateau, changing $\alpha$ and allowing a new convergence.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 9. Evolution of map permittivity over the SQ inversion process, with the coefficients $b_i$ associated to the regularization matrix.</figcaption></figure>

one with a cylinder of road sample buried in sand, as shown in fig. 11.

<figure><img src="image_13.png" /><figcaption>Fig. 10. $L_2$-norm error convergence between $S^{\text{meas}}$ and $S^{\text{sim}}$ for SQ regularization.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. 2D slice of the samples studied</figcaption></figure>

For the first case, Fontainebleau sand has been tested thanks to an EM resonant cavity, considered as a referent method, and gave a permittivity mean of $\bar{\epsilon}_{SDF} = 2.59$. Since the 2D slice is homogeneous, the algorithm converges in a few iterations (fig. 12) and recovers a sample with permittivity $\bar{\epsilon} = 2.57$. Fig. 13 shows the S-parameter from the permittivity recovered and the machine S-parameter.

<figure><img src="image_15.png" /><figcaption>Fig. 12. Evolution of map permittivity for Fontainebleau sand over the inversion process.</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 13. S-parameter measured and calculated in time domain for the Fontainebleau sand.</figcaption></figure>

Second case points out the machine ability to find burried object. Fig. 14 represents the reconstruction of the permittivity, revealing the shape of the road cylinder in the last iterations. The semi-quadratic regularization allows to clearly see the sample's shape, as shown in fig. 15. As for the previous example, fig. 16 plots the S-parameters measured and simulated from the permittivity recovered.

<figure><img src="image_17.png" /><figcaption>Fig. 14. Evolution of map permittivity for the burried sample test over the inversion process.</figcaption></figure>

<figure><img src="image_18.png" /><figcaption>Fig. 15. Map permittivity at the end of the SQ inversion process, with the coefficients $b_i$ associated to the regularization matrix.</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>Fig. 16. S-parameter measured and calculated in time domain for the burried sample test.</figcaption></figure>

## CONCLUSION

### Future work

3D works are in progress, following the domain decomposition method described in [10]. Matrices associated to each subdomains are used as pre-conditioners for a GMRES resolution. First tests have been performed for a homogeneous sample, to validate the algorithm and domain decomposition method. The data set used is a vector of measurement $S^{\text{meas}}$ of size $n_{\text{src}} \times n_{\text{rec}}$, indexed by $l = s+rn_{\text{src}}$. The inversion is done for one frequency without rotation, considering several source point (as shown in fig. 17). Fig. 18 shows the convergence rate between a $S^{\text{meas}}$ from a homogeneous cylinder with $\epsilon = 3$ and $S^{\text{sim}}$ simulated over the iterations.

<figure><img src="image_20.png" /><figcaption>Fig. 17. Instance of different sources positions.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Fig. 18. $L_2$-norm error convergence between $S^{\text{meas}}$ and $S^{\text{sim}}$.</figcaption></figure>
Last work are waiting for samples with known density to study them and validate the whole density measurement procedure since the main objective is to replace the standard nuclear method used in laboratory to control this parameter.

## Assessment

The modelisation of the bench and antennas as source points, the implementation of the inverse problem resolution and a semi-quadratic regularization method have been set up thanks to this work. The 2D cases have been tested and validated, and 3D cases are still in development.

## ACKNOWLEDGMENT

The author would like to thank Brice Delaporte, Hervé Dumont and Christine Leroy from Routes de France and the Cerema for their help and support in making this work possible. Future collaborations are being discussed for national projects (Dephy3Geo, ...) and geophysical applications of the machine.

## REFERENCES

[1] C. Fauchard, B. Li, L. Laguerre, B. Héritier, N. Benjelloun, M. Kadi, "Determination of the compaction of hot mix asphalt using high-frequency electromagnetic methods", NDT & E International, 60, pp.40-51, 2013.

[2] V.-L. Coli, P.-H. Tournier, V. Dolean, I. El Kanfoud, C. Pichot, "Detection of Simulated Brain Strokes Using Microwave Tomography", IEEE journal of electromagnetics, RF and microwaves in medicine and biology, 2019.

[3] S. Huang, J. Lin, Q. Huang and R. Liang, "An Emerging Method Using EM Wave Computed Tomography for the Detection of Karst Caves", Geotechnical and Geological Engineering, 38, pp. 2713-2723, 2020.

[4] I. Boughanmni, "Contribution au développement, à la validation et à l'expérimentation d'un banc de caractérisation diélectrique des matériaux en tomographie EM", Thèse de Doctorat, Université de Rouen Normandie, pp. 175, 2021.

[5] E. Guillanton, "Etude d'un système d'imagerie microonde multistatique- multifréquence pour la reconstruction d'objets enfouis", Thèse de Doctorat, Laboratoire d'Electronique, Antennes et Télécommunications, 2000.

[6] P. Monk, "Analysis of a Finite Element Methods for Maxwell's Equations", SIAM on Numer. Analysis 29 (3), pp. 714-729, 2005.

[7] D. Givoli, "A tutorial on the adjoint method for inverse problems", Computer Methods in Applied Mechanics and Engineering, vol. 380, pp. 113-810, 2021.

[8] P. Charbonnier, L. Blanc-Féraud, G. Aubert, M. Barlaud, "Deterministic Edge-Preserving Regularization in Computed Imaging", IEEE Transactions on Image Processing, 6, pp. 298-311, 1997.

[9] S. Lambot, "Modeling of ground-penetrating Radar for accurate characterization of subsurface electric properties", IEEE Trans. on Geoscience 42, pp. 2555-2568, 2004.

[10] M. El Bouajaji, B. Thierry, X. Antoine and C. Geuzaine, "A Quasi-Optimal Domain Decomposition Algorithm for the Time-Harmonic Maxwell's Equations", Journal of Computational Physics, 2015.