Erreur de contraste normalisée en fonction de N s et N r . (a) CSI et (b) sont pour un CSI basé sur des ondelettes db14. Un bruit gaussien de 20 dB est ajouté aux champs diffractés. . . . xiii 3 (a) : Vérité terrain. Permittivité reconstruite par (b) CSI, (c) GS-CSI avec le premier arrangement de groupement, et (d) CSI ℓ 1 -pénalisée. N s × N r = 16 × 24. Un bruit gaussien de 20 dB est ajouté aux champs diffractés. . . . . .

4.9 Parametric study for the choice of q ℜ and q ℑ . The final error Err (χ) is plotted as a function of q ℜ and q ℑ and compared to the one of CSI (green line). N s × N r = 32 × 32 . . . . . . . . 4.10 "Austria" model -Relative permittivity ϵ(r) and conductivity σ(r) (mS/m) . . . . . . . . . . . 4.11 Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 20 dB is added to the scattered fields.

4.12 Inversion error curves of (a) CSI and (b)-(d) CSI-GS with three grouping arrangements. Gaussian noise of 20 dB is added to the scattered fields. . . . . . . . . . . . . . . . . . . . . . . . 4.13 Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 15 dB is added to the scattered fields. 4.14 Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 10 dB is added to the scattered fields. 4.15 Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 20 dB is added to the scattered fields. [1.3, 2]. Gaussian noise of SNR = 20 dB is added. . . . . . . .

6.4

Relative permittivity of reconstructed profiles of samples #4-#6 from the test dataset. The range of relative permittivity is [1.3, 2]. Gaussian noise of SNR = 20 dB is added. . . . . . . . 6.5 Err χ , Err

Résumé en français

Le contexte de recherche

Cette thèse est une contribution aux problèmes de diffraction électromagnétique inverse (dorénavant notés ISP comme Inverse Scattering Problems). Les problèmes de diffraction directs consistent à déterminer la distribution des champs qui sont présents dans une zone donnée, à partir des informations sur les ondes incidentes sur cette zone, et sur le ou les diffracteurs en son sein. De manière correspondante, comme de fait l'inverse du problème direct, la résolution des ISP vise à retrouver les caractéristiques d'un ou de diffracteurs inconnus, telles que la géométrie ou la distribution des paramètres physiques, à partir de la connaissance des champs qu'il(ils) diffracte(nt) en présence de sources de rayonnement électromagnétique données [START_REF] Chen | Computational Methods for Electromagnetic Inverse Scattering[END_REF].

L'étude des ISP est une préoccupation constante car ceux-ci sont pertinents sous une forme ou sous une autre dans un large éventail d'applications [START_REF] Pastorino | Microwave Imaging Methods and Applications[END_REF], telles, et sans exhaustivité, le contrôle non destructif (c'est-à-dire la détection des défauts d'un matériau, d'un composant ou d'un système sans y causer de dommages) [START_REF] Ida | Handbook of Advanced Nondestructive Evaluation[END_REF][START_REF] Zoughi | Microwave Non-destructive Testing and Evaluation[END_REF], l'imagerie médicale (par exemple, la détection des tumeurs du sein [START_REF] Qin | Early Breast Anomalies Detection with Microwave and Ultrasound Modalities[END_REF][START_REF] Nikolova | Microwave Imaging for Breast Cancer[END_REF])), la télédétection [START_REF] Varotsos | Microwave Remote Sensing Tools in Environmental Science[END_REF][START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF], et la caractérisation de matériaux [START_REF] Alexopoulos | Effective Parameters for Metamorphic Materials and Metamaterials through a Resonant Inverse Scattering Approach[END_REF].

Les problèmes directs de diffraction qui consistent donc à calculer les champs diffractés à partir de la résolution des équations pertinentes qu'ils satisfont sont linéaires et bien posés au sens de Hadamard, par exemple, une solution unique existe et elle dépend continûment des données. Au contraire, les problèmes inverses sont généralement mal posés. En particulier, contrairement aux autres problèmes inverses linéaires, tels la tomodensitométrie [START_REF] Buzug | Computed Tomography[END_REF]) ou la déconvolution [START_REF] Campisi | Blind Image Deconvolution: Theory and Applications[END_REF], les ISP sont non linéaires en raison des phénomènes de diffraction multiple.

Par conséquent, le défi de la mise en oeuvre d'algorithmes fiables et efficaces pour les ISP est double.

D'une part, les ISP souffrent d'un caractère mal posé, ce qui se traduit par des solutions qui ne sont pas uniques et sont très sensibles à toute variation de données ou de paramètres. D'autre part, le modèle inverse étant non linéaire, la fonction d'inadéquation aux données correspondante pour l'inversion est susceptible d'être non convexe, ce qui rend l'inversion plus difficile voire inefficace.

Objectif et organisation de la thèse

Cette thèse se concentre sur l'exploration des ISP bidimensionnels full-wave, où par full-wave on entend qu'aucune approximation de l'interaction électromagnétique n'est effectuée dans le modèle mis en place. Une façon courante de s'attaquer aux ISP consiste à utiliser des approches déterministes, avec certaines informations préalables incorporées dans la fonction de coût à minimiser. La parcimonie est une information a priori largement utilisée à cet effet et la transformée discrète en ondelettes (Discrete Wavelet Transform ou DWT) [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]) est souvent employée pour atteindre une représentation parcimonieuse de l'inconnu. Dans cette thèse, la première partie porte ainsi sur l'utilisation des ondelettes et la régularisation de la parcimonie dans les ISP, basée sur la méthode CSI mentionnée ci-dessus. Une approche multi-échelle basée sur les ondelettes est d'abord explorée. L'intégration de la parcimonie peut être réalisée en ajoutant des termes de régularisation, qui sont généralement sous la forme de normes vectorielles. Une approche pénalisée par la parcimonie par groupe est donc proposée, en ajoutant une norme mixte à la fonction de coût de la méthode CSI. Ensuite, la seconde partie concerne l'application de techniques d'apprentissage en profondeur pour s'attaquer aux ISP. La combinaison de la méthode itérative traditionnelle et du CNN est réalisée par une approche dite déroulée ou unrolled.

La thèse consiste en sept chapitres et est structurée comme suit :

• Le chapitre 1 décrit le contexte du travail et résume la contribution de la thèse. Une introduction générale aux problèmes inverses, aux algorithmes d'inversion pour les ISP et aux difficultés rencontrées est proposée.

• Le chapitre 2 présente la formulation du problème direct. La méthode des moments utilisée pour la modélisation est introduite.

• Le chapitre 3 illustre l'algorithme CSI traditionnel. La décomposition en ondelettes est brièvement introduite, et l'approche CSI basée sur les ondelettes est décrite.

• La transformée en ondelettes complexe à double arborescence est introduite au chapitre 4, et une méthode de solution de source de contraste pénalisée par la parcimonie par groupe est investiguée.

• Au chapitre 5, les détails du réseau de neurones convolutifs (CNN) sont introduits et la méthode déroulée est discutée.

• Au chapitre 6, un réseau de neurones convolutionnels déroulés dit CSI-Net est proposé, combiné avec la méthode CSI. Une comparaison avec l'état de l'art est effectuée. Ensuite, le CSI-Net est étendu au schéma multi-fréquence.

• Le chapitre 7 résume le travail accompli et introduit quelques perspectives.

• La dissémination associée au présent travail est précisée à la fin du manuscrit.

Problème direct de diffraction électromagnétique ( §2)

E tot (r, r s ) = E inc (r, r s ) + D G(r, r ′ )J(r ′ , r s )dr ′ , r ∈ D (1) et E diff (r r , r s ) = D G(r r , r ′ )J(r ′ , r s )dr ′ , r ∈ D, (2)
où J(r, r s ) est la source de contraste J(r, r s ) = χ(r)E tot (r, r s ). E tot sont respectivement les champs totaux, E diff les champs diffractés, et E inc les champs incidents. G représente la fonction de Green pertinente. Le problème direct consiste donc à calculer les champs diffractés en résolvant les équations ci-dessus, ce qui est un problème bien posé. Le problème inverse correspondant consiste à récupérer les informations physiques du diffracteur, e.g., la fonction de contraste χ(r), à partir des données collectées E diff et des informations connues E inc et G.

D C Sources

Receivers

Ensuite, une méthode de moments (MoM, pour Method of Moments) est utilisée afin de transformer les équations intégrales vectorielles continues en équations algébriques scalaires discrètes :

E tot i = E inc i + G D J i , E diff i = G S J i et J i = diag(χ)E tot i , où i = 1, . . . , N s avec N s le nombre de sources, G D et G S étant
les fonctions de Green.

Méthode d'inversion de source de contraste et technique d'ondelettes ( §3)

Ce chapitre présente une méthode d'inversion de source de contraste multi-échelle basée sur la transformée discrète en ondelettes (DWT) et la méthode traditionnelle d'inversion de source de contraste (CSI).

Cette dernière consiste à minimiser une fonctionnelle de coût de la forme suivante

F (J i , χ) = Ns i=1 ∥ξ i -G S J i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + Ns i=1 ∥diag(χ)E inc i -J i + diag(χ)G D J i ∥ 2 Ns i=1 ∥χE inc i ∥ 2
(3) Trois arrangements de groupements sont proposés pour diviser les coefficients d'ondelettes en sousgroupes qui se chevauchent. La stratégie de réplication est utilisée dans le but de découpler les coefficients qui se chevauchent, de sorte que le problème d'optimisation pénalisé par la norme ℓ 2,1 peut être résolu par la méthode de descente de gradient proximal avec un opérateur proximal dit de forme fermée. Des simulations exhaustives sont réalisées pour évaluer les performances de l'approche proposée. La figure 3 illustre les résultats obtenus par CSI, CSI pénalisé par la parcimonie par groupe (noté GS-CSI) et CSI pénalisé par la norme ℓ 1 pour préciser et de fait confirmer l'avantage de la méthode proposée.

ξ i étant
Bases du réseau de neurones convolutifs (CNN) ( §5) En ce qui concerne le CSI-Net, la taille de pas de l'algorithme CSI n'est pas un paramètre apprenable. Faire de la taille du pas un paramètre apprenable est une voie à suivre, qui permet de simplifier la structure du réseau. Les informations physiques peuvent être intégrées dans le modèle non seulement via des entrées et des structures de réseau pré-initialisées, mais également via des modifications de la fonction de perte (loss fnction). La fonction de perte utilisée dans le travail actuel peut être encore modifiée en ajoutant des termes d'inadéquation aux données obtenus à partir de la connaissance du domaine, comme proposé dans [START_REF] Liu | Physics-guided Loss Functions Improve Deep Learning Performance in Inverse Scattering[END_REF].

De plus, toutes les méthodes proposées sont encore limitées aux situations 2-D et alors au cas de polarisation TM standard. Plus d'investigations doivent être menées pour les cas pleinement 3-D pour lesquels la modélisation mathématique est plus compliquée, et la non-linéarité et le caractère mal posé deviennent plus sévères. De plus, les méthodes proposées ne sont testées à ce stade que sur des données de simulation. D'autres tests devraient être effectués sur des données expérimentales, au départ en situation contrôlée de laboratoire.

Introduction

Research background

This thesis is a contribution to the electromagnetic inverse scattering problems (ISPs). The direct scattering problems are about determining the distribution of the fields present in a given area, based on information from incident waves and from the scatterer. Correspondingly, as the inverse to the direct problem, solving the ISPs aims at retrieving the characteristics of an unknown scatterer, such as its geometry or the distribution of its physical parameters, from the knowledge of the fields it scatters when it is illuminated by given electromagnetic sources [START_REF] Chen | Computational Methods for Electromagnetic Inverse Scattering[END_REF]. Investigating the ISPs is a continuing concern since they are relevant to a wide range of applications [START_REF] Pastorino | Microwave Imaging Methods and Applications[END_REF], including non-destructive testing (i.e., detecting defects of a material, component or system without causing damage) [START_REF] Ida | Handbook of Advanced Nondestructive Evaluation[END_REF][START_REF] Zoughi | Microwave Non-destructive Testing and Evaluation[END_REF], medical imaging (e.g., detecting breast tumors [START_REF] Qin | Early Breast Anomalies Detection with Microwave and Ultrasound Modalities[END_REF][START_REF] Nikolova | Microwave Imaging for Breast Cancer[END_REF], remote sensing [START_REF] Varotsos | Microwave Remote Sensing Tools in Environmental Science[END_REF][START_REF] Tsang | Theory of Microwave Remote Sensing[END_REF], and material characterization [START_REF] Alexopoulos | Effective Parameters for Metamorphic Materials and Metamaterials through a Resonant Inverse Scattering Approach[END_REF]).

The direct scattering problems which consist in calculating the scattered fields based on the Maxwell equations are linear and well-posed in Hadamard's sense, e.g.,. a unique solution exists and it depends continuously on its data. On the contrary, inverse problems are generally ill-posed. Especially, different from other linear inverse problems, i.e., computed tomography [START_REF] Buzug | Computed Tomography[END_REF], deconvolution [START_REF] Campisi | Blind Image Deconvolution: Theory and Applications[END_REF], the ISPs are non-linear due to multiple scattering. Hence, the challenge of implementing reliable and efficient algorithms for ISPs is twofold. On the one hand, ISPs suffer from ill-posedness, which results in non-unique and highly sensitive solutions. On the other hand, the inversion model being nonlinear, the corresponding data misfit function for inversion is likely to be non-convex, making the inversion more difficult and inefficient.

Overview of inversion techniques

Inverse problems are essentially optimization problems, which aim at minimizing the discrepancy between estimated results and expected ones

x = arg min x F (O, K(x)) (1.1)
K is the forward map (linear or nonlinear), O the measurements. x means the input elements to be optimized, F is the data misfit function which is usually modeled by the L2 norm (Euclidean norm). Iterative methods are most often employed to optimize the cost function and find the solution, and they are generally categorized into two types.

The deterministic methods consist in minimizing the cost function and searching for the optimal solution step by step according to the search direction and step size which are determined at each iteration. Deterministic methods are basically local optimization methods, which are efficient in finding local minima for non-convex optimization problems. Given an initial guess, the deterministic methods are capable of searching the local minimum that is closest to the initial point. However, it is difficult for them to find the global optimal solution unless provided with a reasonable initial guess which is close to the global minimum. Moreover, deterministic methods can realize the local linearization of a nonlinear problem at each iterative step. There are some commonly used deterministic methods, e.g., steepest descent method [START_REF] Curry | The Method of Steepest Descent for Non-linear Minimization Problems[END_REF], conjugate gradient descent method [START_REF] Dai | A Nonlinear Conjugate Gradient Method with a Strong Global Convergence Property[END_REF], and Gauss-Newton algorithm [START_REF] Hartley | The Modified Gauss-Newton Method for the Fitting of Non-linear Regression Functions by Least Squares[END_REF].

The nonlinear inverse problems can be tackled by linear or quasi-linear algorithms with a reasonable initial guess provided. Yet, such methods usually fail when the problem has a strong nonlinearity which increases the non-smoothness of the cost function, and as a result, a more complicated relationship between data and solution.

Another type of iterative optimization methods is called stochastic methods, e.g., colony optimization [START_REF] Pastorino | Stochastic Optimization Methods Applied to Microwave Imaging: A Review[END_REF], genetic algorithms [START_REF] Qing | Electromagnetic Inverse Scattering of Two-dimensional Perfectly Conducting Objects by Real-coded Genetic Algorithm[END_REF], which have been proposed to remedy the drawback of local optimization methods. They have a better capability of finding the global minimum however with no global convergence guarantee, with a low convergence rate due to a larger computational burden compared to deterministic methods. In this thesis, the stochastic methods are not envisaged.

A central feature of inverse problems is ill-posedness. It can be interpreted as the lack of information, so that the solution either does not exist or cannot be uniquely determined. For instance, the ill-posedness of ISPs is caused by the presence of non-radiating currents which do not contribute to the scattered fields. To alleviate the ill-posedness of inverse problems, prior information should be added. From a mathematical point of view, regularization terms should be added to the original cost function, which allows the prior information to be integrated into the inversion process.

x = arg min x F(O, K(x)) + ΛG(x) (1.2)
G(•) is the regularization term, and Λ is the corresponding hyperparameter which is related to the trade-off between the data fitting term and regularization terms. Denote x = K -1 (y) with K -1 the inverse operator of K, there are two errors that are produced by the cost function (1.2)

e 1 = K-1 (y) -x (1.3) e 2 = K-1 (O) -K-1 (y) (1.4)
K-1 is the estimation of the inverse operator. e 1 , which is brought by the regularization term, measures the discrepancy between the estimated inverse operator and the real one. e 2 measures the misfit of observation and estimated data, which is represented by the data fitting term. However, e 1 and e 2 cannot both reach a small value. Hence, the hyperparameter Λ controls the trade-off.

To quickly summarize, ill-posed inverse problems usually call for additional regularization terms to stabilize the inversion process. Various regularization schemes have been proposed (Engl et al. 2000) to insert the prior information or constraints, and a simple example is the positivity constraint. In addition, the Tikhonov regularization [START_REF] Tikhonov | Numerical Methods for the Solution of Ill-posed Problems[END_REF] is one of the most commonly used techniques. Note that the regularization terms are also expressed by the ℓ p norm of a vector/ matrix, which can be defined as

∥x∥ p = i |x i | p 1 p (1.5)
A general expression of the Tikhonov regularization term is G(x) = ∥x∥ 2 2 which can penalize the entries with large norm. ∥ • ∥ 2 is the ℓ 2 norm so it is also called ℓ 2 regularization. A major drawback of the Tikhonov regularization is the over-smoothness of the solution, i.e., the loss of non-smooth information.

For example, the use of Tikhonov regularization in image reconstruction problems can lead to the loss of detail information in the image. [START_REF] Rudin | Nonlinear Total Variation Based Noise Removal Algorithms[END_REF]) introduces the total variation (TV) regularization, which is a well-known technique in noise removal. For a one dimensional digital signal x, the TV regularization has the following form: G(x) = i |x i+1 -x i |. It can reduce noise while keeping the detail information, e.g., the edges present in a 2-D image.

Sparsity [START_REF] Bach | Optimization with Sparsity-inducing Penalties[END_REF]) is also a widely used prior information that helps to improve the inversion results. A signal can be considered as sparse if most of its elements are zero, or almost zero. The simplest approach to determining whether a signal is sparse or not is to count its non-zero elements, which is exactly the definition of ℓ 0 norm ∥x∥ 0 = #{i = 1, . . . , N |x i ̸ = 0}

(1.6)

By adding the ℓ 0 norm [START_REF] Horn | Norms for Vectors and Matrices[END_REF] as the regularization term, the number of nonzero elements contained in the data is penalized. Hence, sparsity can be enforced on the data. However, the optimization of a non-convex cost function with ℓ 0 norm is a NP-hard problem (Van Leeuwen 1991), which can be as difficult as trying out all possible non-zero supports of coefficients. A common solution is to replace the ℓ 0 norm with other indirect norms which can also promote the sparsity. A natural approximation of the ℓ 0 norm is the ℓ 1 norm, the solution of which has been verified to be very close to the sparsest solution [START_REF] Donoho | For Most Large Underdetermined Systems of Linear Equations the Minimal ℓ 1 -norm Solution Is Also the Sparsest Solution[END_REF].

The optimization of a cost function with ℓ 1 norm embedded is also called the least absolute shrinkage and selection operator (LASSO) [START_REF] Tibshirani | Regression Shrinkage and Selection Via the Lasso[END_REF]. Furthermore, some mixed norms [START_REF] Kowalski | Sparse Regression Using Mixed Norms[END_REF]), e.g. ℓ 2,1 , ℓ ∞,1 , can also promote the sparsity while preserving structural information of the signal.

In addition to iterative optimization approaches, deep-learning-based approaches have also been intensively investigated in recent years. Deep learning, mostly convolutional neural networks (CNN), has been widely applied in inverse imaging problems [START_REF] Lucas | Using Deep Neural Networks for Inverse Problems in Imaging: Beyond Analytical Methods[END_REF][START_REF] Jin | Deep Convolutional Neural Network for Inverse Problems in Imaging[END_REF], such as image denoising [START_REF] Zhang | Image Denoising Using a Neural Network Based Non-linear Filter in Wavelet Domain[END_REF], restoration [START_REF] Zhang | Learning Deep CNN Denoiser Prior for Image Restoration[END_REF], and super-resolution [START_REF] Cui | Deep Network Cascade for Image Super-resolution[END_REF]. Depending on whether the forward model F is known or not, different deep learning techniques can be applied. In our work, the forward model is fully known and one has full access to the paired ground truth and measurements. Hence, the simplest way to apply deep learning techniques is to construct a network that maps the measurements to the ground truth. As the forward model is known, some domain knowledge can be further incorporated into the network if provided with a reasonable initialization. Moreover, in addition to applying a single network to map the input to the output, the unrolled method [START_REF] Ongie | Deep Learning Techniques for Inverse Problems in Imaging[END_REF]) has been proposed to tackle inverse problems, especially those for which the forward model is fully known.

This approach is capable of integrating traditional iterative solvers into the CNN structure, which provides a fertile ground for future study.

The overview of the application of the above techniques in ISPs is provided in the following chapters.

Inverse scattering methods

As previously classified, there are three categories of methods for solving inverse problems: deterministic methods, stochastic methods, and deep-learning-assisted approaches.

For inverse scattering problems (ISPs), there are non-iterative deterministic methods like Rytov [START_REF] Devaney | Inverse-scattering Theory within the Rytov Approximation[END_REF] and Born [START_REF] Devaney | Inversion Formula for Inverse Scattering within the Born Approximation[END_REF] approximations, which employ first-order approximations to linearize the problem. This type of methods requires less computational time, but they are mostly used for weak-enough scatterers. They can also be used as the initial point for iterative solvers.

To more rigorously deal with high contrast scatterers, a number of iterative methods with many variants have been developed. In the distorted Born iterative method (DBIM) [START_REF] Chew | Reconstruction of Two-dimensional Permittivity Distribution Using the Distorted Born Iterative Method[END_REF]), the scatterer is regarded as a perturbation to the background medium. At each iteration, the incident fields and the Green functions are updated and the scattered fields generated by the perturbation part are minimized. Here, the nonlinear inverse problem is linearized at each iteration. The Contrast Source Inversion method (CSI) (Van den [START_REF] Van Den Berg | Contrast Source Inversion Method: State of Art[END_REF] consists in minimizing the cost functional built with source-type integral equations.

At each iteration, the contrast source and the contrast are alternatively updated, and the solution of the forward problem is avoided. The subspace-based optimization method (SOM) [START_REF] Chen | Subspace-based Optimization Method for Solving Inverse-scattering Problems[END_REF]) is also a widely used method that treats the source-type integral equations. It consists in determining the main part of the contrast source by spectral analysis, and then obtaining the remaining part by an optimization method. The contrast source is thus divided into two parts: signal subspace and noise subspace. Hence, this method enjoys good robustness against noise. Compared with first-order-approximation-based methods, iterative methods suffer from a higher computational burden yet they have a much wider range of applications.

Stochastic approaches, e.g., colony optimization method [START_REF] Pastorino | Stochastic Optimization Methods Applied to Microwave Imaging: A Review[END_REF]) and genetic algorithms [START_REF] Qing | Electromagnetic Inverse Scattering of Two-dimensional Perfectly Conducting Objects by Real-coded Genetic Algorithm[END_REF], have a better capability to find global optimal solutions, yet they are generally less popular in the community due to a large computational burden.

Deep-learning-assisted approaches for solving ISPs have been much investigated. A simple way to apply deep learning techniques is to directly regress the contrast distribution from the scattered fields, which is the input of the neural network, the output being the contrast. In [START_REF] Ran | Subwavelength Microstructure Probing by Binary-Specialized Methods: Contrast Source and Convolutional Neural Networks[END_REF], the CNN is used to diagnose the damaged microstructure and the designed architecture takes the collected data as the input and outputs the contrast distribution.

Nevertheless, learning directly the contrast from the scattered fields requires unnecessary time to learn the known physical process. As a result, the combination of domain knowledge and neural networks has been much exploited. In (Wei et al. 2019a) and (Wei et al. 2019b), several physical-based learning approaches have been proposed such as the back-propagation scheme (BPS) and the induced current learning method (ICLM), which pre-initialize the input of the neural network based on the domain knowledge. In [START_REF] Liu | Physics-guided Loss Functions Improve Deep Learning Performance in Inverse Scattering[END_REF], two physical-guided loss functions are proposed, which on the one hand enhance the noise robustness, and on the other hand, improve the reconstruction accuracy. These methods can reduce the nonlinearity of the relationship between the input and output of the neural network and simplify the learning process.

The neural network can be combined also with traditional optimization approaches. In [START_REF] Sanghvi | Embedding Deep Learning in Inverse Scattering Problems[END_REF], the learning process is followed by a refinement stage, where the SOM is applied in order to refine the result obtained from the previous stage. Apart from the CNN, the Generative Adversarial Network (GAN) [START_REF] Goodfellow | Generative Adversarial Nets[END_REF] and Recurrent Neural Network (RNN) have been applied also for tackling ISPs [START_REF] Ran | Imaging of Subwavelength Microstructures by Time Reversal and Neural Networks, from Synthetic to Laboratory-controlled Data[END_REF][START_REF] Ye | An Inhomogeneous Background Imaging Method Based on Generative Adversarial Network[END_REF][START_REF] Song | Electromagnetic Inverse Scattering with Perceptual Generative Adversarial Networks[END_REF]), yet one will not dwell upon those here as out of our scope.

Objective and outline

As introduced above, a common way to tackle the ISPs is to employ deterministic approaches, with some prior information incorporated in the cost function to minimize. Sparsity is a widely used prior information and the discrete wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]) is often employed to reach a sparse representation of the unknown. In this thesis, the first part is about the use of sparsity in ISPs, based on the CSI method mentioned above. A wavelet-based multiscale approach is first explored. As previously mentioned, the integration of sparsity can be realized by adding regularization terms, which are usually of the form of vector norms. A group sparsity penalized approach is hence proposed, by adding a mixed norm to the cost function of CSI method. Then, the second part is about applying deep learning techniques to tackle the ISPs. The combination of traditional iterative methods and CNN is realized by the so-called unrolled approach.

The thesis contains eight main chapters in full and is structured as follows:

• Chapter 1 (here) describes the background of the work and summarizes the contribution of this thesis. An overall introduction of inverse problems, inversion algorithms for ISPs, and difficulties faced with is given.

• Chapter 2 presents the formulation of the forward problem. The method of moments that is used for modeling is briefly introduced.

• Chapter 3 illustrates the traditional CSI algorithm. The wavelet decomposition is briefly introduced, and the wavelet-based CSI approach is proposed.

• The dual-tree complex wavelet transform is introduced in Chapter 4, and a group sparsity penalized contrast source solution method is proposed.

• In Chapter 5, the details of the convolutional neural network (CNN) are introduced and the unrolled method is presented.

• In Chapter 6, an unrolled convolutional neural network CSI-Net is proposed, combined with the CSI method. A comparison with the state of the art is proposed. Then, the CSI-Net is extended to the multi-frequency scheme.

• Chapter 7 summarizes the achieved work and provides some perspectives.

• The dissemination associated with the present work is at the end of the manuscript.

Forward problem

In this chapter, the forward model is firstly built, from which the inverse model can be derived. The direct electromagnetic scattering problem is about modeling the physical interaction between the object that is supposed to be well-known, and the collected data, i.e., scattered fields.

Formulation

Consider a two-dimensional non-magnetic homogeneous isotropic object, immersed in a homogeneous medium contained in a given test domain. Assume that the object is infinitely long and uniformly distributed along the axis z. The incident fields are uniformly distributed along the axis z, and the polarization direction is parallel to the axis. The object is illuminated by time-harmonic electromagnetic waves E inc (r) with time dependency e -iωt , where ω is the angular frequency. Due to the phenomena of electromagnetic induction and polarization, the object can be regarded as a secondary source that emits secondary fields in the whole space, i.e., the scattered fields E diff (r). r ∈ D is a given observation point.

Maxwell's equations describe the properties of electromagnetic fields and the relationships between each component [START_REF] Tsang | Scattering of Electromagnetic Waves: Theories and Applications[END_REF]. Their formulation in isotropic and homogeneous medium without source or charge reads as

∇ × E(r) = iωµH(r) (2.1) ∇ × H(r) = J (r) -iωε(r)E(r) = [σ(r) -iωε(r)] E(r) (2.2) ∇ • E = 0 (2.3) ∇ • H = 0 (2.4)
where E, H and J represent the electric and magnetic field strength, and the electric current density, respectively. µ = µ 0 µ r is the permeability with µ 0 = 4π × 10 -7 Henry/meter the permeability in free space, µ r = 1. ε(r) = ε 0 ε r is the permittivity with ε 0 = 8.854 × 10 -12 Farad/meter the permittivity in free space and ε r the relative permittivity of the scatterer. σ(r) is the conductivity (σ = 0 Siemens/meter in free space).

By eliminating H in equations (2.1), (2.2), we have

∇ × ∇ × E(r) = ε c (r)k 2 B E(r) (2.5) where ε c (r) = ε r (r) + iσ(r) ωε 0 is the complex relative permittivity, k 2 B = ω 2 ε 0 µ 0 is the propagation constant
of the embedding medium.

The incident fields E inc (r) satisfy the scalar Helmholtz equation

∇ 2 + k 2 B E(r) = 0 (2.6)
and the total fields are the superposition of incident fields and scattered fields

E tot (r) = E inc (r) + E diff (r) (2.7)
From equations (2.5) and ( 2.6), we can get

∇ × ∇ × E diff (r) -k 2 B E diff (r) = k 2 B [ε c (r) -1]E tot (r) (2.8) wherein χ(r) = k 2 B [ε c (r) -1] = k 2 (r) -k 2 B (2.9)
is defined as the contrast function. Also, the scattered fields satisfy the Sommerfield radiation condition, so we have

E diff (r) = D G(r, r ′ )χ(r ′ )E tot (r ′ )dr ′ (2.10)
where G(r, r ′ ) is Green's function, which represents the field at point r produced by a point source at r ′ .

Here, it is in 2-D case and given by

G(r, r ′ ) = -i 4 H (1) 0 (k B |r -r ′ |) = -i 4 (J 0 (k B |r -r ′ |) + iY 0 (k B |r -r ′ |)) (2.11)
where H

(1) 0 is the zero-order Hankel function of the first kind, J 0 and Y 0 are the zero-order Bessel and

Neumann functions, respectively.

In this work, the transverse magnetic (TM) polarization is considered. The object is immersed in a domain of interest (DoI) D contained in free space C. It is illuminated by TM waves generated by N s ideal line sources located outside the scatterer at positions r s . For each illumination, the scattered fields are collected by N r ideal receivers located at position r r along a circle of observation S, not intersecting the scatterer, and not necessarily jointing with the sources. By considering all those conditions, we now establish the forward model that involves two source-type integral equations.

The first one is called the state equation, which describes the interaction between the scatterer and the electromagnetic waves. The second one is the data equation, which describes the scattered fields generated by the scatterer. They have the following formulation

E tot (r, r s ) = E inc (r, r s ) + D G(r, r ′ )J(r ′ , r s )dr ′ , r ∈ D (2.12) E diff (r r , r s ) = D G(r r , r ′ )J(r ′ , r s )dr ′ , r ∈ D (2.13) with J(r, r s ) the contrast source J(r, r s ) = χ(r)E tot (r, r s ) (2.14)
The forward problem consists in calculating the scattered fields by solving the equations (2.12) and (2.13).

The forward model, the physical properties of the scatterer, and the domain knowledge are supposed to be well-known. Hence, it is a well-posed problem. The corresponding inverse problem is to retrieve the physical information of the scatterer, e.g., the contrast function χ(r), from the knowledge of the collected data E diff , and the known information E inc (r) and G(r, r ′ ).

Discretization of the forward model

To derive the discrete form of the equations, the Method of Moments (MoM) [START_REF] Davidson | Computational Electromagnetics for RF and Microwave Engineering[END_REF]) is used, which transforms continuous vector integral equations into discrete scalar algebraic equations through two steps: discretization and testing. For simplicity, a pulse-basis point matching is employed for the discretization, and the delta testing functions are adopted for the testing step.

Let L be a linear operator, f is a known function, g is an unknown function, we have

Lf = g (2.15)
One seeks to determine f knowing L and g, such that arg min

f |g -Lf | (2.16) Expand f with a series of basis functions {u n |n = 1, 2, . . . , N }, with α n decomposition coefficients f = ∞ n=1 α n u n ≈ N n=1 α n u n = f (2.17)
Combining (2.17) and (2.15), we have

g = L N n=1 α n u n = N n=1 α n Lu n (2.18)
In order to minimize the error between f and f , a set of testing functions should be introduced t 1 , t 2 , . . . , t K , then the vector inner product is performed

⟨t k , g⟩ G = N n=1 α n ⟨t k , Lu n ⟩ , ∀k = 1, . . . , K (2.19) where ⟨•⟩ G is the inner product on G such that ⟨t k , g⟩ G = G t k (x)g * (x)dx. The equation (2.19) can be rewritten as ω = Lα (2.20) with L =       ⟨t 1 , Lu 1 ⟩ ⟨t 1 , Lu 2 ⟩ . . . ⟨t 1 , Lu N ⟩ ⟨t 2 , Lu 1 ⟩ ⟨t 2 , Lu 2 ⟩ . . . ⟨t 2 , Lu N ⟩ . . . . . . . . . . . . ⟨t K , Lu 1 ⟩ ⟨t K , Lu 2 ⟩ . . . ⟨t K , Lu N ⟩       , α =     α 1 . . . α N     , ω =     ⟨t 1 , g⟩ . . . ⟨t K , g⟩     , (2.21)
Hence, we have α = L -1 ω, and

f = N n=1 (L -1 ω) n u n (2.22)
By applying the MoM with pulse basis and point matching method, a discretized version of previous integral equations can be obtained. The DoI D is discretized into N = N x ×N y small squares subunits, and the electric fields and permittivity are assumed to be constant within each unit. Every cell is approximated by a small disk with the same surface and equivalent radius R. The total field, the contrast function, and the contrast sources can be approximated as

E tot (r) = N i=1 E tot i C n (r) (2.23) χ(r) = N i=1 χ i C n (r) (2.24) J(r) = N i=1 J i C n (r) (2.25)
where C n (r) is the basis function

C n (r) =    0 r / ∈ D i 1 r ∈ D i (2.26)
where D i represents the i-th cell.

As previously mentioned, the test functions are a set of delta functions δ which are located at the center of the subcells. By applying the test functions on the incident fields

E inc (r) = N i=1 E inc i δ(r i -r), δ(r i -r) = 1 r = r i 0 r ̸ = r i (2.27)
where r 1 , . . . , r N are the centers of the sub cells. Then, decompose the state equation with the basis functions:

E tot (r, r s ) = E inc (r, r s ) + N i=1 D i G(r, r ′ )C n (r ′ )dr ′ J i , (2.28)
and one projects the test functions by Eq. ( 2.27)

E tot j = E inc j + N i=1 D i G(r, r ′ )C n (r ′ )dr ′ J i = E inc j + N i=1 D i G(r, r ′ )C n (r ′ )dr ′ χ i E tot i , j = 1, . . . N (2.29)
Now, we have a linearized system

E inc j = N i [δ ji -G ij χ i ]E tot i with G ij = D i G(r j , r ′ )dr ′ (2.30)
the above equation can be rewritten as

      E inc 1 E inc 2 . . . E inc N       =       1 -G 11 χ 1 G 12 χ 2 . . . G 1N χ N G 21 χ 1 1 -G 22 χ 2 . . . G 2N χ N . . . . . . . . . . . . G N 1 χ 1 1 -G N 2 χ 2 . . . 1 -G N N χ N       ×       E tot 1 E tot 2 . . . E tot N       , (2.31)
So, we have

E tot i = E inc i + G D J i , i = 1, . . . , N s (2.32)
with [START_REF] Richmond | Scattering by a dielectric cylinder of arbitrary cross section shape[END_REF])

G D (r, r ′ ) =        ik B πR 2 J 1 (k B R)H (1) 0 (k B |r -r ′ |) , r ̸ = r ′ ik B πR 2 H (1) 1 (k B R) -1 otherwise (2.33) where R = ∆ x ∆ y π , ∆
x and ∆ y are the length of small cell ∆.

Similarly, we have the observation equation for i r -th receiver

E diff ir (r) = N i=1 D i G(r, r ′ ) χ i E tot i dr ′ , (2.34) for N r receivers,       E diff 1 E diff 2 . . . E diff Nr       =       G 11 G 12 . . . G 1N G 21 G 22 . . . G 2N . . . . . . . . . . . . G Nr1 G Nr2 . . . G NrN       ×       χ 1 0 χ 2 0 . . . χ N       ×       E tot 1 E tot 2 . . . E tot N       (2.35)
Thus, we have (2.37) are the discretized version of the source-type integral equations, which will be used for all the approaches developed in this thesis.

E diff i = G S J i (2.36) J i = diag(χ)E tot i (2.37) where i = 1, . . . , N s , and G S (i r , r) = ik B πR 2 J 1 (k B R)H (1) 0 (k B |i r -r|). Equations (2.32), (2.36),

Contrast source inversion and wavelet technique

This chapter presents a multi-scale contrast source inversion method, based on the discrete wavelet transform (DWT). First, we present the details of Contrast Source inversion method (CSI). Then, the twodimensional DWT is introduced based on which a wavelet-based CSI method is developed. Finally, numerical results are displayed.

Contrast source inversion method

The CSI method is one of the most used methods to tackle the ISPs. It is based on the modified gradient method [START_REF] Kleinman | A Modified Gradient Method for Two-Dimensional Problems in Tomography[END_REF]. The cost function to minimize is built on source-type integral equations, which is the linear combination of a set of normalized mismatches. It is composed of two terms, one measures the error in the data equation, and the other measures the error in the object equation. The unknowns in the cost function are the contrast source J and the contrast χ. At each iteration, contrast source and contrast are updated alternatively, and the computation of the forward problem is avoided.

The cost functional stands as

F (J i , χ) = Ns i=1 ∥ξ i -G S J i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + Ns i=1 ∥diag(χ)E inc i -J i + diag(χ)G D J i ∥ 2 Ns i=1 ∥diag(χ)E inc i ∥ 2 (3.1)
ξ i being the measured scattered fields. In the following, diag(χ)E inc/tot i is written as χE inc/tot i for simplicity.

The contrast source J and the contrast χ are alternatively retrieved with a conjugate gradient iterative method. At each iteration, the contrast source J is first updated by minimizing the entire cost functional. Then, the contrast χ is obtained from the updated J , meanwhile minimizing the object error. The minimization at each iteration is linear. First, J and χ are initialized by back-propagation as

J bp i = ∥G * S ξ i ∥ 2 ∥G S G * S ξ i ∥ 2 G * S ξ i (3.2) χ bp = Ns i=1 J bp i E tot * i Ns i=1 ∥E tot i ∥ 2 . (3.3)
Here,

E tot i = E inc i + G D J bp i .
* is the conjugate transpose operator. At the k-th iteration, the data equation error ρ i,k and the state equation error r i,k are first determined

ρ i,k = ξ i -G S J i,k-1 (3.4) r i,k = χ k-1 E tot i,k-1 -J i,k-1 .
(3.5)

Then, the contrast source J i,k is updated by the conjugate-gradient method as

J i,k = J i,k-1 + α J k ν i,k (3.6) with ν i,k the search direction ν i,k = g J i,k + ℜ i g J i,k , g J i,k -g J i,k-1 i g J i,k-1 , g J i,k-1 ν i,k-1 (3.7)
where g J k is the Fréchet derivative evaluated at J k-1 and χ k-1 .

g J i,k = - G * S ρ i,k i ∥ξ i ∥ 2 - r i,k -G * D χ k-1 r i,k i ∥χ k-1 E inc i ∥ 2 (3.8)
and the step size α J k is defined as

α J k = -i g J i,k , v i,k i ∥G S ν i,k ∥ 2 i ∥ξ i ∥ 2 + i ∥ν i,k -χ k-1 G D ν i,k ∥ 2 i ∥χ k-1 E inc i ∥ 2 (3.9)
Then, the total field is updated by

E tot i,k = E inc i + G D J i,k (3.10)
The contrast χ is updated from J by an analytical solution

χ k = Ns i=1 J i,k E tot * i,k Ns i=1 ∥E tot i,k ∥ 2 (3.11)
Note that a positivity constraint is imposed on the contrast via the projection method (Van den [START_REF] Van Den Berg | Contrast Source Inversion Method: State of Art[END_REF]) by imposing ε r (r) ≥ 1 and σ(r) ≥ 0. The algorithm is sketched in Algorithm 1.

Algorithm 1 CSI

1: Initialize J and χ by back-propagation

Review of wavelet transform

The Fourier Transform (FT) [START_REF] Bracewell | The Fourier Transform and Its Applications[END_REF]) can transform a time-domain signal into the frequency domain. However, the frequency information can only be retrieved in a vague manner, i.e., the time when each frequency component appears is unknown, which can be explained by the Heisenberg uncertainty principle [START_REF] Busch | Heisenberg's Uncertainty Principle[END_REF]. The wavelet transform is one of the solutions to overcome this shortcoming of FT, as it can simultaneously yield time and frequency information with different resolutions.

To introduce the Wavelet transform (WT), we first define the decomposition of a given signal into the linear combination of a basis as

x(t) = k a k b k (t) (3.12)
with k an integer that could be finite or infinite, a k the decomposition coefficients, and b k (t) the basis. Consider bk (t) the dual function of b k (t), then we can find the decomposition coefficients by the orthonormal relation between the dual pairs

a k = x(t), bk (t) = x(t) b * k (t)dt (3.13)
where * is the complex conjugation.

Continuous wavelet transform

The WT was first introduced to provide a novel mathematical tool for seismic wave analysis. The Continuous wavelet transform (CWT) has the following formulation:

y(s, τ ) = x(t)Ψ * s,τ (t)dt (3.14)
where s is the continuous scaling parameter, and τ is the translation parameter. Ψ s,τ (t) represents the wavelets, which is a family of functions generated by applying translations and dilations on a "mother wavelet", defined as

Ψ s,τ (t) = 1 √ s Ψ t -τ s .
(3.15) Equation (3.14) signifies the decomposition of x(t) by the wavelets Ψ s,τ (t), and the inverse wavelet transform is given by

x(t) = y(s, τ )Ψ s,τ (

Discrete wavelet transform

In practice, the collected signals are discrete signals, which are generally non-continuous and contained within a limited interval. Discrete signals, e.g., digital images, should be tackled by the Discrete wavelet transform (DWT). To make the WT discrete, the wavelets should be discretized. They are not time-discrete, but discrete with respect to the translation and scale steps. Define the discretized wavelet function and scaling function as

Ψ j,k (r) = 2 j/2 Ψ(2 j r -k) (3.17) Φ j,k (r) = 2 j/2 Φ(2 j r -k) (3.18)
where j is the dilation parameter, and k is the position parameter. The scaling function Φ(r) is used for approximating the signals and the wavelet function Ψ(r) is used to retrieve the details. Consider a discrete signal x(r), it can be decomposed by

x(r) = k W Φ (j 0 , k)Φ j 0 ,k (r) + ∞ j=j 0 k W Ψ (j, k)Ψ j,k (r) (3.19)
and the wavelet coefficients can be obtained by

W Φ (j 0 , k) = ⟨x(r), Φ j 0 ,k (r)⟩ (3.20) W Ψ (j, k) = ⟨x(r), Ψ j,k (r)⟩ , j ≥ j 0 (3.21)
W Φ (j 0 , k) are approximation coefficients and W Ψ (j, k) are detail coefficients. In practice, the scaling function and wavelet function are known, e.g., Haar, Daubechies, .etc. (Chui 1992).

The Haar wavelet is a very simple wavelet, and it is widely adopted to analyze the discrete signal, its mother wavelet function is

Ψ(x) =          1 0 ≤ x < 1 2 -1 1 2 ≤ x < 1 0 otherwise. (3.22)
and its scaling function can be described as

Φ(x) =    1 0 ≤ x < 1 0 otherwise. (3.23)
The Daubechies wavelets are a family of orthogonal wavelets dedicated to the DWT, which are characterized by a higher number of vanishing moments A (which is also called tap). The length of the Daubechies wavelets is related to their vanishing moments. For a Daubechies wavelet with A vanishing moments, its wavelet filter and scaling filter are both of length 2A. When A = 1, the corresponding Daubechies wavelet is denoted as D2 or dB1, it is the Haar wavelet. As shown in Fig. 3.1, an image X to be decomposed is firstly passed through two filters: low-pass filter and high-pass filter. The coefficients of the low-pass filters for D2 -20 are detailed in Appendix A.

The corresponding wavelet coefficients can be derived from the scaling coefficients. The low-pass filter is used to approximate the image, and the other can extract information on edges. Then, the image can be decomposed into two parts: low-frequency components and high-frequency components. Next, the two sub-images repeat the previous operation, respectively. Finally, we can get 4 sub-images of level 1, corresponding to 4 spatially oriented frequency channels:

• X LL : A re-scaled approximation of the original image • X HL : Horizontal details information • X LH : Vertical details information • X HH : Diagonal details information By repeating the 1-level decomposition on the LL components, we can obtain wavelet coefficients of higher levels, as shown in Fig. 3.2. The frequency resolution is increased with the decomposition level increasing.

Figure 3.3 shows the location of wavelet coefficients in the transformed matrix. For the coefficients of level

Multi-scale wavelet based CSI

As mentioned above, applying the 2-D DWT on a given image produces wavelet coefficients of different levels corresponding to different frequency resolutions. Incorporating the DWT into classical inversion methods is attracting much attention due to its multiresolution representation. Our work aims to build up a CSI method in the wavelet domain, combined with the multiresolution property, to see which configuration might outperform the classical CSI approach in the spatial domain.

In [START_REF] Li | A Contrast Source Inversion Method in the Wavelet Domain[END_REF], the authors have transformed the cost function of CSI into the wavelet domain.

ρ W i,k = ξ i -G S W * {γ i,k-1 } (3.26) r W i,k = W{W * {β k }E inc i } + W{W * {β k }G D W * {γ i }} -γ i (3.27)
Thus, with all elements defined above, the cost function of CSI can also be transformed into the wavelet domain

F W (γ i , β) = Ns i=1 ∥ρ W i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + ∥r W i ∥ 2 Ns i=1 ∥W{W * {β}E inc i }∥ 2 (3.28) Compared to the original cost function F (J i , β), F W (γ i , β
) is a function of γ i and β instead of J i and χ.

At each iteration, the contrast source and the contrast are first transformed into the wavelet domain, so we can operate on their wavelet coefficients. Then, the spatial domain coefficients are reconstructed from the modified wavelet coefficients.

As introduced above, on the one hand, the wavelet coefficients include two main parts: approximation coefficients and detail coefficients. Approximation coefficients contain low-frequency components, while the detail coefficients contain high-frequency components which are expected to be more contaminated with noise. Therefore, the approximation coefficients are more robust to noise and they concentrate most of the energies. On the other hand, different decomposition levels L yield decomposition coefficients at different scales.

Based on these two properties, a new strategy can be established. Here, only approximation coefficients are concerned. We start from the reconstruction of approximation coefficients at the finest scale (the maximum level of the decomposition L max ). Once the reconstruction at the current level L is finished, the result will be used as the initialization for the next step, where the approximation coefficients at scale L = L -1 are updated. This step is repeated until L reaches 1. Here, only the wavelet coefficients of the contrast χ are concerned, so we can simplify the equation (3.28) as

F W (J i , β) = Ns i=1 ∥ρ i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + Ns i=1 W * {β}E inc i + W * {β}G D J i -J i 2 Ns i=1 W * {β}E inc i 2 (3.29)
Compared to equation (3.28), the above equation is a function of J i and β. At each iteration, only the contrast has to be decomposed, which reduces the computational burden. The algorithm of the WCSI is summarized in Algorithm 2.

Numerical results

Test with "Austria" profile

In the numerical simulation, the so-called "Austria" profile is used, which contains two disks and one ring, as shown in Fig. 3.5(a). The scatterer is lossless, the relative permittivity ε(r) is 2 for the object, and for the embedding medium, ε b = 1. The region of interest D is l = 2 m sided. For the direct problem, D is discretized into 50 × 50 cells, and for the inverse one, it is discretized into 64 × 64 pixels. The scattering system is shown in

Update the contrast β k by β k = W      Ns i=1 J i,k E tot * i,k Ns i=1 E tot i,k 2      12:
Keep only the approximation coefficients As the region of interest is discretized into 64 × 64 pixels, the maximum level of decomposition that can be reached is L max = 6, which means that a single WCSI result is obtained after 6 successive inversions. For CSI, the algorithm stops when the maximum number of iterations is reached. Here, the maximum number of iterations is set as 500. For WCSI, 200 iterations are conducted for each inversion.

Permittivity

-1 0 1 x[m] -1 -0.5 0 0.5 1 y[m]
To compare the reconstruction performance of CSI and WCSI, maps of the normalized contrast error as a function of the number of sources and receivers are given in Fig. 3.6 -3.7. The number of sources N s and number of receivers N r vary from 6 to 48.

Moreover, we compare the reconstruction performance of wavelets with various vanishing moments within the same Daubechies family. Three wavelets are chosen: Haar, db4, and db14. Fig. 3.6 shows the relative contrast error obtained by CSI and WCSI with three wavelet bases. Here, the Gaussian noise of 20 dB is added to the scattered fields. From Fig. 3.6, we can observe that both CSI and WCSI with Haar wavelets work well when there is a sufficient number of sensors. However, WCSI achieves better performance as N s and N r decrease. The results obtained by WCSI with db4 and db14 wavelets are displayed in Fig. 3.6(c) and Fig. 3.6(d). We can observe that using db4 and db14 wavelets can further improve the reconstruction quality of WCSI. The more the vanishing moments of the wavelet, the higher the order of its filters. With such filters, the energy of the decomposition coefficients becomes more compact and concentrated. Also, the scaling function and wavelet function become smoother, which improves the ability of the wavelet to approximate a signal.

To investigate the robustness of the proposed approach, we further test CSI and WCSI for cases where 10 dB noise affects the collected data; the results are shown in Fig. 3.7. We can observe from the results that the reconstruction performance of CSI strongly degrades with stronger noises when the number of sensors is reduced. On the contrary, WCSI is more robust to noise, especially with the db4 and db14 wavelets (Fig. 3.7(c) and Fig. 3.7(d) ). Illustrative results obtained by CSI and WCI applied to cases with different N s , and N r are shown in Fig. 3.9 and Fig. 3.10. The relative reconstruction errors are given in Table 3.1. Four combinations of N s and N r are considered: N s × N r = 48 × 48, 18 × 12, 12 × 9, and 6 × 18.

In Fig. 3.9, the results are obtained from measurement data affected by a Gaussian noise of 20 dB. The obstacle can be well reconstructed by all methods for the first configuration with enough sensors. As shown in Table 3.1, the relative errors obtained by CSI, WCSI-db4, and WCSI-db14 are less than 10%. Note that the WCSI with Haar wavelet basis has the highest reconstruction error because of the loss of the detail information.

When N s and N r decrease, the reconstruction quality of the results obtained by CSI degrades significantly. For the configuration 18×12, CSI cannot obtain meaningful results. For the configuration N s ×N r = 12 × 9, the obstacle cannot be well retrieved by CSI; one of the disks and a part of the ring are missing. In comparison, the results obtained by WCI have a much better reconstruction quality. When N s ×N r = 6×18, the obstacle can be roughly reconstructed by CSI. Still, WCSI performs better than CSI for all three wavelet basis. The db14 basis outperforms the other two bases, as illustrated in Table 3

Test with circular dataset

The proposed method is further tested on a synthetically generated circular dataset (Wei et al. 2019a).

Each sample contains disks of random quantities: number (between 1 and 3), relative permittivity (between 1.5 and 2.5), radius (between 0.15 m and 0.4 m), and location. The scatterers may overlap and are lossless. The dataset contains 250 samples. The configuration of the scattering system is the same as in the previous section. The number of sensors is of N s × N r = 6 × 18.

N s × N r = 48 × 48 N s × N r = 18 × 12 N s × N r = 12 × 9 N s × N r = 6 × 18
N s × N r = 48 × 48 N s × N r = 18 × 12 N s × N r = 12 × 9 N s × N r = 6 × 18

Conclusion

From all the above results, WCSI can improve the quality of reconstruction compared to CSI, especially when the information is reduced. Also, WCSI is more robust than CSI. However, as WCSI has multiple optimization stages corresponding to different decomposition levels, it is more time-consuming than CSI.

Introduction

Due to the strong ill-posedness of ISPs, additional penalty terms are usually needed to regularize the optimization faced. To that intent, various regularization techniques can be applied, such as positivity, total variation (TV), and sparsity, as previously introduced.

Traditional computational methods to tackle ISPs, e.g., CSI and DBIM, have already been combined with various regularization techniques. In [START_REF] Berg | Extended Contrast Source Inversion[END_REF], the TV regularization is implemented within the CSI framework as a multiplicative regularization term that leads to a better reconstruction quality and improves robustness to noise. In the DBIM algorithm, Tikhonov regularization [START_REF] Tikhonov | Solution of Incorrectly Formulated Problems and the Regularization Method[END_REF]) is employed by default in order to circumvent the instability of the problem. Edge-preserving regularization has also been introduced into ISPs [START_REF] Lobel | A New Regularization Scheme for Inverse Scattering[END_REF][START_REF] Qin | Joint Inversion of Electromagnetic and Acoustic Data with Edge-preserving Regularization for Breast Imaging[END_REF] to improve the reconstruction quality with respect to noisy corrupted data while preserving the edges of the images. In [START_REF] Qin | Joint Inversion of Electromagnetic and Acoustic Data with Edge-preserving Regularization for Breast Imaging[END_REF], Huber regularization is incorporated with the CSI method as a smoothness constraint.

Among different regularization techniques, sparsity is a widely investigated choice.

In [START_REF] Taşkin | Sparsity Regularized Nonlinear Inversion for Microwave Imaging[END_REF], the fast iterative shrinkage-thresholding algorithm (FISTA) [START_REF] Beck | A Fast Iterative Shrinkage-thresholding Algorithm for Linear Inverse Problems[END_REF]) is combined with a modified gradient method (MGM) to enforce the sparsity in both non-sparse and sparse domains, where DWT is employed to reach a sparse representation of the unknown.

In [START_REF] Zaimaga | Sparsity-enforced Microwave Inverse Scattering Using Soft Shrinkage Thresholding[END_REF], an iterative algorithm incorporating a Tikhonov functional with a sparsity promoting ℓ 1 term is proposed, which enables a sharp reconstruction in the sparse domain.

Iterative shrinkage thresholding algorithm (ISTA) [START_REF] Beck | A Fast Iterative Shrinkage-thresholding Algorithm for Linear Inverse Problems[END_REF]) is implemented within the framework of a boundary integral method in [START_REF] Desmal | Shrinkage-thresholding Enhanced Born Iterative Method for Solving 2D Inverse Electromagnetic Scattering Problem[END_REF] to solve the optimization problem with zeroth/firstnorm penalty term. The proposed method is proved to be able to yield a sharper and more accurate reconstruction while benefiting from a faster convergence speed.

In [START_REF] Sanghvi | Iteratively Reweighted ℓ 1 -ℓ 2 Norm Minimization Using Wavelets in Inverse Scattering[END_REF], a joint ℓ 1 -ℓ 2 norm regularized BIM approach in the wavelet domain is proposed, where appropriate weights adaptively updated at each iteration are allocated to wavelet coefficients at different levels. This scheme enables the performance to be independent of the initial weights. In [START_REF] Winters | A Sparsity Regularization Approach to the Electromagnetic Inverse Scattering Problem[END_REF]) and [START_REF] Anselmi | Microwave Medical Imaging Based on Sparsity and an Iterative Method with Adaptive Thresholding[END_REF], DBIM is incorporated with ℓ 2 -ℓ 1 and ℓ 1 -ℓ 0 joint norm, respectively, solved by various optimization techniques [START_REF] Zou | Regularization and Variable Selection Via the Elastic Net[END_REF][START_REF] Azghani | Iterative Algorithms for Random Sampling and Compressed Sensing Recovery[END_REF]. In [START_REF] Bevacqua | Non-linear Inverse Scattering Via Sparsity Regularized Contrast Source Inversion[END_REF], two sparsity-enhanced approaches within the CSI framework are proposed, the ℓ 1 norm being used as weighted penalty term and constraint. Therein, the stationary wavelet transform (SWT) [START_REF] Nason | The Stationary Wavelet Transform and Some Statistical Applications[END_REF]) is adopted to project the unknown onto a redundant dictionary.

Compressive sensing (CS) based reconstruction algorithms have also been exploited [START_REF] Oliveri | Compressive Sensing As Applied to Inverse Problems for Imaging: Theory, Applications, Current Trends, and Open Challenges[END_REF][START_REF] Massa | Compressive Sensing in Electromagnetics -a Review[END_REF]. In [START_REF] Anselmi | Microwave Medical Imaging Based on Sparsity and an Iterative Method with Adaptive Thresholding[END_REF], the commonly used DWT offers a sparse representation of the profile to be reconstructed, then the linearized inverse scattering problem within the first-order Born approximation framework is solved with the Bayesian version CS-based procedure.

The aforementioned state-of-the-art generally uses ℓ 0 , ℓ 1 , ℓ 2 norms or their combinations to add sparsity constraint to the optimization problem faced with. These vector norms can regularize each element individually regardless of its spatial position and the potential connections between different elements. When the elements are not only sparse but also have potential structure features, the structure information cannot be extracted using these single norms.

Mixed norms [START_REF] Kowalski | Sparse Regression Using Mixed Norms[END_REF]) such as ℓ 2,1 or ℓ ∞,1 are usually employed to impose both sparsity and structural prior information in the optimization problem [START_REF] Tomassi | Wavelet Shrinkage Using Adaptive Structured Sparsity Constraints[END_REF]. In [START_REF] Zaimaga | Sparse Reconstruction Algorithms for Nonlinear Microwave Imaging[END_REF], as the row vectors of the contrast source share a common non-zero support, a two-step sparse recovery process is investigated where the ℓ 2,1 norm is to enforce joint sparsity on the contrast source.

In this chapter, a group sparsity penalized CSI method is proposed, which employs the mixed ℓ 2,1 norm to impose the group sparsity on the contrast in the wavelet domain. Herein, we use the dual-tree complex wavelet transform (DT-CWT) [START_REF] Selesnick | The Dual-tree Complex Wavelet Transform[END_REF]) instead of DWT because of a better wavelet tree structure provided by CWT. In the following, we first detail the DT-CWT. Then, the wavelet tree structure [START_REF] Baraniuk | Optimal tree approximation with wavelets[END_REF]) is introduced, which exists in the wavelet domain, and led by the parent-child relationship [START_REF] Liu | Information-theoretic Analysis of Interscale and Intrascale Dependencies between Image Wavelet Coefficients[END_REF] between wavelet coefficients across different scales. Next, the proposed approach is introduced and numerical illustrations are presented.

Dual-tree Complex Wavelet Transform

The DT-CWT is a relatively recent enhancement of DWT, with some important additional properties. As previously introduced, decomposing a signal X by DWT is to recursively apply the discrete-time low pass and high pass filters, as well as up-sampling and down-sampling operations. DWT has been widely used for signal analysis and reconstruction, which provides a sparse representation of a given signal, and enables multiscale analysis. However, there are four main shortcomings as follows:

1. Oscillations: The wavelet transform can provide a sparse representation for a given signal because the wavelets oscillate locally, and large wavelet coefficients are produced when the wavelets overlap a singularity. Otherwise, small wavelet coefficients are produced. However, the wavelet coefficients tend to oscillate in a positive and negative manner around the singularities, which makes the waveletbased processing more complicated.

Shift variance:

The wavelet coefficients of a given signal are sensitive to translations of the input signal. A small shift of the input signal results in an obvious change in the wavelet domain, i.e., the distribution of energy between DWT coefficients at different scales. This property also complicates the wavelet-domain signal processing.

3. Aliasing: Substantial aliasing appears during the decomposition process, and the inverse DWT cancels this aliasing. The time-domain signal can be perfectly reconstructed from its wavelet coeffi-cients thanks to the specially designed decomposition and reconstruction filter banks. However, any wavelet coefficients processing" e.g. thresholding, filtering, and quantization, upset the balance between the direct and inverse transforms, resulting in artifacts in the reconstruction.

4. Lack of directionality: When employing DWT to decompose a 2-D image, the filters and down-sampling operations are applied separately along row and column directions and four subimages can be obtained: LL, LH, HL and HH. The subbands LH, HL and HH contain detail information in three directions.

However, this weakens the detail information in other directions, making it difficult to represent nonhorizontal or non-vertical singular features. This drawback greatly complicates the modeling and processing of geometric image features such as ridges and edges.

Inspired by the FT based on complex-valued oscillating sinusoids, the DT-CWT was introduced by Kingsbury in 1998 [START_REF] Selesnick | The Dual-tree Complex Wavelet Transform[END_REF]. Unlike the conventional DWT which uses real-valued oscillating wavelets, the DT-CWT uses complex-valued scaling function and wavelet

Φ c (r) = Φ r (r) + jΦ i (r) (4.1) Ψ c (r) = Ψ r (r) + jΨ i (r) (4.2)
Projecting the signal onto 2 j/2 Ψ(2 j r -k), we can obtain the complex wavelet coefficients

W c (j, k) = W r (j, k) + jW i (j, k) (4.3)
If Ψ r (r) and Ψ i (r) form a Hilbert transform pair, then Ψ c (r) is an analytic signal. The complex scaling function is defined similarly. The complex wavelet and scaling functions are expected to be as close as possible to analytic in order to perfectly overcome the four DWT shortcomings. The DT-CWT uses a dualtree structure to obtain the real and imaginary parts of the complex wavelet coefficients. Two real DWTs are employed and each tree corresponds to a DWT. The first DWT gives the real part of the complex coefficients and the second one gives the imaginary part. show the analysis and synthesis filter banks of the DT-CWT. h 0 (n), h 1 (n) denote the low pass/high pass filter pair for the upper filter bank (real branch). g 0 (n), g 1 (n) denote the low pass/high pass filter pair for the lower filter bank (imaginary branch). The wavelet functions Ψ h (r) and Ψ g (r) associated with the two DWTs satisfy Equation (4.2), and Ψ g (r) is designed to be approximately the Hilbert transform of Ψ h (r) so that the corresponding complex wavelet Ψ c (r) is approximately analytic. The inverse of DT-CWT is also the synthesis of two real DWTs, which correspond to the inverse of the two real DWTs used in the forward transform. g(n) is the inverse of g(n) and h(n) is the inverse of h(n). It should be noted that the implementation of the DT-CWT requires that the filters used in the first stage should be different from the succeeding stages, in order to make the filter banks approximately analytic.

Comparison of 2-D DWT and 2-D DT-CWT

As previously introduced, 2-D DWT can produce three sub-bands which contain detail information: LH, HL and HH, obtained by three wavelets characterized by

Ψ LH (x, y) = Φ(x)Ψ(y) (4.4) ℎ 0 1 (𝑛) ℎ 1 1 (𝑛) 2 2 ℎ 0 (𝑛) ℎ 1 (𝑛) 2 2 ℎ 0 (𝑛) ℎ 1 (𝑛) 2 2 𝑔 0 1 (𝑛) 𝑔 1 1 (𝑛) 2 2 𝑔 0 (𝑛) 𝑔 1 (𝑛) 2 2 𝑔 0 (𝑛) 𝑔 1 (𝑛) 2 2
(a) Analysis filter bank Ψ HL (x, y) = Ψ(x)Φ(y)

෨ ℎ 0 1 (𝑛) ෨ ℎ 1 1 (𝑛) ෨ ℎ 0 (𝑛) ෨ ℎ 1 (𝑛) ෨ ℎ 0 (𝑛) ෨ ℎ 1 (𝑛) 2 𝑔 0 1 (𝑛) 𝑔 1 1 (𝑛) 𝑔 0 (𝑛) 𝑔 1 (𝑛) 𝑔 0 (𝑛) 𝑔 1 (𝑛)
(4.5) Ψ HH (x, y) = Ψ(x)Ψ(y) (4.6)
The LH wavelet is the product of the scaling function Φ along the x-direction and the wavelet function Ψ along the y-direction. The HL and HH wavelets are similarly defined. Fig. 4.2 shows the three wavelets. We can observe that the LH and HL wavelets are oriented vertically and horizontally. However, a checkerboard artifact appears in the HH wavelet, which is due to the use of real-valued separate filters. This makes the wavelet unable to distinguish +45 • and -45 • orientations.

Differently from the 2-D DWT, the 2-D DT-CWT produces six sub-bands for both real and imaginary parts.

The transformation result is four times more expansive than that of DWT. The six wavelets for the real part are defined in the following way

Ψ i (x, y) = 1 √ 2 (Ψ 1,i (x, y) -Ψ 2,i (x, y)) (4.7) Ψ i+3 (x, y) = 1 √ 2 (Ψ 1,i (x, y) -Ψ 2,i (x, y)) (4.8) where i = 1, 2, 3, 1 √ 2
is the normalization factor to make the sum/difference operation constitutes an orthonormal operation. The wavelet bases are defined in the same manner as equations (4.4)-(4.6)

Ψ 1,1 (x, y) = Φ h (x)Ψ h (y), Ψ 2,1 (x, y) = Φ g (x)Ψ g (y) (4.9) Ψ 1,2 (x, y) = Ψ h (x)Φ h (y), Ψ 2,2 (x, y) = Ψ g (x)Φ g (y) (4.10) Ψ 1,3 (x, y) = Ψ h (x)Ψ h (y), Ψ 2,3 (x, y) = Ψ g (x)Ψ g (y) (4.11)
Likewise, the six wavelets for the imaginary parts are defined as

Ψ i (x, y) = 1 √ 2 (Ψ 3,i (x, y) -Ψ 4,i (x, y)) (4.12) Ψ i+3 (x, y) = 1 √ 2 (Ψ 3,i (x, y) -Ψ 4,i (x, y)) (4.13)
where i = 1, 2, 3 and the wavelet bases are defined as

Ψ 3,1 (x, y) = Φ g (x)Ψ h (y), Ψ 4,1 (x, y) = Φ h (x)Ψ g (y) (4.14) Ψ 3,2 (x, y) = Ψ g (x)Φ h (y), Ψ 4,2 (x, y) = Ψ h (x)Φ g (y) (4.15) Ψ 3,3 (x, y) = Ψ g (x)Ψ h (y), Ψ 4,3 (x, y) = Ψ h (x)Ψ g (y) (4.16)
With the twelve oriented wavelets defined above for the real part and imaginary part, a set of complex wavelets can be constructed.

2-D separable wavelets (db4) --LH, HL, HH

Wavelet quadtree structure

The wavelet transform has several attractive characteristics [START_REF] Mallat | A Wavelet Tour of Signal Processing, Third Edition: The Sparse Way[END_REF]):

• Edge detection and sparsity: the wavelet transform is an efficient tool to capture singularities of signals, in 2-D cases, the extraction of image features like ridges and edges. As wavelets oscillate locally, only wavelets overlapping a singularity have large wavelet coefficients and all other coefficients are small. This is also the key to the sparsity of wavelet coefficients because the edge is only a small part of the image content.

• Multiscale analysis: the wavelet transform allows an image to be analyzed at different resolution scales since the transformation domain has a nested set of scales.

These attractive characteristics can be regarded as the primary properties, from which a secondary property can be derived:

• Persistence: the wavelet transform naturally yields a dependent quadtree structure. This dependency tends to propagate across different scales. Consequently, wavelet coefficients across different scales can be regarded as statistically dependent. If a wavelet coefficient is large/small, the coefficients at the same orientation in the next adjacent scale are likely to be large/small. This is the so-called parent-child relationship. [START_REF] Crouse | Wavelet-based Statistical Signal Processing Using Hidden Markov Models[END_REF]).

This secondary property has been widely exploited in the statistical signal processing community [START_REF] Crouse | Wavelet-based Statistical Signal Processing Using Hidden Markov Models[END_REF][START_REF] Romberg | Bayesian Tree-structured Image Modeling Using Waveletdomain Hidden Markov Models[END_REF], generally modeled as a hidden Markov tree (HMT). be regarded as the parent of its children coefficients in the next scale. This leads to a hierarchical quadtree structure. More precisely, for a detail coefficient located at (x, y) at the decomposition level L, its four children will be in the next decomposition level L -1, located respectively at (2x -1, 2y -1), (2x -1, 2y), (2x, 2y -1) and (2x, 2y), occupying a 2 × 2 region. The arrows in Fig. 4.6(a) identify the parent-children

dependencies.

An example of a wavelet quadtree is given in Fig. 4.6(b). Consider that the node 1 is a coefficient in {HL, LH, HH} 3 and make it a root node. A quadtree can be constructed based on this root node, where its four children are placed in the second layers. The four children coefficients lead to sixteen "grandchildren" of the root node in the third layer. Emphasis should be put on the parent-child relationship led by the wavelet quadtree structure. Since dependency relations exist between wavelet coefficients, i.e., if a parent coefficient is negligible, then its children are also generally negligible. The wavelet coefficients can therefore be divided into different groups, based on which the group sparsity technique can be applied.

Nevertheless, the 2-D DWT suffers from a strong shift-variance, which greatly affects the persistence property of wavelet coefficients. It has been proved in [START_REF] Choi | Hidden Markov Tree Modeling of Complex Wavelet Transforms[END_REF] that the near shift invariance of the oriented DT-CWT leads to better persistence of the magnitudes around the edgy regions, which ensures stronger dependence in inter-scale adjacent coefficients. In other words, the oriented DT-CWT strengthens the parent-child relationship, which is the key to implementing group-sparsity-based technique.

Sparsity inducing norms

In the previous section, the DT-CWT and the wavelet quadtree structure are introduced, which are two key points to achieve the group sparsity penalized CSI method. In this section, the sparsity regularization terms that enable the incorporation of prior information into optimization problems are sketched.

Sparsity through single norm

A vector x of dimension M can be considered as sparse if most of its elements are zero, or almost zero.

For an optimization problem referring to x = arg min

x F (4.17) the sparsity can be directly achieved by adding a regularization term that penalizes the number of non-zero elements of the vector, in the form of the vector norm ∥x∥ q . Thus the combinatorial optimization arises

x = arg min x F + ∥x∥ q (4.18)

To impose sparsity as the prior information into the optimization problem, there are several choices of q:

• When q = 0, the regularization term is the ℓ 0 norm, defined as ∥x∥ 0 = #{i = 1, . . . , M |x i ̸ = 0}, which counts the non-zero elements of the vector. However, solving optimization problems with the non-convex ℓ 0 norm is as difficult as trying out all possible non-zeros supports of coefficients. To avoid this NP-hard problem, a tractable approach to solving the ℓ 0 minimization problem is to replace it with a reasonable and convex approximation.

• When q = 1, the regularization term is the ℓ 1 norm, which is a natural approximation of ℓ 0 norm, defined as ∥x∥ 1 = M i=1 |x i |. Unlike ℓ 0 norm, the optimization of ℓ 1 norm is a convex problem and can be solved efficiently [START_REF] Boyd | Convex Optimization[END_REF]).

Sparsity through mixed norm

In the previous formulation with ℓ 1 norm, every element of the vector is regarded as independent from the others. In some situations, however, the elements of a vector share the same support, or dependency relationships exist between different elements, which allows them to be divided into subgroups. It is then easy to see that the elements within a given group can be selected or deleted simultaneously, benefiting from this prior structure information. A regularization norm exploiting such structural prior information proposed in [START_REF] Yuan | Model Selection and Estimation in Regression with Grouped Variables[END_REF]) has the following form

∥x g ∥ 2,1 = s i=1 ∥x g i ∥ 2 (4.19)
g represents the subgroups, x g i are the sub-vectors of x indexed by g i , i = 1, . . . , s where s is the number of subgroups. This is the so-called mixed ℓ 2,1 norm, acting like an ℓ 1 norm on the vector ∥x g i ∥ 2 . It should be noted that the mixed ℓ 2,1 norm does not yield sparsity within a group. That is, if most of the elements in a group are non-zero, then they will all be imposed to be non-zero. ℓ ∞,1 norm is also a popular one to induce group sparsity. In this work, the ℓ 2,1 is used.

Solution to optimization problems with regularization terms

To solve a non-constrained optimization problem where the cost function to minimize is differentiable, a common way is to use the iterative gradient descent method. At each iteration, x is updated according to

x k = x k-1 -t k ∇F(x k-1 ) (4.20)
k represents the k-th iteration, t k is the step size and ∇F(x k-1 ) is the gradient of F at the point k -1. However, due to the presence of the regularization norm, the whole cost function becomes non-differentiable, so the existence of the gradient is not always guaranteed.

Subgradient method One possible manner to deal with this kind of optimization problem is to use the subgradient technique [START_REF] Shor | Minimization Methods for Non-Differentiable Functions[END_REF], which allows approximating the gradient of the non-differentiable function at non-derivable points. For a given function F, ∀y, g is the subgradient of F at point x if it satisfies

F(y) ≥ F(x) + g T (y -x) (4.21)
Convex functions always have sub-gradients. Provided that the left-hand derivative of a convex function at

point x 0 is a = lim x→x - 0 f (x) -f (x 0 ) x -x 0 (4.22)
and the right-hand derivative is b = lim

x→x + 0 f (x) -f (x 0 ) x -x 0 (4.23)
the subgradient of the convex function can be any value in [a, b]. Non-convex functions do not necessarily have sub-gradients even if they are differentiable. However, the convergence rate of the subgradient method is very slow.

Proximal gradient descent method What one adopts in the present work is another well-known approach called the proximal gradient descent method [START_REF] Rockafellar | Convex Analysis[END_REF]. Now consider the function can be decomposed into two functions as follows

F(x) = G(x) + H(x) (4.24)
where G is convex, differentiable, H is convex and not necessarily differentiable. For the optimization problem

x = arg min x {G(x) + H(x)} (4.25)
x can be updated at iteration k according to

x k = prox tH( ) x k-1 -t k ∇G x k-1 (4.26)
wherein t k is the descent step size at k-th iteration. The proximal operator of H( ) at point x has the following expression and it is a function of H and t

prox tH( ) (t) = arg min x 1 2t ∥x -t∥ 2 2 + H(x) (4.27)
The key advantages of the proximal gradient descent method are as follows:

• Proximal operator can be computed analytically for plenty of widely used H functions. For example, when H(x) = ∥x∥ 1 , the corresponding prox tH( ) ( ) is the soft-thresholding operator. • Proximal operator only depends on H, this means that G can be a complicated function.

Solution to ℓ 2,1 norm penalized optimization problems

When H(x) = ∥x∥ 2,1 , this is the so-called sparse group Lasso [START_REF] Simon | A sparse-group lasso[END_REF]) -least absolute shrinkage and selection operator-problem, and the solution depends on the grouping methods.

For non-overlapping groups Non-overlapping means that each element belongs to only one group. To tackle the non-overlapping group Lasso problem where the ℓ 2,1 norm is applied, the proximal gradient descent method is usually employed because the proximal operator of ℓ 2,1 has a closed-form solution [START_REF] Kowalski | Sparse Regression Using Mixed Norms[END_REF] with λ the hyperparameter linked to the regularization term

x ′ g i = x g i × 1 - λ ∥x g i ∥ 2 +
, ∀i For overlapping groups In overlapping group cases, an element can appear in different subgroups. For instance, the parent-child relationship leads to a hierarchical tree structure, which is a typical overlapping grouping arrangement. Take the example of the wavelet quadtree structure displayed in Fig. 4.7(b), the nodes {1, 2, 3, 4, 5} constitute a subgroup, and the nodes {5, 18, 19, 20, 21} form another subgroup. Here, the node 5 appears in two subgroups.

For the overlapping group Lasso problem, where an element can appear in different subgroups, the problem becomes more complicated since the standard group Lasso optimization strategies cannot be used directly.

Much attention has been paid to exploiting overlapping situations, especially those with a hierarchical tree structure. In [START_REF] Liu | Moreau-Yosida Regularization for Grouped Tree Structure Learning[END_REF], an approach to deal with tree-structured group Lasso problems is offered, in which the Moreau-Yosida regularization is employed, and its analytical solution is given. In (Jenatton et al. 2011), the dual approach is used, and the proximal operator of the norm associated with overlapping hierarchical groups proved to be exactly computable, which can be regarded as the composition of elementary proximal operators.

The strategies for overlapping group optimization problems mentioned above are complicated to implement. A smart way to transform the overlapping group Lasso problems into non-overlapping ones is proposed in [START_REF] Jacob | Group Lasso with Overlap and Graph Lasso[END_REF][START_REF] Rao | Convex Approaches to Model Wavelet Sparsity Patterns[END_REF][START_REF] Chen | The Benefit of Tree Sparsity in Accelerated MRI[END_REF]. The authors propose the so-called replication strategies, which consist in decoupling the overlapping groups from each other by replicating the elements appearing repeatedly. For example, by applying the replication strategy, the two aforementioned subgroups {1, 2, 3, 4, 5} and {5, 18, 19, 20, 21} are decoupled, and the node 5 is replicated. Then, the groups become {1, 2, 3, 4, 5 a } and {5 b , 18, 19, 20, 21}. The decoupled groups can tackle this problem using the same proximal operator as in the non-overlapping case.

Group sparsity penalized CSI

In this section, the proposed group sparsity penalized CSI approach is sketched.

ℓ 2,1 norm regularized cost function

Here, the contrast function χ is decomposed by the 2-D oriented DT-CWT. Owing to the two properties of wavelet coefficients: sparsity and persistence, the group sparsity regularization is imposed in the cost function of the CSI method in the wavelet domain.

Take the cost function of CSI in the wavelet domain

F W (J i , β) = Ns i=1 ∥ρ i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + Ns i=1 W * {β}(E inc i + G D J i ) -J i 2 Ns i=1 ∥W{W * {β}E inc i }∥ 2 (4.29)
Here, W{ } represents the 2-D oriented DT-CWT instead of conventional 2-D DWT, and W * { } is the inverse transform. Employing the DT-CWT, the contrast χ is decomposed in the complex wavelet domain at each iteration. β = W{χ} are the complex wavelet coefficients of the contrast function. Now, one adds the group sparsity regularization term to the cost function (4.29), i.e., the ℓ 2,1 norm, and the cost function becomes

F W (J i , β) = Ns i=1 ∥ρ i ∥ 2 Ns i=1 ∥ξ i ∥ 2 + Ns i=1 ∥W * {β}(E inc i + G D J i ) -J i ∥ 2 Ns i=1 ∥W{W * {β}E inc i }∥ 2 + Λ∥β g ∥ 2,1 (4.30)
Subgroups can be constituted based on the dependency between parent and children coefficients. Thus, the group sparsity regularization can be applied to the contrast in the wavelet domain, which promotes the persistence of magnitudes of wavelet coefficients across scales. If most of the coefficients in a subgroup are almost zero or nearly zero, then the whole group will be regularized by the sparsity constraint.

Grouping arrangement

To properly optimize the ℓ 2,1 norm regularized cost function, a key issue is to clarify the grouping arrangement. Three grouping choices are proposed. An example of wavelet quadtree structure is given in Fig. 4.8. Given that the elements surrounded by boxes of the same color belong to the same group, three groups appear:

• {1, 5} (in red) • {1, 2, 6} (in blue) • {1, 2, 3, 4, 5} (in yellow)
The three subgroups represent three ways to group the wavelet coefficients. For a wavelet coefficient A, the subgroups can be constituted as follows:

• {A, parent of A} (in red)

• {A, all ancestors of A} (in blue)

• {A, all children of A} (in yellow)

All three grouping choices form overlapping subgroups. Thus, the regularized cost function must be tackled with specific strategies. Herein, one adopts the replication approach proposed in [START_REF] Chen | The Benefit of Tree Sparsity in Accelerated MRI[END_REF]. By replicating the crossing elements, the aforementioned subgroups become

• {1 a , 5 a } (in red) • {1 b , 2 a , 6} (in blue) • {1 c , 2 b , 3, 4, 5 b } (in yellow)
The nodes (1 a , 1 b , 1 c ), (2 a , 2 b ), and (5 a , 5 b ) are replications of nodes 1, 2 and 5, respectively. When applying the group sparsity regularization, the replicated elements can be considered as independent from each other.

Optimization of the group sparsity penalized cost function

By replicating the crossing elements in all subgroups, the subgroups are decoupled. Therefore, the optimization of the cost function 4.30 can be treated as a non-overlapping group Lasso problem where the corresponding proximal operator has an analytic expression.

Denote f d as the normalized data error in (4.30) and f o the normalized object error. Enforcing the replication strategy, the previous cost function becomes

F W (J , β) = f d (J ) + f o (J , β) + Λ 1 ∥z g ∥ 2,1 + Λ 2 2 ∥z -G βz β∥ 2 2 (4.31)
z is a "flattened" version of β, where all overlapping groups in β become non-overlapping by replicating the crossed elements. Λ 1 = {λ 1,1 , λ 1,2 , . . . , λ 1,K } and Λ 2 = {λ 2,1 , λ 2,2 , . . . , λ 2,K } are two hyperparameter sets with K the total number of iterations. G βz is the mapping matrix that yields z from β such that z = G βz β (4.32)

Although the elements appearing repeatedly are replicated and become independent from each other, they are essentially the same coefficient and should be consistent. The second regularization term aims to maintain consistency between β and its replicated and flattened version z.

The analytic solution of ℓ 2,1 norm for non-overlapping group Lasso problems is introduced in equation (4.28). Since there is a second regularization norm in the cost function (4.31), the solution becomes

x ′ g i = x g i × 1 - λ 1 λ 2 ∥x g i ∥ 2 + , ∀i (4.33)
The optimization procedure is sketched in Algorithm 3. Update the contrast χ k by 10:

β k-0.5 = W{ Ns i=1 J i,k E tot * i,k Ns i=1 ∥E tot i,k ∥ 2 } (β k-0.5 means the preliminary result in k-th iteration) 11: Update z g k = (G βz β k-0.5 ) g × 1 - λ 1,k /λ 2,k ∥(G βz β k-0.5 ) g ∥ 2 + 12: Update β k = β k-0.5 -∇ β λ 2,k 2 ∥z -G βz β k-0.5 ∥ 2 2 13:
Update hyperparameters λ 1,k+1 and λ 2,k+1

14:

Reconstruct the spatial-domain contrast by χ k = W * β k 15:

Enforce positivity constraint on χ k 16: until a stopping criterion is satisfied.

Implementation details

Wavelet settings

At each iteration, the contrast χ is decomposed by the DT-CWT. The filter bank (FB) of the first stage of the DT-CWT must be different from those of the succeeding stages to make the DT-CWT approximately analytic for every stage. Let us denote l max the maximum length of the decomposition filters and N x and N y the row and column dimensions of the image, respectively. The level of decomposition L then should be a positive integer holding that:

L ≤ log 2 min (N x , N y ) l max + 1 (4.34)
Hereafter, Farras filters are used for the first stage, and Kingsbury Q-shift 6-tap filters for the following stages. So, here, N x = N y = 64, l max = 10 (due to the use of the Kingsbury Q-shift 6-tap filters) and according to (4.34) the largest possible decomposition level available is L = 3 which will be used in the following. The latter allows a finer description of the reconstructed obstacle thanks to the largest number of detail coefficients to be retrieved.

As seen in Fig. 4.8, three grouping strategies are adopted to introduce the structure information of the wavelet quadtree. Denote the three strategies as CSI-GS-1 (red), CSI-GS-2 (blue), and CSI-GS-3 (yellow). Here, the possible group sizes of three grouping arrangements are discussed. For a wavelet quadtree structure with a level of decomposition L:

• CSI-GS-1: If A belongs to the subbands {HL, LH, HH} L , then A has no parent coefficient, and there is only one element in the subgroup. Otherwise, the subgroup contains A and the parent of A. The group size of CSI-GS-1 is 1 or 2. q ℑ 85 95 68

• CSI-GS-2: If A belongs to the subbands {HL, LH, HH} L , as previously analyzed, the corresponding group size is 1. If A belongs to the subbands {HL, LH, HH} 1 , the group size is L. Therefore, the group size of CSI-GS-2 varies according to the level of decomposition, ranging from 1 to L. • CSI-GS-3: If A belongs to the subbands {HL, LH, HH} 1 , then A have no children coefficients and in this case the group size is 1. Otherwise, A and its four children are included in the same subgroup.

Consequently, the group size of CSI-GS-3 is 1 or 5.

Choice of hyperparameters

In equation (4.31), two hyperparameter sets Λ 1 and Λ 2 have to be determined. Following [START_REF] Taşkin | Sparsity Regularized Nonlinear Inversion for Microwave Imaging[END_REF], they are chosen adaptively according to the magnitudes of wavelet coefficients and updated at each iteration.

Let us denote P i% ( ) the i-th percentile of a vector. At k-th iteration, λ 1,k is set to P q% (β k ) and λ 2,k is associated with λ 1,k where the ratio of λ 2,k to λ 1,k is set to P q% (∥z g ∥ 2,1 ). q, which remains the only hyperparameter to be chosen, is fixed via a "trial-error" scheme according to the applied grouping strategy.

Since χ is a complex-valued function, q is also a complex-valued hyperparameter:

q = q ℜ + iq ℑ (4.35)
q has to be a priori fixed and the choice is made by running the inversion for a set of (q ℜ , q ℑ ) (varying from 50 to 95). In Fig. 4.9 the final error Err (χ) is plotted as a function of q ℜ and q ℑ and compared to the one of CSI (green line).

The above exhibits that CSI-GS methods outperform the CSI one for a large range of q ℜ and that their effectiveness does not depend too much on the chosen value. It can be noticed that the range for q ℜ and q ℑ is narrower for CSI- ) than for CSI-GS-1 (Fig 4.9(a)) and CSI- ). The choices of q ℜ and q ℑ for different methods are in Table 4.1 and are drawn as a point in each subfigure of Fig. 4.9.

Numerical results

In the numerical simulation, the model used is still the "Austria" profile. The settings of the scattering system are the same as in Chapter 3. To remind, the region of interest D is l = 2 m sided. For the direct problem, D is discretized into 50 × 50 cells, and for the inverse one, it is discretized into 64 × 64 cells. The region of interest is contained in free space C. N s sources and N r receivers are evenly distributed on an observation circle of radius r = 3.75 m. The relative permittivity of the embedding medium is ε b = 1. The collected data are corrupted by Gaussian noise. The frequency of operation is 600 MHz.

In this Chapter, the algorithm is tested on lossy scatterers. The scatterer has constant permittivity and conductivity. ε(r) = 2 and σ(r) =17 mS/m. Fig. 4.10 shows the relative permittivity and conductivity of the model. .9: Parametric study for the choice of q ℜ and q ℑ . The final error Err (χ) is plotted as a function of q ℜ and q ℑ and compared to the one of CSI (green line). N s × N r = 32 × 32 When there are 32 sources and 32 receivers, the real part of the model retrieved using CSI is roughly reconstructed. The general shape of the disks and rings can be recognized but with obvious discontinuities.

As for the imaginary part, only the external boundary of the profile is preserved with also discontinuities, whereas the internal pixels cannot be found.

With the group sparsity-guided regularization, CSI-GS-1, CSI-GS-2, and CSI-GS-3 allow better reconstructions. The model is smoother, and the edges are well preserved. The reconstructed permittivity and conductivity both have relatively smaller values. When the number of sensors is reduced to N s × N r = 18 × 26 and 16 × 24, more and more discontinuities appear in the results of CSI. Still, the proposed methods are capable of better reconstructing the model.

It should be noted that the results obtained by CSI-GS-3 are not as smooth as with the other two methods.

Rectangle artifacts appear in both reconstructed permittivity and conductivity. This is due to the grouping method. The additional group sparsity regularization can be seen as an image restoration step. The pixels are supposed to be restored from their ancestors or children at adjacent scales. The first two grouping methods connect each pixel to its ancestor(s), making it possible to smooth pixels with the same ancestor(s).

However, this implicit regularization disappears in the third grouping method as it places each pixel and its Fig. 4.12 shows the inversion curves of CSI, CSI-GS-1, CSI-GS-2 and CSI-GS-3, respectively. The data and object errors decrease smoothly for all four methods.

Analysis of robustness

To evaluate the robustness of proposed methods, noises of higher levels are considered with SNR = 15 dB and 10 dB. Simulations are conducted with CSI and the proposed methods with the same parameter settings. Relative errors are given in Table 4.3 andTable 4. When the noise is of 15 dB, the relative error of χ with CSI, CSI-GS-1, CSI-GS-2, and CSI-GS-3 for the configuration N r × N s = 32 × 32 increases by 21%, 2%, 3% and 2%, compared with a noise of 20 dB. For the configuration N r × N s = 18 × 26, the increase is of 14%, 3%, 3% and 2%. For the configuration N r × N s = 16 × 24, the increase is of 7%, 5%, 6% and 4%.

For a 10 dB noise, the increase is of 55%, 7%, 13% and 7% for the configuration N r × N s = 32 × 32. For the configuration N r × N s = 18 × 26, the increase is of 52%, 13%, 22% and 11%. For the configuration N r × N s = 16 × 24, the increase is of 21%, 9%, 13% and 8%.

To conclude, the quality of reconstruction degrades if higher and higher noise level for all algorithms.

Yet, the proposed methods are much more robust compared with CSI.

Comparison of DT-CW T and DWT filters

To further illustrate the advantages of the DT-CWT, experiments are conducted with DWT using the proposed methods.

Before that, the Kingsbury Q-shift filters used in our work have been explored. In previous work, the Kingsbury Q-shift 6-tap filters are used in the filter banks of DT-CWT. Here, the 14-tap and 18-tap filters are applied to the proposed methods to compare their performance. Kingsbury Q-shift 14-tap and 18-tap filters Relative errors of reconstructed contrast obtained with Kingsbury Q-shift 14-tap and 18-tap filters are given in Table 4.5 and Table 4.6. There is no noticeable improvement compared to the results obtained by 6-tap filters (Table 4.2).

Various wavelet filters Knowing that a large set of wavelets exists, only some more common ones are used. Relative errors obtained by CSI-GS-1, CSI-GS-2 and CSI-GS-3 using Haar, Daubechies 4 (db4), Daubechies 8 (db8), Coiflets 2 (C2), and Coiflets 4 (C4) are given in Table 4.7, 4.8 and 4.9. It can be observed that the relative errors obtained using DWT with different wavelet bases do not vary much. The results also show that the reconstruction quality obtained with DWT is not as good as that obtained with DT-CWT (Table 4.2).

Illustrative results obtained with db4 wavelet are shown in Fig. 4.15. Suffering from lack of directionality, only details of three basic directions (horizontal, vertical, and diagonal) are well retrieved. Also, the weaker persistence between parent and children coefficients reduces the retrieval quality.

Comparison of ℓ 2,1 norm and ℓ 1 norm

To enforce sparsity, the most popular regularization term is the ℓ 1 norm. To further demonstrate the benefit of the structure information provided by the ℓ 2,1 norm, simulations are led with ℓ 1 -penalized CSI, in which the group sparsity regularization term λ∥β g ∥ 2,1 in (4.30) is replaced by λ ′ ∥β∥ 1 .

To tackle the new cost functional with the ℓ 1 norm, lines 11-12 in the algorithm dedicated to the ℓ 2,1 norm are replaced by β k = S λ ′ {β k-0.5 } with S soft-thresholding operator. The hyperparameter λ ′ is determined like λ 2 . λ 1 = P q% (β), and q is chosen in accord with Table 4.1.

The results of ℓ 1 -penalized CSI are shown in Fig. 4.16. The permittivity shape can be roughly identified, though obvious discontinuities appear when the number of sensors is reduced. As for the conductivity, only edges are found. The relative error is in Table 4.10. To conclude, the performance of ℓ 2,1 norm is better than that of ℓ 1 norm thanks to the incorporation of structure information provided by the parentchild relationship between wavelet coefficients.

As for respective CPU times of the methods, they are illustrated on Model 1, with maximum iteration number of 500, on an Intel Core i5-8365U CPU (1.60 GHz) with 16 GByte memory, refer to Table 4.11.

Test with circular dataset

The proposed method is further tested on a synthetically generated circular dataset, whose configuration and implementation details are the same as in section 3.4.2. Here, CSI and CSI-GS-1 are further tested using the dataset and the reconstructed results of four samples from the circular dataset are given in Fig. 4.17.

The mean and variance of Err(χ) obtained by CSI and CSI-GS-1 on 250 samples are given in Table 4.12.

Conclusion

To conclude on the work discussed in the present chapter, the proposed group sparsity penalized CSI can improve the quality of reconstruction compared to CSI with reduced information. However, it can bring meaningful improvement only if CSI provides a satisfactory result as a good starting point.

Basics of convolutional neural networks

Deep learning, mostly convolutional neural networks (CNN), have achieved great success for tasks as computer vision, e.g., image segmentation [START_REF] Litjens | A Survey on Deep Learning in Medical Image Analysis[END_REF], object recognition [START_REF] Liang | Recurrent Convolutional Neural Network for Object Recognition[END_REF].

Deep learning is also successful for inverse imaging problems [START_REF] Mccann | Algorithm Unrolling: Interpretable, Efficient Deep Learning for Signal and Image Processing[END_REF][START_REF] Jin | Deep Convolutional Neural Network for Inverse Problems in Imaging[END_REF][START_REF] Lucas | Using Deep Neural Networks for Inverse Problems in Imaging: Beyond Analytical Methods[END_REF], like image denoising [START_REF] Zhang | Image Denoising Using a Neural Network Based Non-linear Filter in Wavelet Domain[END_REF], restoration [START_REF] Zhang | Learning Deep CNN Denoiser Prior for Image Restoration[END_REF], superresolution [START_REF] Cui | Deep Network Cascade for Image Super-resolution[END_REF], with worthy application in MRI-based medical imaging [START_REF] Aggarwal | Modl: Model-based Deep Learning Architecture for Inverse Problems[END_REF] as well.

Deep-learning-based methodologies have recently been much explored in the ISPs community, which will be further discussed in the next chapter. This chapter introduces some basics of CNN and the unrolled method, based on which we will construct our network structure.

Introduction

Deep learning [START_REF] Lecun | Deep Learning[END_REF]) is a set of learning methods that model data with complex network architectures and nonlinear transformations. In the last decade, several types of neural networks have been developed. The multilayer perceptrons, composed of multiple layers of perceptrons, are the oldest and simplest ones. Recurrent neural networks (RNN) are usually used to process sequential data such as text or time series signals.

Among these different network structures, a Convolutional Neural Network (CNN) is designed to process data with a grid pattern, such as images. Note that the data to be processed by a traditional artificial neural network (ANN) should be in vector form. Grid-like data such as images must be expanded into a onedimensional vector to be processed by ANN, which brings three limitations. Due to the high dimension of the resulting vector, the ANN tends to struggle with the computational complexity required to deal with image data. Moreover, the spatial information contained in the images is also lost when the images are expanded into vectors. The third limitation is the presence of overfitting phenomena due to a massive number of parameters contained in traditional ANN structures. CNN can alleviate these problems. CNN acts directly on 2-D matrices or 3-D tensors, with the third channel containing extra information, e.g., three RGB color channels for digital images, real and imaginary parts in our case. Besides, the convolution operation can extract the spatial features from the original images. Furthermore, fewer parameters are produced by CNN than by traditional ANN by replacing the fully-connected layers with convolutional layers. The fewer parameters required to train the network, the less likely the network will overfit, and the predictive performance of the model will be improved. A simple CNN architecture is illustrated in Fig. 5.1. The convolution operation is usually followed by a nonlinear activation function (usually a rectified linear unit (ReLU)), forming a convolution block together.

CNN layers

Pooling layers are used to downsample the feature maps to reduce their dimensionality. The above layers are the core part of CNN and are usually repeated several times to improve the performance of the model.

Finally, the fully-connected layers are typically used to interfere with all extracted features and produce the final output. More details on each layer are given below.

Convolutional layer

The convolutional layer is the core part of the CNN structure. Let us recall the definition of the convolution between two functions f and g

(f * g)(x) = t f (t)g(x -t) (5.1)
Likewise, the 2-D convolution for images is defined as

(K * I)(i, j) = m,n K(m, n)I(i -m, j -n) (5.2)
K is the convolution filter, also called the kernel. Consider an input sensor I, at each location of the tensor, the convolution calculates an element-wise product between each element of the kernel in the first step.

The second step is to sum up, the previous results to obtain the final output value in the corresponding location of the output tensor. Then, the kernel moves by a number s of pixels, s is called the stride. In other words, the kernel slides over the whole input tensor.

Zero-padding

The above convolution operation does not allow the center of each kernel to overlap the outermost element of the input tensor. Furthermore, each convolution operation reduces the height and width of the feature map. To effectively extract the boundary information of the image and control the dimension of the output feature map, we usually use padding techniques, especially zero-padding, where rows and columns of zeros are added on each side of the input tensor. We can then alter the dimension of the output of convolutional layers.

Suppose the input tensor is of dimension W i × H i × C i (W i and H i denote the width and height of the 2-D image, respectively. C i is the number of channels.) If we apply C o kernels, then the output tensor is of dimension W o × H o × C o , the output size after the convolution and padding operation can be calculated by the following formula.

W o = W i -k + 2p s + 1 (5.3) H o = H i -k + 2p s + 1 (5.4)
where k is the kernel size, p represents the amount of zero padding set, and s is the stride.

Stride and kernel size The stride and kernel size are two key parameters to characterize the convolution operation. Define the dimensionality of the region overlapping with the kernel as the receptive field size of the neuron.

• Stride The stride of the kernel represents the extraction accuracy of the features. If the stride is small, more features will be extracted, and less information will be lost. Yet, it may result in a higher computational burden because more parameters are left after the convolution operation. In addition, the choice of the stride also affects the spatial size of the output feature maps.

• Kernel size The kernel size represents the size of the receptive field, and it is generally odd number× odd number i.e., 1 × 1, 3 × 3, 5 × 5, 7 × 7, for the following reasons:

-Benefit the padding operations As previously mentioned, we sometimes need the size before and after convolution to remain unchanged. Then, we need to use padding techniques. Following equations (5.3) and (5.4), the input and output sizes are the same given that the padding amplitude is set to (k -1) 2 while s is set to 1. It is one possible choice and the most common choice to ensure that the spatial dimension remains unchanged before and after the convolution. If k is a even number, (k -1) 2

is not an integer.

-Facilitate finding "anchor" point In CNN, the convolution operation is generally performed based on the position of the kernel. The "anchor" point is used to determine the position of the kernel with respect to the image. If the kernel size is an odd number, the "anchor" point is easy to find, which is naturally the center of the kernel. However, if the convolution kernel is an even number, the "anchor" point is difficult to find.

Nonlinear activation function

The outputs of convolutional layers are then passed through a nonlinear activation function. Generally, nonlinear activation layers are employed after all learnable layers. The main role of the activation function is to provide the neural network with nonlinear modeling capabilities. Convolution is a linear operation, so non-linear elements must be introduced to solve problems that linear models cannot. The neural network has hierarchical nonlinear mapping capability only after adding the activation function.

Activation functions (especially ReLU) can also build sparse matrices, which can remove redundancy in the data and preserve the characteristics of the data as much as possible. The ReLU is a simple thresholding operator ReLU(x) = max(x, 0)

(5.5)

Pooling layer

The pooling layer is usually placed after the convolution block. Pooling is essentially a down-sampling operation. For a given input feature map, the pooling operation takes the maximum or mean value on a small region. Similarly to the convolutional layer, we should choose a stride value and the small patches on which the pooling acts.

Among the pooling operators, max-pooling is the most commonly used. Fig. 5.3 illustrates an example of a max-pooling operation. The pooling process is similar to the convolution process. As shown in Fig. 5.3, the value in the neighborhood of a feature map is "scanned" with a 2 × 2 filter with a step size of 2. Here, the maximum value in the neighborhood is taken and output to the next layer. Unlike convolutional layers, the pooling layer has no parameters to update during the training process. It has two main roles:

• To retain the main features while reducing the dimensionality of feature maps to prevent overfitting and reduce the computational burden.

• To make the network less sensitive to small translations of the input tensor.

Fully-connected layer

After several repeated convolutional blocks and pooling layers, a traditional CNN structure generally ends with several fully-connected layers to map the extracted features to the final output. Namely, the neurons of a fully-connected layer are directly connected to the neurons in adjacent layers. This is why the dimension of the input tensor of fully-connected layers should be fixed.

The parameters of the fully connected layer are redundant, thus increasing the computational burden.

This also leads to overfitting of the network. Therefore, it is common to replace fully connected layers with global average pooling layers or convolutional layers with kernel size 1 × 1.

Batch normalization

The training of a neural network is difficult because CNN contains many hidden layers, and the parameters of each layer will change with training. Therefore, the input distribution of the hidden layer will always vary, and each hidden layer will face the problem of covariate shift. This phenomenon is called internal covariate shift, which makes the input of each layer no longer independent and identically distributed. Due to this phenomenon, the data input into the activation function will fall into the saturation region, reducing the learning efficiency and even leading to the vanishing gradient problem.

Batch normalization (BN) [START_REF] Ioffe | Batch Normalization: Accelerating Deep Network Training by Reducing Internal Covariate Shift[END_REF] can be regarded as a data pre-processing step at every layer of the network. BN is to normalize the input of each layer by mini-batch into a standard normal distribution with a mean of 0 and a variance of 1, thereby avoiding the vanishing gradient problem and greatly accelerating the convergence rate.

Training of CNN

In the previous part, different layers of CNN are described. Here, we discuss the training process of CNN by detailing several key concepts: loss function ( § 5.3.1), parameter initialization ( § 5.3.2), and optimizer choice ( § 5.3.3).

Loss function

Training a neural network is basically to find the parameters in the network which minimize the discrepancy between the output of the network and the ground truth in a training dataset. A loss function, also referred to as a cost function, is used to measure the data misfit between estimated results and the true one.

Different loss functions exist for different types of tasks. In ISPs, the CNN is usually used to obtain the contrast function (or contrast sources) from given inputs (scattered fields, pre-initialized contrast sources, etc.). It can be regarded as a regression task. Therefore, the Euclidean loss (or mean squared error (MSE))

is used.

The formula for calculating MSE loss between the output of the network ŷ ∈ R N and the ground truth y ∈ R N is as follows:

L(ŷ, y) = 1 N (ŷ -y) 2
(5.6)

The loss function value is calculated through forward propagation, where the input data is fed in the forward direction through the network and gets the output. Once we have the loss function value, this error will back-propagate from the last layer to the first layer based on the chain rule, and the gradient of each learnable parameter in the network can be obtained. These gradients will then be fed back to optimization methods (such as stochastic gradient descent), which are used to update the learnable parameters (kernels in convolutional layers and weights in fully-connected layers).

Parameter initialization

The initialization of the network parameters is essential because an inadequate initialization can result in a low convergence rate or even failure to converge. In this work, Kaiming initialization [START_REF] He | Delving Deep into Rectifiers: Surpassing Human-level Performance on Imagenet Classification[END_REF] is used, which is a recently proposed initialization method taking into account the non-linearity of activation functions, such as ReLU activations.

Optimization algorithm

The gradient descent method and its variants are commonly used in traditional optimization problems. It is also the usual selection for CNN learning problems. In our work, the Adaptive Moment Estimation (Adam) method is used, which is one of the latest trends in deep learning optimization. Adam has two key benefits:

it is more memory efficient than other algorithms and needs less computational power. The algorithm is summarized in Algorithm 4.

During the learning process, the learning rate, which is the step size of the parameter updating, should be pre-defined. It should be adequately chosen to train the network effectively. Learning rate decay is a mod-Algorithm 4 ADAM optimization method 1: Initialize the biased first moment estimate and biased second raw moment estimate m t = 0, v t = 0; 2: Initialize the hyper-parameter β 1,0 = 0.9, β 2,0 = 0.999, δ = 10 -8 , β ern training technique that starts the training with a relatively high learning rate and then decreases/decays it slowly. It has been empirically observed to benefit the learning process and improve the generalization ability of the trained network.

Due to the limitation of available GPU memory, the training data is usually divided into mini-batch data.

The number of samples in a batch is called batch size, a hyperparameter to be chosen. By employing the mini-batch training, the update of network parameters is not obtained by a single sample but is a weighted result from a batch of data.

The number of epochs, which defines the number of times the training algorithm will run on the entire training dataset, is also a hyperparameter to be pre-defined. It is a critical issue. Too few training epochs will lead to insufficient complexity of the trained network. Therefore, the network can capture the underlying data's features. This is called "underfitting". "Overfitting" is the opposite of "underfitting". Too many epochs lead to a too complex trained model that can capture even the noise in the data, and the training error of the model on the training data set will gradually decrease. However, the model does not work on data sets other than the training set.

The early stopping method is usually employed to address the problem of "overfitting". Early stopping consists in truncating the number of epochs in advance, that is, stopping epochs before the model con- With the components of CNN introduced above, we can construct and train a neural network. In the next section, the unrolled method is introduced, on which our network is designed.

Unrolled method

Background

Deep learning has achieved great success in some problems for computer vision tasks. In high-level vision problems such as face recognition and target detection, the performance of deep networks has been far superior to traditional methods. However, in some low-level vision problems, such as inverse imaging problems, data-driven methods cannot completely replace traditional model-driven methods. networks has become a new research trend in the field of deep learning. Among them, the unrolled method (or unrolling method) [START_REF] Mccann | Algorithm Unrolling: Interpretable, Efficient Deep Learning for Signal and Image Processing[END_REF] provides effective means to address this issue.

Model-driven methods consist in establishing a mathematical model for the problem based on domain knowledge and then employing optimization tools to obtain the optimal solution of the model. Correspondingly, data-driven methods directly establish a nonlinear mapping relationship between input data and output results, and their performance is often closely related to the quality of training data and the network structure. To sum up, we find that the above two approaches mainly focus on using "domain knowledge" or are based on "data (and experience)", respectively. Then, if we can effectively combine these two different strategies, it will be very promising to design more robust and interpretable deep network models.

Designing an unrolled algorithm consists of two main steps:

• Take a traditional iterative algorithm and unroll it.

• Make some elements of the algorithm learnable.

In the following, we will give some examples of unrolled algorithms for solving inverse imaging problems.

Unrolled method for inverse imaging problems

The general framework and variants of the unrolled method have been much discussed to solve the inverse problems of interest. It replaces the components of iterative optimizers with learnable parameters or networks while integrating the knowledge of the physical model. Unlike the commonly used deep learning schemes that mainly adopt a single network to map the measurements or the approximate reconstructed results to the ground truth, the unrolled method uses a recurrent block architecture where the number of repeated blocks corresponds to the number of iterations.

In [START_REF] Ongie | Deep Learning Techniques for Inverse Problems in Imaging[END_REF], learning methods for solving different inverse imaging problems are discussed according to whether or not a forward model is known and to what extent it is used in training and testing, and whether or not the learning is supervised or unsupervised. The supervised case where the forward model is fully known is presented here. A general unrolled scheme is provided, based on the proximal gradient descent algorithm, where a trainable CNN replaces the proximal operator. In addition to the proximal operator, other free parameters, such as the step-size of descent, are also transformed into learnable parameters. Denote K the forward model and K-1 an approximate inverse of K, and we have

x = K-1 (y)

(5.7)

When applying the learning approach to tackle inverse problems, the goal is to estimate a network that maps measurements y to images x. Suppose K is a fully known linear operator; it can be embedded into the network structure in two ways. First, the incorporation of knowledge of K can be realized by mapping the measurements back to the image domain using the approximate inverse operator K-1 . The result of the previous step is then used as the initialization of the network. In addition to the initialization step, the forward model K can also be incorporated into multiple layers of an unrolled network.

Consider an optimization problem of the following form arg min

x 1 2 ∥K(x) -y∥ 2 2 + H(x) (5.8)
where H(x) represents the regularization term. The proximal gradient descent method is commonly used to optimize it. Recall the formulas of proximal algorithm (4.26) and (4.27) introduced in the previous chapter,

x can be updated at iteration k according to

x k = prox tH( ) x k -t k ∇ x 1 2 ∥K(x) -y∥ 2 2
(5.9) wherein t k the descent step size at k-th (5.10) Now, we unroll the above proximal gradient descent steps and take each iterative step as a sub-block for the network. The recurrent network architecture is sketched in Fig. 5.5. An initialization layer is placed at the beginning of the network by applying an approximate inverse operator to map measurements to the image domain. Then, the gradient descent step is performed. At the end of each sub-block, the proximal operator prox tH( ) ( ) is replaced with a trainable neural network mapping from images to images. The unrolled reconstruction network can be interpreted as learning a proximal operator.

In [START_REF] Mardani | Neural Proximal Gradient Descent for Compressive Imaging[END_REF] and [START_REF] Diamond | Unrolled Optimization with Deep Priors[END_REF], authors both use the residual network to model the proximal map in traditional iterative algorithms. In [START_REF] Zhang | ISTA-net: Interpretable Optimization-inspired Deep Network for Image Compressive Sensing[END_REF], the architecture is designed based on the Iterative shrinkage-thresholding algorithm (ISTA), which aims at solving sparse linear inverse problems with ℓ 1 -norm. Instead of using the traditional sparsifying transform, such as wavelet transform, two symmetric convolutional blocks are applied to construct a learnable non-linear sparsifying transform and its inverse. In [START_REF] Jin | Deep Convolutional Neural Network for Inverse Problems in Imaging[END_REF], the relationship between the general CNN structure and the iterative optimization methods is explored. Given a forward model where the forward operator is the convolution, the process of the corresponding proximal gradient descent algorithm can be interpreted as iteratively applying the convolutional filters, adding a bias, and using a point-wise nonlinearity, which has a similar form to the CNN. Based on this connection, modifications have been made to the widely used U-net architecture [START_REF] Ronneberger | U-net: Convolutional Networks for Biomedical Image Segmentation[END_REF] to reconstruct X-ray CT image. Likewise, in [START_REF] Li | Deepnis: Deep Neural Network for Nonlinear Electromagnetic Inverse Scattering[END_REF], the authors investigate the connection between iterative methods for nonlinear ISPs and the deep neural network (DNN).

A complex-valued DNN is designed based on this connection.

Unrolled convolutional neural networks for ISPs

In this chapter, an unrolled deep learning scheme for solving full-wave non-linear ISPs is presented. An iterative neural network structure combining the CSI method and the residual network (ResNet) is designed, which is inspired by the so-called unrolled method. In the following, the application of CNN in ISPs is first discussed. Then, the proposed unrolled scheme combined with the CSI method is presented. Thorough numerical tests are carried out to evaluate the performance, stability, robustness, and reliability of the proposed approach.

Introduction

Recently, CNN have shown great potential in solving ISPs [START_REF] Massa | DNNs As Applied to Electromagnetics, Antennas, and Propagation-a Review[END_REF][START_REF] Chen | A Review of Deep-learning Approaches for Inverse Scattering Problems (invited Review)[END_REF][START_REF] Salucci | Artificial Intelligence: New Frontiers in Real Time Inverse Scattering and Electromagnetic Imaging[END_REF]. Unlike traditional iterative methods that get the solutions from measurements and the domain knowledge, they benefit from the large data set to learn the unknown solution while the domain knowledge is less involved.

The advantages of applying deep learning tools to solve ISP are twofold. On the one hand, it allows realtime reconstruction whereas traditional iterative methods require computational time. On the other hand, by benefiting from the information contained in large-scale data, deep learning tools can sometimes solve highly non-linear ISPs for which traditional methods fail. However, collecting large-scale datasets may not be straightforward in real-world applications.

Deep learning approaches for solving ISPs can be roughly categorized into three types as follows, as already briefly introduced in Chapter 1.

• The contrast distribution can be directly derived from the measured data, the input of the neural network being the scattered fields, and the output being the contrast [START_REF] Ran | Subwavelength Microstructure Probing by Binary-Specialized Methods: Contrast Source and Convolutional Neural Networks[END_REF]. However, learning contrasts directly from scattered fields takes unnecessary time to learn known physical processes. Therefore, we need to combine domain knowledge and learning methods.

• Approaches of the second type consist in incorporating physical knowledge into deep learning (Wei et al. 2019a;Wei et al. 2019b;[START_REF] Liu | Physics-guided Loss Functions Improve Deep Learning Performance in Inverse Scattering[END_REF]). These methods can somewhat reduce the nonlinearity of the relationship between the input and output of the neural network and simplify learning.

• A neural network can be combined also with traditional optimization methods. This third type still employs the traditional iterative approach to solve the ISPs but uses deep learning to complete the optimization procedure [START_REF] Sanghvi | Embedding Deep Learning in Inverse Scattering Problems[END_REF].

To achieve sound retrievals, incorporating domain knowledge into the deep model is of great importance.

In addition to initializing the network with approximate results based on physical information as often made, the unrolled method, which has been introduced in the previous chapter, is a promising tool to address this issue, as exemplified in the recent contribution [START_REF] Liu | Physical Model-inspired Deep Unrolling Network for Solving Nonlinear Inverse Scattering Problems[END_REF].

𝐸 diff , 𝐸 inc , 𝐺 𝑆 , 𝐺 𝐷

Proposed unrolled CNN scheme

In this section, a novel unrolled model which combines updating steps of CSI and a classical CNN structure is detailed. Denote the proposed unrolled network as CSI-Net.

The entire network has a recurrent block architecture that contains a series of sub-blocks. Each includes a set of CSI update steps and a CNN block, which corresponds to an update of contrast and contrast source, as it will be seen.

Then, several strategies are designed to improve the model. Due to the recurrent network structure, a combined loss function is set up, the weight parameters being determined automatically during training.

Modifications are also made compared to a general unrolled method to adapt the network to nonlinear inverse problems.

The overall structure of the proposed network is sketched in Fig. 6.1, where each sub-block of the network contains two modules: Module CSI and Module CNN. In the following, the CSI block is first reviewed as the main brick of Module CSI. Then, the network structure adopted for Module CNN is introduced. Details of the approach follow.

Module CSI

Module CSI is constructed based on the CSI method, which was detailed in Section 3.1. As shown in Fig. 6.1, each sub-block contains a series of Module CSI, and each Module CSI contains one single CSI iteration.

In the first sub-block, Module CSI aims at providing an immediate initialization of J and χ for the Module CNN. Then, in the following sub-blocks, Module CSI further refines the results that come from the precedent Module CNN. Note that the step size α J is not a learnable parameter and it is determined by minimizing both data and state error functions.

Module CNN

In the proposed unrolled schema, the network structure in Module CNN is stacked by residual blocks.

Residual learning has been proposed in [START_REF] He | Deep Residual Learning for Image Recognition[END_REF] to deal with network degradation, which means that adding more layers might degrade the learning performance. The entire network structure of Module CNN is illustrated in Fig. 6.2. It thus consists of two residual blocks and two convolutional blocks.

Adding a shortcut between the input x and the output, the network has to learn the residual H(x) between ground truth and input. It should be easier to optimize than the original mapping, and this shortcut ensures that higher layers perform at least as well as lower ones.

The input of Module CNN is the contrast χ obtained by Module CSI, which is a tensor of shape

[Number of samples × Number of channels × H × W ]
where H and W represent image height and width.

The number of input/output channels is 2 for a complex-valued image by splitting real and imaginary parts into two different channels; for a real-valued image, the second channel is empty.

Conv. For all convolutional layers, the kernels used to obtain the feature maps are of size 3 × 3. Inside the residual block, the number of kernels for each layer is 64. The number of kernels of layers Conv.1 and Conv.2 in Module CNN is 64 and 2, respectively.

To keep the spatial size of feature maps unchanged after each convolution, the stride of the convolution operation is set to 1, and a 1-sized zero-padding is added. Each convolutional layer is followed by batch normalization (BN).

Finally, a ReLU layer is added as the nonlinear activation layer used to overcome the vanishing gradient problem and add the nonlinearity to the model.

CSI-based unrolled CNN framework: CSI-Net

Inputs of CSI-Net As depicted in Fig. 6.1, the network contains a series of sub-blocks. For the first sub-block, the input is the pre-initialized contrast χ and contrast source J , obtained by standard backpropagation.

A closed-form contrast source is obtained as:

J 0 = ∥G * S E diff ∥ 2 ∥G S G * S E diff ∥ 2 G * S E diff (6.1)
and the corresponding contrast χ 0 derived from J 0 :

χ 0 = Ns is=1 J is,0 E tot * is,0 Ns is=1 ∥E tot is,0 ∥ 2 (6.2)
Additionally, the inputs include as well auxiliary data E diff , E inc , G S , G D to perform the CSI iterations.

Forward path Denote N p the total number of sub-blocks and N k the total number of CSI iterations performed inside each sub-block. In the p-th sub-block, there are two main steps.

• Firstly, the input passes through a sequence of Module CSI. In the commonly used unrolled deep learning framework, the gradient descent update step is performed only once inside each sub-block, where the step size is considered as a free parameter learned during training. In contrast, the step size α J here is a constant determined at each iteration, thus not a learnable parameter. Moreover, instead of running a single CSI iteration in each sub-block, the CSI iteration is repeated several times.

This modification comes from the following considerations: Data in Module CSI are not normalized, so the initialization of step size is essential and should respect the correct scale. So it is more convenient to determine the step size by following a traditional optimization. As it cannot be learned, the single iteration step is replaced by multiple iterations to speed up convergence and amplify improvements brought by Module CSI. At the end of this step, a preliminary reconstruction χ p-0.5 of the contrast χ is obtained.

• Secondly, the previous result χ p-0.5 is applied as input of Module CNN. Note that it should be normalized before being involved in the network. The output of this module is χ p , which is the final result of the p-th sub-block, and the corresponding contrast source J p is updated by solving the forward problem.

Loss function

λ i = L 0 L i (6.5)
where N e denotes the number of elements in xi , and λ i is the weighting coefficient of i-th sub-block.

Evaluation's criteria & implementation's details

Evaluation criteria To evaluate the reconstruction quality from different points of view, in addition to the normalized reconstructed contrast error defined in equation (3.30), two more criteria are adopted: a proposed network. To this end, a real-valued circular dataset including 5000 samples is used. Five datasets of 2000 data pairs are randomly sampled from the latter as training datasets for the models. Fig. 6.9 illustrates that the performance of the proposed approach remains stable when the dataset changes.

Generalization capability is always an issue, so challenging examples must be considered. The models trained with the circular dataset are now tested on six test samples that are not within the range of the training dataset. The results are in Fig. 6.10 and 6.11. Notice that positivity constraint is enforced in CSI, significantly improving its performance. The number of iterations is set to 200. In particular, samples #7 and #8 are the well-known "Austria" profile with relative permittivity of 1.5 and 2, respectively.

CSI-Net achieves the best reconstruction for most of the samples from #7 to #12. For sample #7, U-net and CC-ResNet only give a rough reconstruction, where the center of the ring is obscured, and the shape is not continuous. For sample #8, where the level of nonlinearity increases, it fails. As for CSI, the retrieved profiles have a smaller size than the true ones. Also, artifacts show up around the edges. For the other samples, U-net and CC-ResNet fail to provide meaningful results and CSI can achieve good results, though with more artifacts than CSI-Net results. The Err χ of the results are in Table 6.2. They show that CSI-Net has a better generalization capability than U-net and can improve the reconstruction quality of several cases compared to CSI.

Test with MNIST dataset

The proposed model is trained with a different dataset where the scatterers are modeled from the MNIST database [START_REF] Deng | The MNIST Database of Handwritten Digit Images for Machine Learning Research[END_REF], which contains handwritten digit samples. Digits are transformed to scatterers of random relative permittivity between 2 and 2.5. To increase diversity, a circular/rectangle of random relative permittivity between 2 and 2.5 is added to each sample. The diameters/side lengths are in between 0.3 m and 0.8 m. multi-frequency data, we can observe that the CSI-Net-MF outperforms CSI-Net in the monochromatic case.

Relative permittivity errors of six samples and the mean value over 200 test samples shown in Table 6.5 confirm it.

Conclusion

To conclude, the proposed CSI-Net can combine the advantages of traditional iterative algorithms and learning approaches. It has a good generalization ability. Moreover, it can improve the confidence level of the results by integrating the domain knowledge into the learning process. The results obtained by the proposed model are more reliable that a black box network. Last but not least, due to the recurrent network structure of the proposed model, it is naturally suitable to be extended to an multi-frequency scheme.

Conclusion and perspectives

Conclusion

In this thesis, two-dimensional full-wave nonlinear electromagnetic inverse scattering problems (ISPs) are investigated. The aim is the retrieval of the characteristics of an unknown scatterer or set of scatterers illuminated by given incident fields from the knowledge of the fields it or they scatter. The proposed methods can mainly be classified into two parts: the traditional optimization methods and the learningbased method. For traditional methods, the discrete wavelet transform and mixed sparse regularization are combined with the contrast source inversion (CSI) method, while the learning-based method consists in combining convolutional neural networks (CNN) and the CSI method through a so-called unrolled scheme.

In the first part, the CSI method is combined with the conventional 2-D discrete wavelet transform. One attractive property of the wavelet transform is that it can provide a multiresolution representation for a given signal. By transforming the contrast function into the wavelet domain with a given wavelet basis, wavelet coefficients of different levels can be obtained, which correspond to different frequency resolutions. Wavelet coefficients can be divided into two parts: approximation coefficients and detail coefficients.

Approximation coefficients contain low-frequency components and thus are more robust to noise. A multiscale wavelet-based CSI is proposed, combining the multiresolution property and the characteristics of approximation wavelet coefficients. Simulations are carried out on different configurations and with diverse wavelet bases to compare the performance of CSI and of the proposed method.

In addition to the multiresolution property of the wavelet coefficients, the sparsity of wavelet detail coefficients also attracts considerable attention. So, to alleviate the ill-posedness of the ISPs, a sparsity penalization is added to the cost function of the CSI method, and a group sparsity penalized CSI method is proposed. Instead of employing the widely-used ℓ 1 norm to impose the sparsity constraint, here, a mixed norm ℓ 2,1 norm is used to impose the so-called group sparsity.

The conventional discrete wavelet transform used in the first method is replaced with the dual-tree complex wavelet transform, which makes it possible to represent the contrast function in a sparse manner and yields the structure information of the contrast. When solving the optimization problem with group sparsity regularization, three grouping strategies are investigated to build overlapping sub-groups, while the replication approach and proximal gradient descent algorithm are used to minimize the cost function. Numerical simulations are carried out on different configurations, and comparisons are conducted to exhibit the pros and cons of the proposed methodology.

In the second part, the deep learning technique is applied to solve ISPs. Deep learning has achieved great success in dealing with computer vision tasks but has not produced a profound impact on the solutions of such ISPS yet. In our work, CNN is combined with the CSI method, and an unrolled network structure is designed. The domain knowledge is incorporated in an effective manner into the learning process by embedding CSI iterations into the network structure. Furthermore, due to the iterative architecture of the network, the proposed model can be easily extended to the multi-frequency scheme. Thorough numerical tests are carried out to evaluate the performance of our model, and the stability, robustness, and reliability of the proposed approach are also analyzed. Comparisons with the widely used U-net structure and CSI illustrate the advantage of the proposed method.

To summarize, one has developed and investigated both traditional optimization methods and learningbased methods for solving 2-D full wavelet ISPs. The multi-scale wavelet-based CSI and group sparsity penalized CSI methods show their advantages when tackling problems with insufficient collected data. As for the learning-based method, it enables real-time reconstruction and it benefits from the information contained in large-scale data. By combining CSI and CNN, the reliability of the model is improved, something which is still a matter of concern for learning-based methods.

Perspectives

In this thesis, several methods including traditional optimization approaches and learning approaches have been investigated. However, further explorations appear of good interest.

Traditional optimization methods

• In the multi-scale wavelet-based CSI method, the inversion starts from the finest decomposition level to the coarsest decomposition level. That is, several successive inversion steps are required throughout the optimization process. The present work sets a fixed number of iterations for each inversion step. The optimization involves only a few coefficients for the first inversion steps, so fewer iterations are required to converge. Therefore, the number of iterations needed for each inversion step is different. Further research should be done to define more flexible stopping criteria to eliminate unnecessary iterations.

Haar and Daubechies' wavelet family are chosen to decompose the unknown object in the present work. However, the optimal choice of the wavelet basis depends on different problems. More research can be carried out to deduce which wavelet basis function can be the most suitable for specific application areas.

• In the group sparsity penalized CSI method, three grouping strategies are proposed to form different sub-groups of wavelet coefficients, based on which group sparsity penalization is applied. Instead of using pre-defined group arrangements, the hidden Markov tree (HMT) model can be employed to dynamically detect the wavelet support at iteration [START_REF] Deka | Wavelet Tree Support Detection for Compressed Sensing MRI Reconstruction[END_REF].

This method focuses on the structural information contained in the wavelet detail coefficients. Images with more detail information should benefit better from this approach. Therefore, it would make sense to apply this approach to other models, such as breast models in biomedical imaging.

• Multi-scale wavelet-based CSI method and group sparsity penalized CSI method concentrate on approximation wavelet coefficients and detail wavelet coefficients, respectively. The combination of the two approaches is of good interest.

Learning method

• As for what concerns the CSI-Net, the step size of the CSI algorithm is not a learnable parameter.

Therefore, several CSI iterations are embedded into each sub-block to speed up the convergence.

Making the step size a learnable parameter is a path ahead, since it makes it possible to simplify the network structure.

• As emphasized already, incorporating physical information into the deep learning model is of great importance. Physical information can be integrated into the model not only through pre-initialized inputs and network structures but also through modifications to the loss function. The loss function used in the current work can be further modified by adding data misfitting terms obtained from domain knowledge, as proposed in [START_REF] Liu | Physics-guided Loss Functions Improve Deep Learning Performance in Inverse Scattering[END_REF].

Other potential work

• All proposed methods are still limited to 2-D configurations and in the standard scalar TM polarization case. Investigations of vector 3-D cases where mathematical modeling is more complicated, and nonlinearity and ill-posedness become more severe, should be carried out.

• The proposed methods are only tested on simulation data. More tests should be carried out on experimental data at first in laboratory-controlled situations.

A Orthogonal Daubechies coefficients 5.100 97 × 10 -3 5.002 268 53 × 10 -4 -5.540 045 49 × 10 -4 -6.054 960 58 × 10 -3 -0.015 00 9.552 297 11 × 10 -4 2.612 967 28 × 10 -3 1.973 36 × 10 -3 -1.661 372 61 × 10 -4 3.258 146 71 × 10 -4 2.817 59 × 10 -3 -3.563 297 59 × 10 -4 -9.699 40 × 10 -4 5.564 551 4 × 10 -5 -1.647 06 × 10 -4 1.323 67 × 10 -4 -1.875 × 10 -5

Table A.1

Titre : Imagerie électromagnétique non linéaire : d'algorithmes basés sur les ondelettes préservant la parcimonie à l'apprentissage en profondeur Mots clés : inversion de source de contraste, transformée d'ondelettes complexes à double arborescence, parcimonie de groupe, apprentissage profond, méthode déroulée Résumé : Ce travail traite de l'imagerie électromagnétique mal-posée non-linéaire à partir de données en régime harmonique dans une expérience de diffraction bidimensionnelle, en mettant l'accent sur deux approches dans le cadre de l'inversion de source de contraste (CSI). La première approche est un CSI pénalisé par la parcimonie de groupe dans le domaine des ondelettes, la seconde est un schéma d'apprentissage profond dit déroulé. Dans la première approche, une dépendance existe entre les coefficients d'ondelettes à différentes échelles, appelée relation parent-enfant, ce qui donne une structure d'arbre quaternaire d'ondelettes, de sorte que les coefficients d'ondelettes sont parcimonieux tant au niveau des pixels qu'au niveau des groupes. L'accent est mis sur l'utilisation de la transformée en ondelettes complexe (CWT) à double arborescence en vue d'obtenir correctement la représentation parcimonieuse par groupe recherchée. Une norme l2,1 ajoutée à la fonction coût standard applique la parcimonie de groupe aux coefficients d'ondelettes du contraste. La stratégie de réplication est combinée avec la méthode proximale afin de résoudre le problème pénalisé de groupes qui se chevauchent. La seconde approche est inspirée de la méthode déroulée. En intégrant les itérations CSI dans le modèle d'apprentissage profond, la connaissance du domaine est intégrée au processus d'apprentissage. Dans les deux cas, des tests numériques approfondis sont effectués afin d'évaluer les performances, la stabilité, la robustesse et la fiabilité, en menant des comparaisons avec des solutions considérées plus standard (comme CSI, la transformée en ondelettes discrète ou DWT, et U-net), et celles-ci illustrent l'avantage des deux approches proposées en de nombreux aspects.

Title: Non-linear electromagnetic imaging: from sparsity-preserving wavelet-based algorithms to deep learning Keywords: contrast source inversion, dual-tree complex wavelet transform, group sparsity, deep learning, unrolled method Abstract:This work deals with nonlinear ill-posed electromagnetic imaging from time-harmonic data within a two-dimensional scattering experiment, the focus being on two approaches in the framework of the contrastsource inversion (CSI). The first approach is a group sparsity penalized CSI in the wavelet domain, the second is an unrolled deep learning scheme. In the first approach, dependency exists between wavelet coefficients at different scales, referred to as parent-child relationship, which yields a wavelet quadtree structure so that wavelet coefficients are both pixel-wise and group-wise sparse. Emphasis is on the dual-tree complex wavelet transform (CWT) to achieve this result. A '2;1 norm added to the standard cost function is to enforce group sparsity onto the wavelet coefficients of the spatially-varying contrast. The replication strategy is combined with the proximal method to solve the overlapping group penalized problem. The second approach is inspired by the unrolled method. By embedding the CSI iterations into the deep learning model, the domain knowledge is incorporated into the learning process. In both cases, thorough numerical tests are carried out to evaluate performance, stability, robustness, and reliability with comparisons with more standard solutions (like CSI, discrete wavelet transform (DWT) and U-net), which exhibit the advantage of the proposed approaches under many aspects.

3

  Contrast source inversion and wavelet technique3.1 Decompose a m × n image X by 2-D DWT of level 1. . . . . . . . . . . . . . . . . . . . . . . 3.2 Decompose a m × n image X by 2-D DWT of level 3. . . . . . . . . . . . . . . . . . . . . . . 3.3 Wavelet coefficients generated by 2-D DWT. Left: 1-level decomposition. Right: 3-level decomposition. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Localization of wavelet coefficients employed in the optimization steps. . . . . . . . . . . . . 3.5 Left: Austria model, relative permittivity ε(r) = 2. Right: scattering system. . . . . . . . . . . 3.6 Normalized contrast error as a function of N s and N r (arbitrary threshold applied at 50%). (a) CSI,(b), (c) and (d) wavelet-based CSI with Haar, db4 and db14 wavelets respectively. Gaussian noise of 20 dB is added to the scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Normalized contrast error as a function of N s and N r (arbitrary threshold applied at 50%). (a) CSI,(b), (c) and (d) wavelet-based CSI with Haar, db4 and db14 wavelets respectively. Gaussian noise of 10 dB is added to the scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . 3.8 Inversion error curves of (a) CSI and (b) WCSI with db14 wavelets. Gaussian noise of 20 dB is added to the scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Reconstructed ϵ(r) of "Austria profile" for various combinations of receiver and source numbers. From top to bottom: results of CSI, WCSI with Haar, db4, db14 basis, respectively. Gaussian noise of 20 dB is added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 Reconstructed ϵ(r) of "Austria profile" for various combinations of receiver and source numbers. From top to bottom: results of CSI, WCSI with Haar, db4, db14 basis, respectively. Gaussian noise of 10 dB is added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11 Reconstructed ϵ(r) of 8 samples in circular dataset for N s × N r = 6 × 18. Gaussian noise of 20 dB is added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . is the filter used for the first stage, and h 0 (n) is for the succeeding stages. The implementation of the DT-CWT requires that the first stage of the filter banks be different from the succeeding stages. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 db4 wavelets (LH, HL, HH) associated with the 2-D DWT in the space domain. Level of decomposition is 5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Typical wavelets associated with the oriented 2-D DT-CWT in the space domain. . . . . . . . 4.4 (a): A "circle" example in spatial domain. (b): Wavelet coefficients obtained by 2-D DWT using db4 wavelet basis. Decomposition level L = 3. . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Wavelet coefficients obtained by oriented 2-D DT-CWT. Farras filters are used for the first stage and Kingsbury Q-shift after. Decomposition level L = 3. . . . . . . . . . . . . . . . . . 4.6 (a): Wavelet quadtree structure of transformed 2-D image. Level of decomposition L = 3. (b): Example of wavelet quadtree structure. . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.7 (a): Non-overlapping grouping arrangement. (b): Overlapping grouping arrangement. . . . . 4.8 Grouping arrangement of wavelet coefficients. . . . . . . . . . . . . . . . . . . . . . . . . .

  The wavelet basis used is db4. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.16 Reconstructed ϵ(r) (top) and σ(r) (bottom) with ℓ 1 -penalized CSI with different N s × N r . The DT-CWT is used. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.17 Reconstructed ϵ(r) of 4 samples in circular dataset for N s × N r = 16 × 24. Gaussian noise of 20 dB is added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Basics of convolutional neural networks 5.1 A typical CNN architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 An example of convolution operation with a kernel size of 3 × 3. The stride s is 1. No zeropadding was added. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 An example of max-pooling operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 An example of early stopping during network training. . . . . . . . . . . . . . . . . . . . . . vi 5.5 Unrolled model derived from the proximal gradient descent algorithm, where the proximal map is replaced with a deep CNN. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6 Unrolled convolutional neural networks for ISPs 6.1 CSI-based unrolled neural network structure. p ∈ {1 . . . N p }, N p number of sub-blocks. k ∈ {1 . . . N k }, N k number of CSI iterations performed inside each sub-block. χ p-0.5 represents a preliminary result of χ at p-th step. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.2 Network structure of Module CNN and its residual sub-block. . . . . . . . . . . . . . . . . . . 6.3 Relative permittivity of reconstructed profiles of Samples #1-#3 from the test dataset. The range of relative permittivity is

  ε and SSIM obtained by applying the trained model on 200 test samples. Gaussian noise of SNR = 20 dB is added to collected fields. . . . . . . . . . . . . . . . . . . . . . . . . 6.6 Inversion error curves of CSI for samples #1-#6. . . . . . . . . . . . . . . . . . . . . . . . . . 6.7 Noise stability: mean and variance of Err χ obtained by CC-ResNet, CSI-Net, and U-net on 200 test samples belonging to the circular database as a function of noise levels. . . . . . . . . . 6.8 Mean and variance of Err χ obtained by CC-ResNet, CSI-Net and U-net on 200 test samples as a function of the number of training pairs. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.9 Mean and variance of Err χ obtained by CC-ResNet, CSI-Net, and U-net on 200 test samples using models trained with five randomly chosen datasets of 2000 data pairs. . . . . . . . . . 6.10 Reconstructed relative permittivity of samples #7-#9 obtained by three models and CSI. The models are trained with the circular dataset. The test dataset lies outside the range of training data. 20 dB Gaussian noise affects the collected scattered fields. . . . . . . . . . . . . . . . . 6.11 Reconstructed relative permittivity of samples #10-#12 obtained by three models and CSI. The models are trained with the circular dataset. The test dataset lies outside the range of training data. 20 dB Gaussian noise affects the collected scattered fields. . . . . . . . . . . . 6.12 Reconstructed relative permittivity obtained by CSI-Net, U-net (CSI and CC-ResNet have been omitted since failing to yield proper reconstruction). The models are trained with the MNIST dataset. The test dataset is outside the range of training data. 20 dB Gaussian noise affects the collected scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.13 Reconstructed relative permittivity obtained by models trained with a complex-valued circular dataset. ℜ{ε} (top) of ground truth is 1.8 and ℑ{ε} (bottom) is 0.5. 20 dB Gaussian noise affects the collected scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.14 Reconstructed relative permittivity obtained by CSI-Net-MF and CSI-Net. 20 dB Gaussian noise affects the collected scattered fields. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . U-net structure used in our work. The dimension of the feature maps passing through each layer is noted as [Number of channels@H × W ]. . . . . . . . . . . . . . . . . . . . . . inversion and wavelet technique 3.1 Normalized contrast error with different noise level and different N s × N r . . . . . . . . . . 3.2 Mean and variance of Err(χ) on 250 samples. . . . . . . . . . . . . . . . . . . . . . . . . . . 4 Group sparsity penalized Contrast Source Method 4.1 Choice of q for different methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Relative error with different N s × N r -SNR = 20 dB . . . . . . . . . . . . . . . . . . . . . . . 4.3 Relative error with different N s × N r -SNR = 15 dB . . . . . . . . . . . . . . . . . . . . . . . 4.4 Relative error with different N s × N r -SNR = 10 dB . . . . . . . . . . . . . . . . . . . . . . . 4.5 Relative error with Kingsbury Q-shift 14-tap filter -SNR = 20 dB . . . . . . . . . . . . . . . . 4.6 Relative error with Kingsbury Q-shift 18-tap filter -SNR = 20 dB . . . . . . . . . . . . . . . . 4.7 Relative error obtained with some common wavelets using CSI-GS-1 . . . . . . . . . . . . . . 4.8 Relative error obtained with some common wavelets using CSI-GS-2 . . . . . . . . . . . . . . 4.9 Relative error obtained with some common wavelets using CSI-GS-3 . . . . . . . . . . . . . . 4.10 Relative error with ℓ 1 -penalized CSI . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.11 Average CPU time (in s) for one iteration on an Intel Core i5-8365U CPU (1.60 GHz) with 16 GByte memory. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.12 Mean and variance of Err(χ) on 250 samples. . . . . . . . . . . . . . . . . . . . . . . . . . . 5 Basics of convolutional neural networks 6 Unrolled convolutional neural networks for ISPs 6.1 Reconstruction quality of three models and CSI. SNR = 20 dB. . . . . . . . . . . . . . . . . . 6.2 Err χ of Test samples #7-#12. SNR = 20 dB . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Err χ of Test samples #13-#18, and 200 test samples. SNR = 20 dB . . . . . . . . . . . . . . . 6.4 Mean and variance of Err χ over 200 test samples. SNR = 20 dB. . . . . . . . . . . . . . . . . 6.5 Err ε of Test samples #1-#6, and 200 test samples. SNR = 20 dB . . . . . . . . . . . . . . . . . 7 Conclusion and perspectives A Orthogonal Daubechies coefficients A.1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . B U-net structure viii C List of publications 83 ix

D

  'une part, de nouvelles approches basées sur une méthode populaire d'inversion, dite de la source de contraste (CSI ou Contrast Source Inversion), sont proposées, qui combinent alors la technique des ondelettes et la régularisation par la parcimonie. D'autre part, une approche assistée par apprentissage profond est établie basée sur la méthode CSI et les réseaux de neurones convolutifs (Convolutional Neural Networks x ou CNN).

Figure 1 :

 1 Figure 1 : La configuration du modèle.

2 ∥χ∥ 2 2.Figure 2 :Figure 3 :

 2223 Figure 2 : Erreur de contraste normalisée en fonction de N s et N r . (a) CSI et (b) sont pour un CSI basé sur des ondelettes db14. Un bruit gaussien de 20 dB est ajouté aux champs diffractés.

Figure 4 :

 4 Figure 4 : (a) : Vérité terrain. Permittivité reconstruite par (b) CSI-Net, (c) CSI et (d) U-net. Un bruit gaussien de 20 dB est ajouté aux champs diffractés.

  error ρ i,k and the object error r i,k on χ k 11: until a stopping criterion is satisfied.

  signals can be classified into continuous signals and discrete signals. The CWT can tackle continuous signals such as the recording of a speech signal or a music signal, which can be formulated as functions of time.

  A two-dimensional(2-D) image can be represented as a matrix of pixels. When applying the two-dimensional DWT to process an image, the pixels are decomposed into wavelet coefficients, which depend on the change rate of the pixel values. The areas where the pixel values rapidly change correspond to the high-frequency wavelet coefficients, whereas the areas that gradually change are transformed into low-frequency components.

Figure 3 . 1 :

 31 Figure 3.1: Decompose a m × n image X by 2-D DWT of level 1.

Figure 3 . 2 :

 32 Figure 3.2: Decompose a m × n image X by 2-D DWT of level 3.

  Figure 3.3: Wavelet coefficients generated by 2-D DWT. Left: 1-level decomposition. Right: 3-level decomposition.

  First, define the wavelet coefficients of the contrast source J and χ as γ = W{J } (3.24) β = W{χ} (3.25) the data equation error ρ i,k and the state equation error r i,k are also transformed into wavelet domain

1 Figure 3 . 4 :

 134 Figure 3.4: Localization of wavelet coefficients employed in the optimization steps.

  Figure 3.4 gives an example of the optimization steps. Suppose that the optimization starts from L max = 3, at this step, only the approximation coefficients LL 3 are updated, as displayed in 3.4(a). This first step provides the inversion with a stable starting point. In the next steps, more approximation coefficients at different scales are involved in the inversion until the decomposition reaches the coarsest scale, as shown in 3.4(c).

Algorithm 2

 2 Fig. 3.5(b), the region of interest is contained in free space C. N s sources and N r receivers are evenly distributed along an observation circle of radius r = 3.75 m. The operating frequency is 500 MHz, and the collected data are corrupted by Gaussian noise. Wavelet-based CSI (WCSI) 1: Initialize J and χ by back-propagation 2: Initialize the level of decomposition L = L max 3: while L > 0 do error ρ i,k and the object error r i,k .

2 ∥x∥ 2 2

 22 13:Reconstruct the spatial-domain contrast by χ k = W * β k 14:Enforce positivity constraint on χ k 15:until a stopping criterion is satisfied 16:L = L -1 17:Set current result as the initialization of the next step 18: end while To quantitatively evaluate the quality of reconstruction, we define the normalized reconstructed contrast error asErr(x) = ∥x -x∥ 2 {x = χ, ℜ(χ), ℑ(χ)} (3.30)where x is the parameter to evaluate and x is the reconstruction of x, ℜ(•), ℑ(•) represents real and imaginary part, respectively.

Figure 3

 3 Figure 3.5: Left: Austria model, relative permittivity ε(r) = 2. Right: scattering system.

Figure 3 . 6 :

 36 Figure 3.6: Normalized contrast error as a function of N s and N r (arbitrary threshold applied at 50%). (a) CSI,(b), (c) and (d) wavelet-based CSI with Haar, db4 and db14 wavelets respectively. Gaussian noise of 20 dB is added to the scattered fields.

Figure 3 . 7 :Figure 3 . 8 :

 3738 Figure 3.7: Normalized contrast error as a function of N s and N r (arbitrary threshold applied at 50%). (a) CSI,(b), (c) and (d) wavelet-based CSI with Haar, db4 and db14 wavelets respectively. Gaussian noise of 10 dB is added to the scattered fields.

Figure 3 . 9 :

 39 Figure 3.9: Reconstructed ϵ(r) of "Austria profile" for various combinations of receiver and source numbers. From top to bottom: results of CSI, WCSI with Haar, db4, db14 basis, respectively. Gaussian noise of 20 dB is added.

Figure 3 . 10 :

 310 Figure 3.10: Reconstructed ϵ(r) of "Austria profile" for various combinations of receiver and source numbers. From top to bottom: results of CSI, WCSI with Haar, db4, db14 basis, respectively. Gaussian noise of 10 dB is added.

Figure 3 . 11 :

 311 Figure 3.11: Reconstructed ϵ(r) of 8 samples in circular dataset for N s × N r = 6 × 18. Gaussian noise of dB is

Fig. 4

 4 Fig. 4.1(a) and 4.1(b) show the analysis and synthesis filter banks of the DT-CWT. h 0 (n), h 1 (n) denote the low pass/high pass filter pair for the upper filter bank (real branch). g 0 (n), g 1 (n) denote the low pass/high pass filter pair for the lower filter bank (imaginary branch). The wavelet functions Ψ h (r) and Ψ g (r) associated with the two DWTs satisfy Equation (4.2), and Ψ g (r) is designed to be approximately the Hilbert transform of Ψ h (r) so that the corresponding complex wavelet Ψ c (r) is approximately analytic. The in-

  Figure 4.1: 1-D DT-CWT. h (1) 0 (n) is the filter used for the first stage, and h 0 (n) is for the succeeding stages. The implementation of the DT-CWT requires that the first stage of the filter banks be different from the succeeding stages.

  Fig. 4.3 illustrates the twelve wavelets, the six wavelets displayed on the first (second) row are interpreted as the real (imaginary) part of a set of six complex wavelets.2-D DT-CWT has good shift invariance and directional selectivity, which are two key advantages over DWT. The good directional selectivity allows capturing more oriented singularities of an image, without the presence of the checkerboard artifact in the diagonal wavelet of DWT. The near-shift invariance also brings a key property for our work, which will be further introduced in the next part. Moreover, although the 2-D DT-CWT produces wavelet coefficients that are four times more expansive than that of DWT, its implementation is still very efficient since only additions and subtractions for respective sub-bands of all DWTs are required.

Fig. 4

 4 Fig. 4.4(a) is a spatial profile to be decomposed. Fig. 4.4(b) shows its wavelet coefficients obtained by conventional 2-D DWT with db4 wavelet basis. In Fig. 4.5, wavelet coefficients obtained by 2-D DT-CWT are displayed. Farras filters (Abdelnour et al. 2004) are used for the first stage and Kingsbury Q-shift 6-tap filters (Kingsbury 2000) for the following stages. As previously mentioned, the result of DT-CWT is four times more expansive than the result of DWT.

Figure 4 . 2 :Figure 4 . 3 :

 4243 Figure 4.2: db4 wavelets (LH, HL, HH) associated with the 2-D DWT in the space domain. Level of decomposition is 5.

Figure 4 . 4 :

 44 Figure 4.4: (a): A "circle" example in spatial domain. (b): Wavelet coefficients obtained by 2-D DWT using db4 wavelet basis. Decomposition level L = 3.

Figure 4 . 5 :

 45 Figure 4.5: Wavelet coefficients obtained by oriented 2-D DT-CWT. Farras filters are used for the first stage and Kingsbury Q-shift after. Decomposition level L = 3.

Fig. 4

 4 Fig. 4.6(a) shows the 2-D DWT result with the dependency between wavelet coefficients visualized. This dependency can be modeled as a set of trees. Denote three sub-bands of detail coefficients at level L as {HL, LH, HH} L . Each coefficient in {HL, LH, HH} L possesses four children coefficients at the same orientation in {HL, LH, HH} L-1 . The wavelet coefficient from the coarser scale of decomposition can

Figure 4 . 6 :

 46 Figure 4.6: (a): Wavelet quadtree structure of transformed 2-D image. Level of decomposition L = 3. (b): Example of wavelet quadtree structure.

Fig. 4

 4 Fig. 4.7(a) displays an example of non-overlapping groups, where each row of the 2-D matrix represents a subgroup.

Figure 4 . 7 :

 47 Figure 4.7: (a): Non-overlapping grouping arrangement. (b): Overlapping grouping arrangement.

Figure 4 . 8 :

 48 Figure 4.8: Grouping arrangement of wavelet coefficients.

Algorithm 3

 3 Group sparsity penalized CSI 1: Initialize J and χ by back-propagation 2: repeat 3:Calculate the data error ρ i,k and the object error r i,k .

  Figure 4.9: Parametric study for the choice of q ℜ and q ℑ . The final error Err (χ) is plotted as a function of q ℜ and q ℑ and compared to the one of CSI (green line). N s × N r = 32 × 32

Figure 4 .

 4 Figure 4.10: "Austria" model -Relative permittivity ϵ(r) and conductivity σ(r) (mS/m)

  Fig. 4.14.

Figure 4 . 11 :Figure 4 . 12 :

 411412 Figure 4.11: Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 20 dB is added to the scattered fields.

Figure 4 . 13 :

 413 Figure 4.13: Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 15 dB is added to the scattered fields.

Figure 4 . 14 :

 414 Figure 4.14: Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 10 dB is added to the scattered fields.

Figure 4 . 15 :

 415 Figure 4.15: Reconstructed ϵ(r) (row 1, 3, 5) and σ(r) (row 2, 4, 6) of "Austria profile" for various combinations of receiver and source numbers: N r × N s = 32 × 32 (row 1, 2), N r × N s = 18 × 26 (row 3, 4), N r × N s = 16 × 24 (row 5, 6). Gaussian noise of 20 dB is added to the scattered fields. The wavelet basis used is db4.

Figure 4 . 16 :

 416 Figure 4.16: Reconstructed ϵ(r) (top) and σ(r) (bottom) with ℓ 1 -penalized CSI with different N s × N r . The DT-CWT is used.

Figure 4 . 17 :

 417 Figure 4.17: Reconstructed ϵ(r) of 4 samples in circular dataset for N s × N r = 16 × 24. Gaussian noise of 20 dB is added.

5. 2 . 1

 21 Figure 5.1: A typical CNN architecture

Figure 5 . 2 :

 52 Figure 5.2: An example of convolution operation with a kernel size of 3 × 3. The stride s is 1. No zero-padding was added.

Figure 5 . 3 :

 53 Figure 5.3: An example of max-pooling operation

  verges on the training dataset to prevent overfitting. To perform early stopping, we keep a small portion (usually 20% -30%) of the training dataset as a validation set which is then used for cross-validation purposes. After each training epoch, the trained model will be applied to this validation set to test its performance. When the performance of the model starts to deteriorate in the next epoch, the training will stop because the generalization ability of the model begins to decline at this point. An illustration of the early stopping method is shown in Fig.5.4.

Figure 5 . 4 :

 54 Figure 5.4: An example of early stopping during network training.

Figure 5 . 5 :

 55 Figure 5.5: Unrolled model derived from the proximal gradient descent algorithm, where the proximal map is replaced with a deep CNN.

Figure 6 . 1 :

 61 Figure 6.1: CSI-based unrolled neural network structure. p ∈ {1 . . . N p }, N p number of sub-blocks. k ∈ {1 . . . N k }, N k number of CSI iterations performed inside each sub-block. χ p-0.5 represents a preliminary result of χ at p-th step.

Figure 6 . 2 :

 62 Figure 6.2: Network structure of Module CNN and its residual sub-block.

  Since the proposed model has a cascade architecture, the choice of the loss function is important. The loss scale corresponding to each sub-block's output could significantly differ. To equally train each sub-block, a combined loss function built with the normalized mean square error (MSE) is chosen:

1 .

 1 Relative error of permittivity over the pixelized domain Err ε = 1 N x × N y   Nx ix Ny iy εix,iy -ε ix,iy ε ix,iy

Figure 6 . 3 :

 63 Figure 6.3: Relative permittivity of reconstructed profiles of Samples #1-#3 from the test dataset. The range of relativepermittivity is[1.3, 2]. Gaussian noise of SNR = 20 dB is added.

Figure 6 . 4 :

 64 Figure 6.4: Relative permittivity of reconstructed profiles of samples #4-#6 from the test dataset. The range of relativepermittivity is[1.3, 2]. Gaussian noise of SNR = 20 dB is added.

Figure 6 . 5 :Figure 6 . 6 :Figure 6 . 7 :Figure 6 . 8 :

 65666768 Figure 6.5: Err χ , Err ε and SSIM obtained by applying the trained model on 200 test samples. Gaussian noise of SNR = 20 dB is added to collected fields.

Figure 6 . 9 :

 69 Figure 6.9: Mean and variance of Err χ obtained by CC-ResNet, CSI-Net, and U-net on 200 test samples using models trained with five randomly chosen datasets of 2000 data pairs.

  The scattering system settings, dataset division, and training hyperparameters are un-

Figure 6 . 10 :

 610 Figure 6.10: Reconstructed relative permittivity of samples #7-#9 obtained by three models and CSI. The models are trained with the circular dataset. The test dataset lies outside the range of training data. 20 dB Gaussian noise affects the collected scattered fields.

  1.

						𝑿 𝑳𝑳,𝟑
					𝑿 𝑳𝑳,𝟐	Level 1 DWT	𝑿 𝑳𝑯,𝟑 𝑿 𝑯𝑳,𝟑	Level 3 coefficients
			𝑿 𝑳𝑳,𝟏 𝑿 𝑳𝑯,𝟏	Level 1 DWT	𝑿 𝑳𝑯,𝟐 𝑿 𝑯𝑳,𝟐	coefficients Level 2	𝑿 𝑯𝑯,𝟑
	𝑿[𝒎, 𝒏]	Level 1 DWT	𝑿 𝑯𝑳.𝟏	coefficients Level 1	𝑿 𝑯𝑯,𝟐
			𝑿 𝑯𝑯,𝟏

  1 or higher levels, the approximation coefficients are always located within the upper left corner, and the rest are detail coefficients. Note that the side length of the input image must be a multiple of 2 L where L is the decomposition level. The maximum level that can be reached is L max = log 2

	√	N where	√	N is the
	side length of the input image.

  The data and object error decrease smoothly for CSI. Since WCSI has multiple inversion stages, its error curve decreases again after a period of convergence each time it completes an inversion step. The CPU times for CSI and WCSI are 47.68 s and 168.18 s, respectively, on an Intel Core i7-6700HQ CPU (2.60 GHz) with 32 GByte memory. More iterations are required for WCSI, resulting in longer computation time. Further improvement on the stopping criteria of WCSI can reduce the computation burden.

.

1. Fig. 3.8

shows the inversion curves of CSI and WCSI with db14 wavelets where N s × N r = 6 × 18.

Table 3 . 1 :

 31 Normalized contrast error with different noise level and different N s × N r Noise (in dB) N s × N r CSI WCSI-Haar WCSI-db4 WCSI-db14

		48 × 48 0.0815 0.1375	0.0891	0.0766
	20 dB	18 × 12	NA	0.2102	0.1611	0.1404
		12 × 9 0.8040 0.2077	0.1394	0.1147
		6 × 18 0.3951 0.2107	0.1553	0.1362
		48 × 48 0.1063 0.1486	0.1011	0.0874
	10 dB	18 × 12	NA	0.1548	0.1066	0.0908
		12 × 9 1.0863 0.4037	0.2553	0.2373
		6 × 18 0.9818 0.4366	0.2939	0.2772

Table 3 . 2 :

 32 Mean and variance of Err(χ) on 250 samples.

	CSI WCSI-db14
	Mean 0.5681	0.1476
	Variance 1.9951	0.1118

  Gaussian noise of dB is added. 4 Group sparsity penalized Contrast Source Method One central feature of DWT is that it can provide a sparse representation for a given signal. As mentioned in the previous chapter, the wavelet coefficients can be divided into two main parts: approximation coefficients and detail coefficients. Approximation coefficients which represent the low-frequency components concentrate most of the energy of the signal. Detail coefficients that contain high-frequency components are sparse. In this chapter, a group sparsity penalized CSI method is proposed based on this important property of wavelet detail coefficients.

Table 4 . 1 :

 41 Choice of q for different methods

		CSI-GS-1 CSI-GS-2 CSI-GS-3
	q ℜ	70	80	60

Table 4 . 2 :

 42 Relative error with different N s × N r -SNR = 20 dB N s × N r the same group, and adjacent pixels do not share the same children.The relative errors of the reconstructed contrast are in Table4.2. It can be observed from the results that the proposed methods provide a better quality of reconstruction. When the number of sensors is reduced, the relative error obtained by CSI clearly increases, whereas the reconstruction quality remains relatively stable with them.

		Error	CSI CSI-GS-1 CSI-GS-2 CSI-GS-3
		Err(χ)	0.34	0.14	0.14	0.16
	32 × 32	Err(ℜ{χ}) 0.27	0.12	0.12	0.15
		Err(ℑ{χ}) 0.61	0.21	0.21	0.23
		Err(χ)	0.42	0.15	0.16	0.18
	18 × 26	Err(ℜ{χ}) 0.35	0.13	0.14	0.16
		Err(ℑ{χ}) 0.70	0.23	0.22	0.25
		Err(χ)	0.51	0.17	0.19	0.20
	16 × 24	Err(ℜ{χ}) 0.47	0.16	0.18	0.19
		Err(ℑ{χ}) 0.67	0.22	0.22	0.24
	four children in

Table 4 . 3 :

 43 Relative error with different N s × N r -SNR = 15 dB N s × N r

		Error	CSI CSI-GS-1 CSI-GS-2 CSI-GS-3
		Err(χ)	0.55	0.16	0.17	0.18
	32 × 32	Err(ℜ{χ}) 0.47	0.14	0.15	0.16
		Err(ℑ{χ}) 0.86	0.23	0.23	0.27
		Err(χ)	0.56	0.18	0.19	0.20
	18 × 26	Err(ℜ{χ}) 0.47	0.16	0.18	0.18
		Err(ℑ{χ}) 0.88	0.26	0.26	0.28
		Err(χ)	0.58	0.22	0.25	0.24
	16 × 24	Err(ℜ{χ}) 0.49	0.20	0.24	0.22
		Err(ℑ{χ}) 0.92	0.27	0.28	0.28

Table 4 . 4 :

 44 Relative error with different N s × N r -SNR = 10 dB N s × N r

		Error	CSI CSI-GS-1 CSI-GS-2 CSI-GS-3
		Err(χ)	0.89	0.21	0.27	0.23
	32 × 32	Err(ℜ{χ}) 0.76	0.18	0.24	0.20
		Err(ℑ{χ}) 1.40	0.33	0.35	0.34
		Err(χ)	0.94	0.28	0.38	0.29
	18 × 26	Err(ℜ{χ}) 0.79	0.26	0.39	0.26
		Err(ℑ{χ}) 1.53	0.36	0.36	0.43
		Err(χ)	0.72	0.26	0.32	0.28
	16 × 24	Err(ℜ{χ}) 0.57	0.24	0.31	0.24
		Err(ℑ{χ}) 1.32	0.35	0.36	0.40

Table 4 .5:

 4 Relative error with Kingsbury Q-shift 14-tap filter -SNR = 20 dB N s × N r

		Error	CSI-GS-1 CSI-GS-2 CSI-GS-3
		Err(χ)	0.14	0.14	0.16
	32 × 32	Err(ℜ{χ})	0.13	0.12	0.14
		Err(ℑ{χ})	0.21	0.21	0.23
		Err(χ)	0.15	0.16	0.17
	18 × 26	Err(ℜ{χ})	0.13	0.14	0.15
		Err(ℑ{χ})	0.23	0.22	0.25
		Err(χ)	0.17	0.19	0.19
	16 × 24	Err(ℜ{χ})	0.16	0.18	0.18
		Err(ℑ{χ})	0.22	0.22	0.23

Table 4 . 9 :

 49 Relative error obtained with some common wavelets using CSI-GS-3 N s × N r

		Error	Haar D4 D8 C2 C4
		Err(χ)	NA 0.25 0.23 0.19 0.23
	32 × 32	Err(ℜ{χ}) NA 0.21 0.18 0.16 0.19
		Err(ℑ{χ}) NA 0.42 0.42 0.32 0.40

Table 4 .

 4 10: Relative error with ℓ 1 -penalized CSI 32 × 32 18 × 26 16 × 24

	Err(χ)	0.22	0.27	0.35
	Err(ℜ{χ}) 0.18	0.24	0.34
	Err(ℑ{χ}) 0.37	0.43	0.38

Table 4 .

 4 11: Average CPU time (in s) for one iteration on an Intel Core i5-8365U CPU (1.60 GHz) with 16 GByte memory.

Table 4 .

 4 12: Mean and variance of Err(χ) on 250 samples.

	CSI CSI-GS-1
	Mean 0.6996 0.4978
	Variance 1.9205 0.9530

  1 , β 2 are exponential decay rates; 3: Update t = t + 1; m t = β 1 m t-1 + (1 -β 1 )g t , g t is the gradient of loss function with respect to W;4: Update v t = β 2 v t-1 + (1 -β 2 )g 2 Compute bias-corrected moment estimates, m t = m t /(1 -β 1,t ), v t = v t /(1 -β 2,t ); 6: Update the parameters W t = W t-1 -αm t /( √ v t + δ), α is exponential decaying step size; 7: ∥W t -W t-1 ∥ 2 < ξ, stop, otherwise, repeat steps from 3 to 6.

	t ;
	5:

  iteration. prox tH( ) (x) is the proximal operator of H( ) at point x and it is a function of H and t

	prox tH( ) (t) = arg min x	1 2t	∥x -t∥ 2 2 + H(x)

Table 6 . 5 :

 65 Err ε of Test samples #1-#6, and 200 test samples. SNR = 20 dB

	Model	Sample #1 #2	#3	#4	#5	#6 Mean
	CSI-Net-MF	0.102	0.101 0.118 0.118 0.071 0.062 0.0551
	CSI-Net	0.201	0.208 0.220 0.613 0.236 0.232 0.134

  × 10 -3 7.832 52 × 10 -4 6.756 36 × 10 -3 -1.523 81 × 10 -3 (b) Coefficients of the scaling functions for D12-20 .547 904 72 × 10 -3 -6.887 719 26 × 10 -3 -6.679 620 23 × 10 -3

	1	0.683 012 7	0.470 467 21 0.325 803 43 0.226 418 98	0.157 43
	1	1.183 012 7	1.141 116 92 1.010 945 72 0.853 943 54	0.699 81
		0.316 987 3	0.650 365	0.892 200 14 1.024 326 94	1.062 76
		-0.183 012 7 -0.190 934 42 -0.039 575 03 0.195 766 96	0.445 32
		-0.120 832 21 -0.264 507 17 -0.342 656 71	-0.319 60
			0.049 817 5	0.043 616 3 -0.045 601 13	-0.183 06
				0.046 503 6	0.109 702 65	0.137 09
			-0.014 986 99 -0.008 826 80	0.038 21
					-0.017 791 87	-0.044 75
		4.717 427 93 D14 D16 D18	D20
	0.110 099 43	0.076 955 62	0.053 850 35	0.037 16
	0.560 791 28	0.442 467 25	0.344 834 30	0.266 18
	1.031 148 49	0.955 486 15	0.855 349 06	0.745 07
	0.664 372 48	0.827 816 53	0.929 545 71	0.973 11
	-0.203 513 82	-0.022 385 74	0.188 369 55	0.397 74
	-0.316 835 01	-0.401 658 63	-0.414 751 76	-0.353 20
	0.100 846 7	6.681 940 92 × 10 -4 -0.136 953 55	-0.277 88
	0.114 003 45	0.182 076 36	0.210 068 34	0.180 45
	-0.053 782 45	-0.024 563 90	0.043 452 675	0.131 99
	-0.023 439 94	-0.062 350 21	-0.095 647 26	-0.100 57
	0.017 749 79	0.019 772 16	3.548 928 13 × 10 -4 -0.041 25
	6.075 149 95 × 10 -4	0.012 368 84	0.031 624 17	0.046 81
	-2

× 32

Err(χ) NA 0.18 0.19 0.18 0.21 Err(ℜ{χ}) NA 0.15 0.16 0.15 0.17 Err(ℑ{χ}) NA 0.30 0.34 0.30 0.37

× 26

Err(χ) 0.25 0.22 0.24 0.21 0.27 Err(ℜ{χ}) 0.30 0.18 0.20 0.18 0.22 Err(ℑ{χ}) 0.30 0.36 0.39 0.33 0.43

× 24

Err(χ) 0.33 0.25 0.29 0.28 0.32 Err(ℜ{χ}) 0.29 0.23 0.27 0.27 0.30 Err(ℑ{χ}) 0.49 0.32 0.37 0.35 0.39 [START_REF] Wang | Image quality assessment: from error visibility to structural similarity[END_REF])

where µ x and σ x are the mean and the standard deviation of x, respectively. σ xy is the covariance between x and y. C 1 = (k 1 L) 2 , C 2 = (k 2 L) 2 are to stabilize the division with weak denominator where L the dynamic range of the pixel values. k 1 = 0.01 and k 2 = 0.03 by default.

To check on the advantage brought by Module CSI, one removes all Module CSI blocks from CSI-Net, and the model becomes a cascaded residual neural network, denoted as CC-ResNet. For comparison, results of CC-ResNet, the widely used U-net structure (detailed in Appendix B), and the iterative CSI method with positivity constraint are also displayed. For the latter, the number of iterations has been arbitrarily chosen as N CSI k = 150 and the algorithm is stopped when N CSI k is reached. Computations are run on a server with NVIDIA RTX A3000 GPU (6 GB). The networks are implemented with Pytorch [START_REF] Paszke | Pytorch: An Imperative Style, High-performance Deep Learning Library[END_REF]. It takes 0.63 h for CSI-Net training and 0.03 h for CC-ResNet and U-net. It takes less than 0.4 s to get the reconstruction using CSI-Net, and about 4 s using CSI, U-Net being almost real-time.

Implementation details

CSI-Net results at single frequency

In this section, numerical experiments are carried out to evaluate the performance of the proposed methods. The details of the scattering system are as follows. The operating frequency is 400 MHz, the region of interest D is of 2 m×2 m size and is discretized into N square cells. For the forward problem, N = 32 × 32. N r = 32 receivers and N s = 16 transmitters are evenly distributed on a circle of radius r = 3 m centered at (0, 0) m. The embedding medium is of ϵ b = 1.

Test with circular dataset

A synthetically generated circular dataset (Wei et al. 2019a) is used. Each sample contains disks of random quantities: number (between 1 and 3), relative permittivity (between 1.3 and 2), radius (between 0.15 m and 0.4 m), and location. The scatterers may overlap. The collected fields used for training are noise-free. Unless specified, the scatterers are lossless.

The trained models are tested with 200 samples. To check on robustness, further simulations are carried out with Gaussian noises of different SNR (20, 10, and 5 dB) added to the scattered fields. In Fig. 6.7, the mean and variance of results obtained with the real-valued circular dataset are displayed. For the mean error, the growth rate of CSI-Net at 20 and 10 dB is larger than for the two other approaches, yet CSI-Net always exhibits the smallest mean error, i.e., it is robust to noise contamination. As for the variance change, improvement by the proposed method lies not only in better overall reconstruction quality but also in a strong reduction of poor results, i.e., confidence level increases.

In general, a larger training dataset leads to better neural network training. However, in the real world, collecting datasets may not be straightforward. So, a neural network should perform well with a small amount of data. Here, the models are trained with the real-valued circular dataset, and the numbers of training pairs are 250, 500, 1000, 1500, and 2000, respectively. The mean value and variance of Err χ are displayed in Fig. 6.8. As the training dataset's size increases, the performance of all models improves.

But, for a small dataset, U-net performance strongly degrades while CC-ResNet and CSI-Net remain stable.

Moreover, CSI-Net consistently outperforms the other two methods.

Evaluations of reliability can be performed by showing that the training dataset is not particularly tuned to The trained models are tested on 200 samples within the same range of the training dataset. As shown in Fig. 6.12, samples #13-#16 are four representative samples among the test dataset, and #17-#18 are "Austria" profiles with relative permittivities of 2.5 and 3, respectively. Due to strong nonlinearity, CSI fails to reconstruct all samples, so its results are not displayed. CC-ResNet is not included since all of its results are worse than CSI-Net. U-net can roughly reconstruct the profiles, however, most details are missing. In comparison, CSI-Net achieves better reconstructions for all test samples.

The Err χ of these test samples obtained by different methods are presented in Table 6.3. Note that its last two columns give the mean value and the variance of Err χ across all 200 test samples. We can conclude from above results that CSI-Net works better than other methods for the cases where the nonlinearity increases.

Test with lossy scatterers

All previous results are on lossless scatterers. The proposed model is now trained with a complex-valued circular dataset to represent lossy scatterers. The training dataset shares the same key parameters as the real-valued one except for the relative permittivity value. Real parts of relative permittivities are between 1.3 and 2.4; imaginary parts vary from 0 to 1.2. The trained model is tested on the "Austria" profile, and the results are in Fig. 6.13. Both CSI-Net and CSI exhibit satisfactory results. On the contrary, U-net fails to reconstruct the ring for both real and imaginary parts. Corresponding Err χ are in Table 6.4. The last two rows give the mean error and variance over 200 test samples. ℜ{ε} (top) of ground truth is 1.8 and ℑ{ε} (bottom) is 0.5. 20 dB Gaussian noise affects the collected scattered fields.

Extension of CSI-Net to multi-frequency data

The frequency hopping approach (Chew and Lin 1995) [START_REF] Bucci | Inverse Scattering Problems with Multifrequency Data: Reconstruction Capabilities and Solution Strategies[END_REF] is a popular technique to improve reconstruction performance reached with the monochromatic data.

Using multi-frequency data has the following advantages:

• First of all, even if the problem is still ill-posed, another type of relation is treated, which improves the accuracy of the reconstruction.

• Second, the use of multi-frequency data reduces ill-conditioning (i.e., improves local stability) of the reconstruction process.

• Third, robustness against false solutions is improved.

In this section, we apply the frequency-hopping technique to improve CSI-Net. The recurrent network structure of CSI-Net makes it naturally suitable for applying the frequency hopping technique. For different sub-blocks, we input data of different frequencies. First, a reliable starting guess can be obtained with low-frequency data. Next, the results are used for the following inversion steps performed at higher frequencies. Denote CSI-Net in the multi-frequency scheme as CSI-Net-MF. Take a CSI-Net that contains 4 sub-blocks (exactly the same structure as that in previous sections), for each sub-block, we input the data at a fixed frequency. Here, we use data obtained at 4 frequencies: 100 MHz, 200 MHz, 400 MHz, and 600 MHz.

To evaluate the performance of the CSI-Net-MF, we apply it to a more complicated dataset that contains profiles of higher nonlinearity. The error of permittivity Err ε is used here to quantitatively evaluate the reconstruction quality. A synthetically generated lossless circular dataset is used. Each sample contains disks of random quantities: number (1 ∼ 3), relative permittivity (2.5 ∼ 3.

B U-net structure

The U-net structure, originally designed for medical image segmentation, is widely used for ISPs (as in (Wei et al. 2019a;[START_REF] Liu | Physics-guided Loss Functions Improve Deep Learning Performance in Inverse Scattering[END_REF]) among many). For comparison, a U-net structure adapted to our work is proposed, and the results are displayed. Details of U-net can be found in [START_REF] Ronneberger | U-net: Convolutional Networks for Biomedical Image Segmentation[END_REF].