Prof Amélie Litman

Dr Oliver Dorn

lovely wife Laurence Lesselier who welcome to me Paris from the rst day and made me feel like at home. I would like to thank Aurélia Fraysse, Maître de Conférences, for her advices and co-operation as well. I would like

Résumé en Français

Ce manuscrit de thèse décrit mes travaux de recherche dans le domaine de l'imagerie micro-onde, qui fait l'objet d'un intérêt de recherche important en raison de son potentiel en tant que technique pratique et ecace pour des systèmes médicaux [1,[START_REF] Henriksson | Quantitative microwave imaging for breast cancer detection using a planar 2.45 GHz system[END_REF], des caractérisations de matériaux [3], la caractérisation du sous-sol, la télédétection, et les essais et évaluations non destructifs [4,[START_REF] Pastorino | Microwave Imaging[END_REF]. L'objectif de l'imagerie micro-onde est notamment d'estimer l'emplacement de diracteurs dans une région d'intérêt, en accédant en sus à la distribution de leurs propriétés électromagnétiques.

Ces dernières années, le besoin de méthodes et de techniques de reconstruction ecaces a fortement émergé pour résoudre les problèmes de diraction inverse électromagnétique qui se posent dans l'imagerie micro-onde an d'en réduire au mieux les dicultés théoriques et pratiques. La demande dans diverses applications impose le développement de méthodes ecaces et précises. Cependant, la mise en ÷uvre d'algorithmes de reconstruction qui soient stables, ables et ecaces est un dé en raison de la non-linéarité des équations de la diraction et du caractère mal-posé du problème auquel on se confronte [4,[START_REF] Di Benedetto | Numerical linear algebra for nonlinear microwave imaging[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF].

En mathématique, il y a souvent deux problèmes qui font face l'un à l'autre, le problème direct et le problème inverse. Dans le problème direct, l'état observable d'un système est déterminé en utilisant tous les paramètres nécessaires. En comparaison, le problème inverse se préoccupe de la condition préalable qui crée une donnée observée. Mathématiquement, nous avons y = K(x), où x ∈ X est l'inconnue et y ∈ Y est le vecteur des observations. Habituellement, K() est un opérateur continu bien posé et la solution du problème direct (trouver y donné x) ne rencontre aucun obstacle signicatif. D'autre part, le mapping inverse de y à x dans les problèmes d'intérêt n'est pas aisé à traiter.

Normalement, les problèmes inverses sont tels que les modèles mathématiques ne sont pas bien posés dans le sens de Hadamard, signiant pas de solution unique, pas d'existence de solution et la non-dépendance des données de manière continue [START_REF] Di Benedetto | Numerical linear algebra for nonlinear microwave imaging[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF]. Ces problèmes sont dits mal-posés et ils sont la source de nombreuses dicultés numériques car ils rendent la plupart des algorithmes numériques instables pour des perturbations de données [4,[START_REF] Hanke | A convergence analysis of the Landweber iteration for nonlinear ill-posed problems[END_REF].

Les méthodes numériques qui peuvent traiter ces problèmes sont appelées méthodes de régularisation.

Les problèmes de diraction électromagnétique inverse possèdent donc deux propriétés fondamentales qui font de leur résolution une tâche dicile, la première est la non-linéarité des équations de diraction [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF], la seconde le caractère mal-posé du problème [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF].

Le champ diracté est une fonction intégrale du champ total induit dans le domaine investigué, multiplié par le contraste (ici la diérence entre les permittivités complexes du domaine d'investigation et du milieu environnant). Le champ total induit est également une fonction non linéaire de la permittivité [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. En conséquence, la reconstruction de cette permittivité à partir de champs diractés est un problème non linéaire.

En outre, le fait que le champ diracté soit une fonction intégrale du champ total et de la permittivité dans le domaine sous investigation conduit à un problème inverse malposé [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. L'opérateur intégral supprime la contribution des composantes à variation rapide de la permittivité au champ diracté et rend dicile leur identication [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF].

Cette nature physique du phénomène de diraction est combinée avec le fait que seul un ensemble ni d'échantillons, non nécessairement aussi de distribution optimale par rapport aux diracteurs recherchés et insusamment indépendants les uns des autres, peut être collecté par les mesures, mesures qui ne fournissent pas susamment d'information, et que ces échantillons sont presque toujours corrompus par le bruit, de sorte qu'un modèle approximatif discrétisé doit remplacer la description mathématique exacte du problème.

Ainsi, cela rend le problème inverse sévèrement mal posé [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. La force de la non-linéarité dans le problème de diraction s'accroît avec l'amplitude des champs diractés couplés avec le/les objet/s dans le domaine investigué. Par conséquent, une permittivité supérieure (cette notion étant aussi dépendante de la taille des diracteurs) signie une non-linéarité plus forte. La plupart des méthodes déterministes qui se révèlent ecaces pour résoudre des problèmes inverses électromagnétiques impliquant des faibles diracteurs utilisent des approximations linéaires de premier ordre [START_REF] Devaney | Mathematical Foundations of Imaging, Tomography and Waveeld inversion[END_REF], telles que la tomographie par diraction [START_REF] Zorgati | Eddy current testing of anomalies in conductive materials. i. qualitative imaging via diraction tomography techniques[END_REF], l'approximation de Born de premier ordre [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Rajan | A comparison between the Born and Rytov approximations for the inverse backscattering problem[END_REF], l'approximation de Kirchho [START_REF] Pastorino | Microwave Imaging[END_REF] et l'approximation de Rytov [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Rajan | A comparison between the Born and Rytov approximations for the inverse backscattering problem[END_REF]. Bien que ces méthodes soient moins exigeantes sur le plan informatique, elles ne fournissent pas de solutions précises lorsque de puissants diracteurs sont présents dans le domaine étudié. Dans de tels cas, lorsque la non-linéarité est forte, la linéarisation ne produit pas une approximation précise du problème non linéaire et il faut utiliser des techniques plus rigoureuses telles que le solveur non linéaire de Newton combiné avec des formulations de champ de contraste ou de source de contraste [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Bozza | An inexact Newton-based approach to microwave imaging within the contrast source formulation[END_REF], la méthode de gradient con-jugué non linéaire [START_REF] Abubakar | Contrast source inversion method: State of art[END_REF], l'algorithme de descente de plus grande pente non linéaire (par ex. Landweber non-linéaire) [START_REF] Li | Image reconstruction by nonlinear Landweber iteration for complicated distributions[END_REF][START_REF] Hettlich | The Landweber iteration applied to inverse conductive scattering problems[END_REF], et le schéma de Levenberg-Marquardt [START_REF] Franchois | Microwave imaging-complex permittivity reconstruction with a Levenberg-Marquardt method[END_REF], ainsi que les méthodes de Born distordu [START_REF] Chew | Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method[END_REF]. Ces méthodes sont de plus exigeantes en terme de calcul comparativement à celles qui utilisent les schémas de linéarisation du premier ordre. D'autres méthodes qui bénécient du meilleur des deux mondes ont également été développées. Celles-ci utilisent des schémas de linéarisation d'ordre supérieur ou une application itérative de premier ordre. Les exemples incluent l'approximation de Born étendu [START_REF] Zhang | Two nonlinear inverse methods for electromagnetic induction measurements[END_REF] et les approximations de Born d'ordre supérieur [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Estatico | An inexact-Newton method for short-range microwave imaging within the second-order Born approximation[END_REF], la méthode itérative de Born (BIM) [START_REF] Wang | An iterative solution of the two-dimensional electromagnetic inverse scattering problem[END_REF] et la méthode itérative variationnelle [START_REF] Zaiping | Variational Born iteration method and its applications to hybrid inversion[END_REF], respectivement. Les méthodes de ce dernier groupe sont plus ecaces en terme de calcul par rapport aux méthodes d'inversion non linéaires, et s'appliquent également à une plus grande gamme de force des diracteurs que les méthodes de linéarisation de premier ordre.

Indépendamment de la méthode utilisée pour s'attaquer à la non-linéarité, le caractère mal-posé du problème inverse de diraction doit être pris en compte [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF][START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. À cette Donoho [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF], et celui de Rao [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF], entre autre, à propos de l'approximation de fonctions et la sélection optimale de bases.

La régularisation par la parcimonie a été un domaine de recherche fort de ces dernières années. Le caractère bien posé et la convergence des méthodes ont été analysés pour des problèmes inverses linéaires dans [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] et pour des problèmes inverses non linéaires dans [START_REF] Grasmair | Sparse regularization with lp penalty term[END_REF]. Il a été démontré que la régularisation par parcimonie est simple à utiliser et très ecace pour des problèmes inverses avec des solutions parcimonieuses. Cette méthode a été appliquée aux problèmes d'imagerie compressive et de tomographie par impédance électrique (EIT) [3032]. Cette méthode de régularisation conduit à des minimiseurs parcimonieux lorsque p = 1 pour la norme l p et favorise la parcimonie pour 1 < p < 2 [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

Dans ce travail, nous avons exclu p = 0 puisque l 0 est une pseudo-norme (la propriété de distance des normes n'est pas satisfaite), qui compte le nombre d'éléments non nuls.

Même si la parcimonie est le plus favorisé lorsque p = 0, notre problème d'optimisation devient NP-hard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. De plus, ce problème est non-convexe, ce qui crée un problème de minima multiples et la recherche des minima globaux optimaux nécessite une recherche intensive très exigeante en calcul [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. Pour éviter l'apparition de minima multiples, la norme l 1 (quand p ≤ 1) est utilisée comme la meilleure approximation convexe du problème d'optimisation en norme l 0 . Ce concept est connu sous le nom de relaxation convexe [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF] en tant qu'alternative au problème NP-Hard. En outre, la norme l 1 propose en général une solution avec d'excellentes caractéristiques de parcimonie. Par conséquent, nous choisissons une pénalité l 1 an de promouvoir la connaissance a-priori de parcimonie tout au long de ce travail de thèse.

L'utilisation de normes l 1 pour assurer la parcimonie est connue depuis près d'une décennie. L'opérateur LASSO ou Least Absolute Shrinkage and Selection Operator a été introduit dans la littérature statistique [START_REF] Tibshirani | Regression shrinkage and selection via the LASSO[END_REF] et l'algorithme dit de Basis Pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] a été proposé dans la communauté du signal à peu près au même moment. Les avantages les plus importants des systèmes de pénalisation l 1 sont leur convexité, et la forte parcimonie des résultats. Diérentes versions de problèmes de régularisation l 1 peuvent être reformulées en tant que programmation quadratique linéaire ou convexe, contrainte ou non, ou de second ordre (SOC), ce qui autorise des algorithmes ecaces et globalement convergents.

Un autre avantage important de l'utilisation de la pénalisation l 1 est l'existence d'un certain nombre de résultats théoriques récents (par exemple [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF]) montrant que, sous certaines conditions de parcimonie sur le signal inconnu sous-jacent, le signal peut être récupéré exactement.

En fait, la fonction inconnue (par exemple, la fonction de contraste ou le courant équivvi Résumé en Français alent dans les problèmes inverses électromagnétiques) peut être considérée comme parcimonieuse (ou compressible) dans une base appropriée, du fait qu'elle puisse n'occuper qu'une petite partie d'un scénario autrement connu ou qu'elle puisse représenter une anomalie étendue presque homogène contenue dans un milieu de fond. En d'autres termes, un diracteur peut être parcimonieux par rapport à une base et non parcimonieux par rapport à une autre. En outre, dans de nombreux cas pratiques, seul un nombre limité de mesures sont permises, de sorte que des outils de traitement appropriés pour gérer ce type de complexité sont nécessaires.

La question de choisir une base appropriée pour une famille de signaux a fait l'objet de beaucoup d'attention cette dernière décennie, et de nombreuses nouvelles bases ont été introduites, telles les bases d'ondelettes, les ridgelets et curvelets, parmi beaucoup d'autres [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Bien que toute base d'extension minimale pour un espace à dimension nie puisse représenter parfaitement tout signal dans l'espace, lorsque seul un sous-ensemble de signaux possibles est intéressant, certaines bases possèdent de meilleures propriétés de représentation parcimonieuse que d'autres. Certaines applications qui bénécient grandement de la parcimonie de la représentation sont la compression du signal, le débruitage et l'estimation des paramètres [3739]. En compression pour la transmission de l'information, si la représentation du signal n'est pas parcimonieuse, nous devons transmettre le signal complet. Cependant, si, sous un changement de base, la représentation devient parcimonieuse, des économies substantielles sont possibles. La plupart des coecients de la représentation sont très petits (par dénition de la parcimonie) et si nous les mettons à zéro, la qualité perceptuelle du signal sera très peu aectée. Par conséquent, on ne transmet que les grands coecients, qui sont peu nombreux.

Dans cette thèse, nous avons proposé des méthodes qui exploitent la parcimonie en tant que brique principale an de résoudre le problème inverse non linéaire. Initialement, nous avons proposé une méthode où la contrainte de parcimonie est directement appliquée au problème de la reconstruction des propriétés diélectriques internes complexes d'un objet en fonction de la connaissance du champ diracté extérieur qui est généré par l'interaction entre l'objet et un champ incident connu. Le problème d'optimisation non linéaire est résolu par un algorithme itératif de contraction douce an d'appliquer la contrainte de parcimonie. Cette parcimonie est appliquée à chaque itération par une fonction de seuillage doux. En outre, nous utilisons une sélection de pas adaptative selon la règle de Barzilai et Borwein (BB) et une projection sur le contraste inconnu de telle sorte que les contraintes de positivité soient prises en compte.

La base adoptée joue un rôle clé, car celle-ci doit être tant précise qu'ecace. L'ecacité vii nécessite un nombre réduit de coecients dans la représentation, tandis que la précision implique une faible erreur de représentation, qui quantie la désadaptation entre la fonction réelle et sa projection sur la base considérée. Étant donné que le problème inverse de diraction électromagnétique que nous considérons n'est pas linéaire et est mal posé, la réduction du nombre de paramètres inconnus est importante pour réduire le caractère mal posé [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Candes | Curvelets: A surprisingly eective nonadaptive representation for objects with edges[END_REF][START_REF] Li | A contrast source inversion method in the wavelet domain[END_REF]. En fait, la décomposition en ondelettes permet de réduire le nombre de paramètres inconnus par rapport à la représentation en pixels habituellement exploitée an d'accroître la abilité de l'inversion [START_REF] Isernia | On the local minima in a tomographic imaging technique[END_REF]. À cet égard, après notre première méthode proposée avec un seuillage doux, nous avons étendu l'application du paradigme d'imagerie micro-onde aux fonctions non-pixelisées, élargissant ainsi l'ensemble des distributions de diracteurs pouvant être considérés aux cibles de forme et taille arbitraires, sous réserve qu'elles soient parcimonieuses par rapport à une base appropriée. Dans ce but, les ondelettes sont utilisées et nous avons étudié si la base d'ondelettes est un bon choix pour tenir compte ou non du compromis entre l'ecacité et la précision de la représentation dans notre cas.

La transformée en ondelettes permet d'avoir des informations tant sur la variation spatiale que sur le contenu fréquentiel de la fonction, en la représentant en tant que somme pondérée de la version dilatée et translatée de l'ondelette dite mère [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. En particulier, les fonctions de base à grande échelle codent le contenu à basse fréquence de la fonction traitée, tandis que celles à petites échelles prennent compte des détails de haute fréquence ou des détails ns, [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Ceci réduit le nombre de paramètres inconnus par rapport à la représentation habituellement exploitée des pixels, améliorant ainsi la abilité de l'inversion [START_REF] Isernia | On the local minima in a tomographic imaging technique[END_REF][START_REF] Bucci | Electromagnetic inverse scattering: Retrievable information and measurement strategies[END_REF]. En conséquence, de telles bases sont intrinsèquement capables de conduire à aborder le problème de diraction inverse dans un cadre multi-échelle, ce qui s'est révélé être un moyen ecace d'améliorer les résultats d'imagerie [START_REF] Caorsi | A new methodology based on an iterative multiscaling for microwave imaging[END_REF][START_REF] Bucci | Subsurface inverse scattering problems: quantifying, qualifying, and achieving the available information[END_REF].

Une grande variété d'études considère dorénavant l'information structurelle des solutions an de faciliter leur estimation. Par exemple, la structure dite de group sparsity ou de parcimonie structurée [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF] a été mise en ÷uvre dans de nombreuses applications, où les composantes des solutions sont susceptibles d'être soit toutes nulles soit non nulles dans un groupe. Ainsi, on vise à diminuer la dispersion an d'améliorer la solution en tenant compte de ce regroupement préalable. Dans ce cadre, nous nous sommes concentrés sur la parcimonie conjointe, ce qui constitue un cas particulier de la parcimonie structurée.

Nous avons proposé une approche d'inversion en deux étapes an de résoudre un problème inverse non linéaire en appliquant une parcimonie commune pour obtenir les sources équivalentes puis le contraste inconnu. Plus précisément, la parcimonie conjointe signie viii Résumé en Français que des vecteurs parcimonieux inconnus partagent le même support non nul [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF][START_REF] Fan | Enhanced joint sparsity via iterative support detection[END_REF].

Pour conclure, deux algorithmes d'inversion sont proposés pour une imagerie micro-onde ecace et précise de domaines d'investigation parcimonieux dans ce manuscrit. Le premier dépend du seuil de contraction douce avec des contraintes diérentes et le second exploite une parcimonie conjointe an de retrouver les propriétés électromagnétiques inconnues des diracteurs. Nous avons travaillé sur diérents diracteurs an de montrer la abilité des méthodes et nous avons proposé des résultats numériques qui illustrent leur ecacité, tout en proposant des améliorations complémentaires dans certains cas. ix List of Figures

General Context

This thesis presents my research work in the area of microwave imaging which has been attracting momentous research interests due to its potential as a convenient and ecient technique for medical systems [1,[START_REF] Henriksson | Quantitative microwave imaging for breast cancer detection using a planar 2.45 GHz system[END_REF][START_REF] Azghani | Microwave medical imaging based on sparsity and an iterative method with adaptive thresholding[END_REF], material characterization [3], subsurface probing, remote sensing, and non-destructive testing and evaluation [4,[START_REF] Pastorino | Microwave Imaging[END_REF]. The main goal of microwave imaging is to retrieve the distribution of the dielectric properties in a region of interest, with usual options of quantitative imaging (wherein the values of the electromagnetic parameters are sought within this region) and qualitative imaging (wherein what in ne matters is to identify certain zones in which such properties dier from those of the embedding medium).

It is well known that the electromagnetic inverse scattering problem is the basic formulation for microwave imaging methods. Although this problem has been suitably studied from a theoretical point of view, new methods are continuously developed to face new and challenging applications. To inspect dielectric targets, several techniques can be adopted [START_REF] Abubakar | Contrast source inversion method: State of art[END_REF][START_REF] Caorsi | A new methodology based on an iterative multiscaling for microwave imaging[END_REF]4953]. There are diverse research studies in this area and we only mentioned few of them here.

The need of ecient reconstruction methods and techniques has widely emerged for solving inverse electromagnetic scattering problems arising in microwave imaging in order to mitigate theoretical and practical diculties. High demand of such methods in various applications enforces the importance and the development of eective and accurate methods. However, implementation of stable, reliable, and ecient reconstruction algorithms is challenging because of the nonlinearity of the scattering equations and ill-posedness of the problem [4,[START_REF] Di Benedetto | Numerical linear algebra for nonlinear microwave imaging[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF].

In mathematics, there are often two problems which are opposite to each other. One of them is called direct (forward ) and the other one is called inverse problem. In the forward problem, the observable state of a system is determined by using all needed parameters. However, the inverse problem is about the precondition which creates observed data.

Mathematically, we have y = K(x), where x ∈ X is the unknown and y ∈ Y is the vector of observations. Usually, K() is a well-behaved continuous operator, and the solution of the forward problem (nd y given x) meets no signicant obstacles. On the other hand, the inverse mapping from y to x in the problems of interest is not easy to deal with. Normally, inverse problems cause mathematical models not to be well-posed in the sense of Hadamard such as no unique solution, no solution existence and non dependence on the data continuously [START_REF] Di Benedetto | Numerical linear algebra for nonlinear microwave imaging[END_REF][START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF]. These problems are called ill-posed and they lead to many numerical diculties as it makes most numerical algorithms unstable under data perturbations [4,[START_REF] Hanke | A convergence analysis of the Landweber iteration for nonlinear ill-posed problems[END_REF]. Numerical methods that can deal with these problems are called regularization methods.

Regularization is used to solve ill-posed problems by incorporating a-priori knowledge about x to stabilize the problem and to provide reasonable and useful solutions. For example, if it is known that the solution should be a discretization of a continuous function, this knowledge allows us to discard the wildest looking candidates, and to considerably reduce the set of possible solutions. The task is to minimize some measure J 1 (x) of proximity of y (observation data), as well as to satisfy as much as possible the a-priori information about x (the unknown parameter), by minimizing some appropriate measure J 2 (x).

The two objectives typically cannot be both reached at the same time, so we need a compromise, which can be simply obtained by taking a linear combination of the two:

J(x) = J 1 (x) + αJ 2 (x) (1.1)
Scalar α is the regularization parameter balancing the tradeo between the delity to the data, J 1 (x), and the delity to the prior information, J 2 (x). There is a whole family of solutions indexed by α, with the non-regularized (least squares) solution if α = 0, and a solution strongly favoring the a-priori information when α is large. In general, choosing an appropriate α is problem-dependent, and is a nontrivial task. With an appropriate choice for J(x), regularization eectively deals with all the three aspects of ill-posedness. Also, proper choice of J 2 (x) deals with lack of uniqueness and can dramatically reduce sensitivity to noise (improve the Lipschitz constant of the inverse function), making it continuous enough for practical applications.

The selection of a proper regularizer intimately depends on the property of (x) that one wishes to enforce, and that depends on the particular application. In many mathematical inverse problems, priors of choice are dierent forms of smoothness or constraints, and the corresponding regularizers are the l 2 norms of (x) or its derivatives. Sparsity prior is useful 1. Introduction 1.1. General Context when signals x that we look for have to be sparse. Sparsity of a vector (x) can be dened by the presence of a number of non-zero elements and zeros elsewhere, with respect to an appropriate basis, introducing the notion of L-sparse signal wherein this number is L. An appropriate numerical measure of sparsity is the count of non-zero elements. The base for this discussion is the work of Mallat [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], Donoho [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF], Rao [START_REF] Gorodnitsky | Sparse signal reconstruction from limited data using focuss: A re-weighted minimum norm algorithm[END_REF] and others on function approximation and optimal basis selection.

The problem of choosing an appropriate basis for a family of signals has received a great deal of attention over the past decade, and many new bases were introduced, such as wavelet bases, ridgelets, and curvelets, among many others [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF]. Despite the fact that any minimal spanning basis for a nite-dimensional space can represent perfectly any signal in the space, when only a subset of possible signals is of interest, some bases have better representational properties than others. Some applications which greatly benet from sparsity of representation are signal compression, denoising, and parameter estimation.

In compression for information transmission, if the representation of the signal is not sparse then we need to transmit the whole signal. However, if under a change of basis the representation becomes sparse, then substantial savings are possible. Most coecients of the representation are very small and if we set them to zero the perceptual quality of the signal will be aected very little. Hence, we are left with transmitting only the large coecients, which are few in number. This idea found use in commercial compression algorithms.

Another application where sparsity plays a key role is denoising. If the signal is sparse then separating it from the noise requires considerably less eort than when signal power is evenly distributed along the support of the signal. Therefore, for the purpose of facility of denoising of a class of signals, it is worthwhile to nd a basis in which the representation of all signals belonging to this class is as sparse as possible.

An inverse problem is often formulated in order to compute an approximation to a solution of the operator equation such as

K(x) = y, (1.2)
where K is an ill-posed operator and with the case of noisy data y δ yy δ ≤ δ

(1.3)
is available. For the stable approximation of a solution (1.2), a sparsity regularization method is used to minimize the functional

Θ(x) := 1 2 K(x) -y δ 2 + α k∈Λ w k | x, ϕ k | p (1 ≤ p ≤ 2) (1.4)
or more generally

Θ(x) := F (K(x), y δ ) + α k∈Λ w k | x, ϕ k | p (1.5)
where α > 0 is a regularization parameter, ϕ k k∈Λ is a basis and w k ≥ w min > 0, ∀k is a weight parameter. The functional F (K(x), y δ ) measures the error between K(x) and y δ .

If we let J 1 (x) = 1 2 K(x) -y δ 2
, and J 2 (x) = k∈Λ w k | x, ϕ k | p , then we have nothing but a regularized inverse problem of the form in (1.1). When p = 1 this method is called basis pursuit [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] (or LASSO [START_REF] Tibshirani | Regression shrinkage and selection via the LASSO[END_REF] in the statistical literature). The prior term, J 2 (x) has an eect of enforcing sparsity. Figure 1.1 gives some insight into why p-regularization with p ≤ 1 favors sparse x. In Fig. 1.1b, we show the level sets of p norms to the p-th power ( x p ) for p = 0.5, p = 1, and p = 2 of a two-dimensional vector. For a xed l 2 -norm, i.e, for all vectors that lie on a circle with xed radius, p norms with p ≤ 1 are minimized on the coordinate axes, i.e. preferring that some of the coecients are set exactly to zero, while others are large. In other words, p norms with p ≤ 1 maintains sparse solutions. This argument can be generalized to vectors in higher dimensions. Figure 1.1a shows l p norms for the same p's in one dimension. It shows that the penalty on large features (large x i ) is less for smaller p. Strong features are penalized much less severely in l p penalization with p ≤ 1 than in l 2 penalization (Tikhonov regularization). This motivates the smoothing eect of l 2 -penalization, and the feature-preserving behavior of p for p ≤ 1. It is well known that l 2 -norm regularizers are known to promote the smoothness in the solution and hence do not eciently produce accurate solutions when applied in domains with sharp variations, discontinuities, or sparse content (i.e., scatterers occupy much smaller volumes/areas in comparison to the whole investigation domain) [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. Such domains exist in many practical applications, such as see-through-the-wall imaging, hydrocarbon reservoir detection, radar imaging, and crack detection.

Another observation from Figure 1.1 is that l 1 -norm is convex, whereas when p < 1, l p -norm is no longer convex. The computational complexity for the minimization of some non-convex cost functions (p in particular) can be improved by using the half-quadratic regularization method [54]. The key idea is to introduce a supplementary vector s, and an extended cost function, Q(x, s), which is quadratic in x for a xed s, and argmin s Q(x, s) = J(x), for any x. If Q(x, s) is also easy to minimize in s (or even better if there is a closed- Sparsity regularization has been the foremost research area for the last years. The well-posedness and convergence rates of the method have been analyzed for linear inverse problems in [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] and for nonlinear inverse problems in [START_REF] Grasmair | Sparse regularization with lp penalty term[END_REF]. It has been shown that sparsity regularization is simple for use and very eective for inverse problems with sparse solutions. This method has been applied to Compressive Imaging and Electrical Impedance Tomography (EIT) problems [3032]. This regularization method leads sparse minimizers of (1.4) for p = 1 and promotes sparsity for 1 < p < 2 [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF].

In this manuscript we have excluded p = 0 since l 0 is a pseudo-norm (the distance property of norms is not satised), which counts the number of non-zero elements. Even if sparsity is most favored when p = 0, the above optimisation problem will become NPhard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF]. In other words, this creates a problem of multiple inma and solving it for the optimal global minima requires an intensive search that is computationally a burden [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. To avoid the occurrence of multiple inma, the l 1 -norm (when p ≤ 1) is used as the best convex approximation to the l 0 -norm optimization problem. This concept is known as convex relaxation [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF] as an alternative to the NP-hard problem. Furthermore, the l 1norm has been veried to elect a solution with excellent sparsity features. Therefore, we choose l 1 penalty to promote a-priori knowledge of the sparse representation throughout the work in this thesis.

Introduction

The use of l 1 -norms to achieve sparsity has been known for almost a decade as already mentioned. The Least Absolute Shrinkage and Selection Operator (LASSO) has been introduced in the statistics literature [START_REF] Tibshirani | Regression shrinkage and selection via the LASSO[END_REF], and Basis Pursuit algorithm [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF] for choosing a sparse basis has been proposed in the signal representation community at around the same time. The most important advantages of l 1 penalization schemes are their convexity, and the strong sparsity of the results (most indices of the result are set exactly to zero). Dierent versions of l 1 -regularization problems can be reformulated as linear, convex constrained or unconstrained quadratic, or second order cone (SOC) programming, all of which allow ecient and globally convergent algorithms. Another signicant benet of using l 1 penalization is a number of theoretical (e.g. [START_REF] Donoho | Uncertainty principles and ideal atomic decomposition[END_REF]) results showing that under certain sparsity conditions on the underlying unknown signal, the signal can be recovered exactly. This is quite surprising since the direct combinatorial formulation of the problem requires comparing solutions with all possible permutations of non-zero indices, which is very hard.

As a matter of fact, the unknown function (e.g., the contrast function or the equivalent current) can be considered as sparse (or compressible) in some suitable basis, due to the fact it may occupy a small portion of an otherwise known scenario or that may represent an almost homogeneous extended anomaly hosted in a background medium. In other words, a scatterer can be sparse with respect to a basis and not sparse versus another one. Moreover, in many practical cases, only a limited number of measurements is allowed, so that proper processing tools to handle this kind of complexity are needed.

One of the aims of this thesis is to extend the application of the microwave imaging paradigm to nonpixel basis functions, thus enlarging the set of retrievable scatterer distributions to targets of arbitrary shape and size provided that they are sparse with respect to a suitable chosen basis. Toward this end, wavelet basis functions will be used.

The wavelet transform stores information on both the spatial variation and the frequency content of the processed function, by representing it as a weighted sum of dilated and translated version of the so-called mother wavelet [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. In particular, large scale basis functions encode the low frequency content of the processed function, while small scale ones account for high frequency or ne details. Such an encoding can be replicated by considering mother wavelets at dierent scales, thus decomposing the functions into a coarse approximation and detail coecients at dierent levels [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Moreover, dierently from other representation bases, such as Fourier harmonics, the wavelet expansion enables a selective allocation of the unknown coecients within the investigated domain, thus allowing to concentrate them only where the object is located 1. Introduction 1.2. Related works among others and to consider a ner texture only where needed. It allows to reduce the number of unknown parameters, as compared to the usually exploited pixel representation, thus improving inversion reliability [START_REF] Isernia | On the local minima in a tomographic imaging technique[END_REF][START_REF] Bucci | Electromagnetic inverse scattering: Retrievable information and measurement strategies[END_REF]. As such, they are intrinsically capable of approaching the inverse scattering problem within a multiscale framework, which has been shown to be an eective way to improve imaging results [START_REF] Caorsi | A new methodology based on an iterative multiscaling for microwave imaging[END_REF][START_REF] Bucci | Subsurface inverse scattering problems: quantifying, qualifying, and achieving the available information[END_REF].

After years of intensive research on l 1 -problem solving, it would appear that most relevant algorithmic ideas have been either tried or, in many cases, re-discovered. Yet interestingly, additionally to the existing methods we also studied the classic idea of soft shrinkage in the application of nonlinear inverse problem (considering wavelet domain as well) where the inverse problem is solved without linearization and the application of the joint sparsity within microwave imaging.

Related works among others

As stated in above, inverse electromagnetic scattering problems have two fundamental properties that make their solution a challenging task, the rst one is the nonlinearity of the scattering equations [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] and the second one is the ill-posedness of the problem [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. The scattered eld is an integral function of the total eld induced inside the investigation domain multiplied by the contrast (the dierence between the permittivities of the investigation domain and the background medium). The induced total eld is also a nonlinear function of the permittivity [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. Consequently reconstruction of the permittivity from scattered elds is a nonlinear inverse problem.

Additionally, the fact that the scattered eld is an integral function of the total eld and the permittivity inside the investigation domain leads to an ill-posed inverse problem [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. The integral operator suppresses the contribution of the fast-varying components of the permittivity to the scattered eld and makes them dicult to retrieve from the measurements [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF]. This physical nature of the scattering phenomena when combined with the fact that only a nite set of samples, possibly not at optimal location with respect to the sought scatterers and possibly as well not suciently independent from one another, can be collected by the measurements, which overall does not provide sucient information. These samples are almost always corrupted by noise, and a discretized approximate model has to replace the exact mathematical description of the problem. Thus, this makes the inverse problem severely ill-posed [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF].

The strength of the nonlinearity in the scattering problem increases with the strength 1. Introduction of the scattered elds due to the object(s) in the investigation domain. Hence, higher permittivity (this notion being size-dependent as well) means stronger non-linearity. Most of the deterministic methods shown to be eective in solving electromagnetic inverse problems involving weak scatterers make use of linear rst-order approximations [START_REF] Devaney | Mathematical Foundations of Imaging, Tomography and Waveeld inversion[END_REF], such as diraction tomography [START_REF] Zorgati | Eddy current testing of anomalies in conductive materials. i. qualitative imaging via diraction tomography techniques[END_REF], rst-order Born approximation [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Rajan | A comparison between the Born and Rytov approximations for the inverse backscattering problem[END_REF], Kirchho approximation [START_REF] Pastorino | Microwave Imaging[END_REF] and Rytov approximation [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Rajan | A comparison between the Born and Rytov approximations for the inverse backscattering problem[END_REF].

Even though these methods are computationally less demanding, they fail to provide accurate solutions when strong scatterers are present in the domain being investigated. In such cases, where the nonlinearity is strong, linearization does not produce an accurate approximation of the nonlinear problem and one needs to use more rigorous techniques such as inexact Newton nonlinear solver combined with either contrast-eld or contrastsource formulations [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Bozza | An inexact Newton-based approach to microwave imaging within the contrast source formulation[END_REF], nonlinear conjugate gradient method [START_REF] Abubakar | Contrast source inversion method: State of art[END_REF], nonlinear steepest descent (e.g. nonlinear Landweber) algorithm [START_REF] Li | Image reconstruction by nonlinear Landweber iteration for complicated distributions[END_REF][START_REF] Hettlich | The Landweber iteration applied to inverse conductive scattering problems[END_REF], and Levenberg-Marquardt scheme [START_REF] Franchois | Microwave imaging-complex permittivity reconstruction with a Levenberg-Marquardt method[END_REF], and distorted Born methods [START_REF] Chew | Reconstruction of two-dimensional permittivity distribution using the distorted Born iterative method[END_REF]. These methods are computationally more demanding when compared to those that make use of the rst order linearization schemes. Besides these deterministic methods, there are also nonlinear approaches based on stochastic (multiple/single agent) techniques such as simulated annealing [START_REF] Van Laarhoven | Simulated annealing[END_REF], genetic algorithms [START_REF] Chiu | Image reconstruction of a perfectly conducting cylinder by the genetic algorithm[END_REF], dierential evolution [START_REF] Qin | Dierential evolution algorithm with strategy adaptation for global numerical optimization[END_REF], etc. Other methods that benet from the best of the both worlds have also been developed. These make use of higher order linearization schemes or iterative application of the rst-order ones. Examples include extended Born approximation [START_REF] Zhang | Two nonlinear inverse methods for electromagnetic induction measurements[END_REF] and the higher-order Born approximations [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Estatico | An inexact-Newton method for short-range microwave imaging within the second-order Born approximation[END_REF], Born iterative method (BIM) [START_REF] Wang | An iterative solution of the two-dimensional electromagnetic inverse scattering problem[END_REF] and the variational Born iterative method [START_REF] Zaiping | Variational Born iteration method and its applications to hybrid inversion[END_REF], respectively. The methods in this last group are computationally more ecient when compared to nonlinear inversion methods but also applicable for a wider range of scatterer strength when compared to the rst-order linearization methods.

Regardless of the method used for tackling the nonlinearity, the ill-posedness of the inverse scattering problem should be accounted for [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF][START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. For this purpose, linear and nonlinear regularization methods, which minimize a cost function weighted between measurement/data mist and a penalty term, have been developed. The most popular choice of the penalty term is the l 2 -norm of the solution [START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. The resulting minimization problem can be solved using the well-known Tikhonov scheme [START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. Additionally, using truncated Landweber or conjugate gradient iterations leads to a similar type of regularization [START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Tveito | Parameter estimation and inverse problems[END_REF]. In linear programming, all these methods eectively lter singular values of the discretized minimization problem, which are smaller than a specied threshold level or a regularization parameter to alleviate its ill-posedness.

As stated before, one common alternative is the convex relaxation of this combinatorial problem which consists in replacing the l 0 -norm by the convex l 1 -norm to take advantage of well-known algorithms in convex optimization. The basis pursuit (BP) principle or its denoising (BPDN) adaptation, proposed in [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], deals with the l 1 relaxation of problems using linear and quadratic programming algorithms (e.g., simplex, active set, or interiorpoint methods). Iterative thresholding algorithms (ISTA) [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF] and its accelerated version FISTA [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF], belonging to the family of forward-backward algorithms [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], are now well known to be very ecient for such l 1 relaxed problems. The fact that this convex relaxation provides sparse solutions comes from the singularity of the l 1 -norm at zero [START_REF] Fan | Variable selection via nonconcave penalized likelihood and its oracle properties[END_REF].

On the other hand, the Compressive Sensing (CS) theory has been introduced for models relating the data to the unknowns that are linear [START_REF] Donoho | Compressed sensing[END_REF]. As a consequence, a number of dierent CS toolboxes are available to tackle linear problems, whereas the development of such tools for the nonlinear case is ongoing and not yet fully assessed [START_REF] Ohlsson | Nonlinear basis pursuit[END_REF][START_REF] Blumensath | Gradient pursuit for non-linear sparse signal modelling[END_REF]. For such a reason, the application of CS to inverse scattering problems has been rst explored in those cases that can be handled through linear scattering approximations. These include the Born [4,[START_REF] Ambrosanio | A compressive-sensing-based approach for the detection and characterization of buried objects[END_REF][START_REF] Bevacqua | A compressive sensing approach for 3D breast cancer microwave imaging with magnetic nanoparticles as contrast agent[END_REF] and Rytov [START_REF] Oliveri | Bayesian compressive optical imaging within the Rytov approximation[END_REF] approximations for weak scatterers, as well as the approximation based on the recently proposed virtual experiments framework [START_REF] Bevacqua | Microwave imaging of nonweak targets via compressive sensing and virtual experiments[END_REF][START_REF] Bevacqua | Exploiting sparsity and eld conditioning in subsurface microwave imaging of nonweak buried targets[END_REF].

However, inverse scattering problems are non-linear, so that such an issue has to be addressed to gain a full advantage of CS in this framework. To this end, several authors are developing methods to take advantage of CS within inverse scattering approaches that are not limited to the range of validity of linear approximations.

For instance, the authors of [START_REF] Poli | MTBCS-based microwave imaging approach through minimum-norm current expansion[END_REF][START_REF] Oliveri | Bayesian compressive optical imaging within the Rytov approximation[END_REF] have proposed an approach that is based on the solution of a number of linear inverse source problems. Such an approach does not consider approximations in modelling the non-linear scattering interactions and it can be applied to all cases in which the contrast source is sparse.

Another possible way to extend the applicability of CS for the solution of inverse scattering problems is to exploit distorted wave methods or distorted-iterated methods, in which the solution is iteratively achieved through a succession of linear inversion steps [START_REF] Azghani | Microwave medical imaging based on sparsity and an iterative method with adaptive thresholding[END_REF]6769]. Because of the linearity of the models considered at each step, these latter methods naturally lend themselves to be paired with CS. CS-based method based on l 1 norm minimization [START_REF] Donoho | Compressed sensing[END_REF] was used to nd the sparse solution in a sparse domain. In [START_REF] Poli | MTBCS-based microwave imaging approach through minimum-norm current expansion[END_REF][START_REF] Oliveri | Bayesian compressive optical imaging within the Rytov approximation[END_REF], a method based on Bayesian framework [START_REF] Gharsalli | Inverse scattering in a Bayesian framework: application to microwave imaging for breast cancer detection[END_REF] in compressive sensing [START_REF] Ji | Bayesian compressive sensing[END_REF] was utilized to solve an inverse scattering problem. The Bayesian framework has been suggested for some time for microwave medical imaging [START_REF] Gharsalli | Inverse scattering in a Bayesian framework: application to microwave imaging for breast cancer detection[END_REF]. Compared with the traditional l 1 norm minimization [START_REF] Donoho | Compressed sensing[END_REF], BCS searches for the sparse solution from a Bayesian probability 1. Introduction perspective.

Two popular non-quadratic cost functions are total variation and entropy [START_REF] Karl | Regularization in image restoration and reconstruction[END_REF]. Total variation puts a penalty on the sum of variations of the signal J 2 = Cx 1 , where C is a discrete approximation to the gradient operator. Total variation is most frequently used in image processing applications, such as image restoration. In comparison to the Tikhonov regularization with L = C, the penalty on strong features is less severe, and the reconstruction can contain sharp edges. It works very well in practice with images that can be described as piecewise-smooth. Even if the regularizing function with l 1 -penalty is similar with the total-variation we take the l 1 norm of the values of x instead of their derivatives. Total variation allows sparse jumps of the gradient of x, whereas the l 1 penalty favors sparse values of x. To lower the penalty on strong features even further, several non-convex functions have also found use [START_REF] Chartrand | Exact reconstruction of sparse signals via nonconvex minimization[END_REF].

It is also worthwhile to mention that unlike inexact Newton [START_REF] Desmal | A preconditioned inexact Newton method for nonlinear sparse electromagnetic imaging[END_REF] and Born iterative [START_REF] Desmal | Shrinkage-thresholding enhanced Born iterative method for solving 2d inverse electromagnetic scattering problem[END_REF][START_REF] Bagci | Sparsity-regularized born iterations for electromagnetic inverse scattering[END_REF] methods with sparsity constraints, the proposed scheme in Chapter 3 avoids generation of a sequence of linear sparse optimization problems and requires only one regularization parameter, which directly penalizes the nonlinear problem, to be set. Consequently, it simplies the task of heuristic parameter tweaking, which is oftentimes very cumbersome for existing inversion algorithms.

Outline of thesis

While describing the contents of the thesis chapter by chapter, we briey summarize our main contributions. The rst major contribution is the development of a sparse signal reconstruction framework for microwave imaging. In this framework we formulate various optimization problems for microwave imaging scenarios. We adapt and use two paradigms for the numerical solution of the optimization problems. Finally, we carry out an extensive performance analysis of the proposed methods.

Chapter 2: Problem Statement

In this chapter, we presented the formulation of the direct problem modeling where incident wave and object are known by using a Method of Moments as discussed in [START_REF] Richmond | Scattering by a dielectric cylinder of abritrary cross-section shape[END_REF].

We show numerical results in order to illustrate the sensitivity of the model to various noise levels and how the same level of noise has dierent randomization. Next, we introduce our nonlinear inversion method in wavelet domain so that it takes the advantage of the adaptive multiresolution features of the wavelet basis to accommodate the trade o between spatial resolution and inversion stability. Specically, we represent the unknown contrast sources in terms of the wavelet basis functions to reduce the number of nonzero coecients whereas this proposed approach augments the work that has been done previously in this research area.

Chapter 4: Two-Step Inversion Method

This chapter is another contribution of our thesis. We address the analysis of the joint sparsity, a special case of group sparsity, which gives multiple sparse solutions that share a common nonzero support. We oer numerical results which demonstrate the potential of the proposed two-step inversion approach both in 2D and 3D. Modeling the associated direct problem in electromagnetic scattering is an initial yet essential step to describe how the parameters of the model are translated into observable eects. The direct problem consists in modeling the physics of the interaction between a known interrogating wave and the object which is supposed to be known as well. This interaction is described here by the Helmholtz wave equation. Applying the Green's theorem to this equation and taking into account the conditions of continuity of the elds and of radiation at innity [START_REF] Chew | Waves and Fields in Inhomogeneous Media[END_REF], we are led to an integral representation of the electric eld consisting of two coupled integral equations, observation and coupling equation (or state), respectively. The solution of the direct problem that requires the discrete counterparts of these integral equations which are obtained in an algebraic framework using the method of moments [START_REF] Richmond | Scattering by a dielectric cylinder of abritrary cross-section shape[END_REF][START_REF] Gibson | The Method of Moments in Electromagnetics[END_REF].

Validation of the direct model is a necessary step before switching to inversion. This is to verify that the direct model well describes the phenomena (in our case the scattered elds) to be observed in a controlled situation. Note that, in this thesis, we have no experimental data. This is why, when solving the inverse problem, we will generate synthetic data using a dierent direct model than the one used during the inversion. Therefore, we will avoid committing an `inverse crime' in the sense of [START_REF] Colton | Inverse Acoustic and Electromagnetic Scattering Theory[END_REF] which would consist in testing the inversion algorithm on the data obtained by means of a model closely related to the one used in the inversion. The validation of the direct model will thus consist in verifying that the dierences between the scattered elds produced by these two models remain weak.

Formulation

In the following we consider homogeneous, non-magnetic and isotropic media and we place ourselves in a 2D imaging model in a transverse magnetic conguration (TM) where we consider that the object is of innite extension and invariant according to one of its axes and is illuminated by a wave whose electric eld (E inc (r) = E inc z) is polarized parallel to this axis (perpendicular to the plane of the gure), which generates a scattered eld perpendicular to the plane of the gure.

We now establish the equation of propagation of the electromagnetic wave on the basis of the Maxwell equations which constitute the basis of electromagnetism and which are written at an observation point r and at time t such as

∇ • D(r, t) = ρ(r, t) (2.1) ∇ × E(r, t) = - ∂B(r, t) ∂t (2.2) ∇ × H(r, t) = J (r, t) + ∂D(r, t) ∂t (2.3) ∇ • B(r, t) = 0, (2.4)
where D(r, t)

is the electric displacement eld in [Cb/m 2 ], E(r, t) is the electric eld in V /m, H(r, t) is the magnetic intensity or magnetic eld in [A/m], J (r, t) is the current density in [A/m 2 ] and B(r, t) is the magnetic induction eld in [W b/m 2 ], t ∈ R, r = (x 1 , x 2 ,
x 3 ), ∇× denotes curl operator and ∇• denotes divergence operator. To solve the system of equations (2.2)-(2.4), one needs to introduce two additional equations, i.e., the so-called constitutive relations. Such equations enable to express the magnetic ux density B(r, t) and the electric displacement D(r, t) as a function of the electromagnetic parameters of the medium where Maxwell equations are set. The constitutive relations read as

B(r, t) = µH(r, t) (2.5) D(r, t) = E(r, t), (2.6)
where is the dielectric permittivity and µ is the magnetic permeability of the medium at r. We now nd ourselves in the framework of a time-harmonic regime where a source emits a wave of pulsation ω and an implicit time dependence of the elds in e -iωt is chosen so that it will be omitted in the remainder of the discussion leading replacement of the derivations with respect to time appearing in the equations of Maxwell by a factor -iω.

To establish the equation of propagation in electric eld, we now place ourselves in a medium without sources or charges. From Maxwell's equations and taking into account constitutive relations, we arrive at

∇ × ∇ × E -ω 2 µ E = 0.
(2.7)

By developing the vector operator appearing in this equation, remembering that the electric eld has only one component along the z-axis that we write E and introducing the propagation constant k of the medium considered (k 2 = ω 2 µ), we obtain the equation of scalar Helmholtz waves veried by the electric eld

∇ 2 E(r) + k 2 E(r) = 0. (2.8)
We shall now suppose that the object which we are interested in is immersed in a homogeneous medium D and is contained in a test domain C (C ⊂ D). The dierent media are non-magnetic (magnetic permeability equal to the one of vacuum µ 0 = 1, 256×10 -6 Hm -1 ) and are characterized by their propagation constant such that k(r) 2 = ω 2 ε 0 ε r (r)µ 0 + iωµ 0 σ(r), where ε 0 and µ 0 are the permittivity and the permeability of air, respectively ε r (r) and σ(r) are the relative permittivity and conductivity of the medium as r ∈ D is an observation point. The dielectric properties of D are described by the inhomogeneous contrast function dened as

χ(r) = (k(r) 2 -k 2 B ),
where k 2 B = ω 2 ε 0 µ 0 is the propagation constant of the embedding medium D. We assume that sources and receivers are located at the positions r s and r r , respectively.

The scattered electric eld E di (r r , r s ) measured via a receiver placed at r r due to the incident wave emitted by a source placed at r s adheres to the following domain integral

equation [5] E di (r r , r s ) = D G(r r , r )χ(r )E(r , r s )dr (2.9)
with E(r, r s ) being the total electric eld induced within the object by the incident wave, and G(r, r ) is a Green`s function which represents the electromagnetic response to a line source radiating in free-space. In a homogeneous medium, in the case of two dimension that concerns us, it is given by:

G(r, r ) = -iωµ 0 4 H (1) 0 (k B r -r ) (2.10)
and H

(1) 0 is the zero-th order Hankel function of the 1st kind for the 2D case. Furthermore, E(r, r s ) is obtained as

E(r, r s ) = E inc (r, r s ) + D G(r, r )χ(r )E(r , r s )dr ∀r ∈ D.
(2.11)

The direct problem is dened as the calculation of E di (r r , r s ) by solving (2.9) and (2.11) when χ(r), G(r, r ) and E inc (r, r s ) are known whereas the inverse problem (or imaging problem) is dened as the determination of χ(r) within a prescribed domain D from the knowledge of E di (r r , r s ), G(r, r ) and E inc (r, r s ) for N s sources and N r receivers.

Formulation of Contrast Source

We now rewrite the observation and coupling equations for the contrast source which is induced inside the object by the incident wave such as

E di (r r , r s ) = D G(r r , r )J(r , r s )dr (2.12) and J(r, r s ) = J inc (r, r s ) + χ(r ) D G(r, r )J(r , r s )dr ∀r ∈ D.
(2.13) the contrast source being dened as J(r, r s ) = χ(r)E(r, r s ).

(2.14) Solving equations (2.12) and (2.13) is done from their discrete counterparts obtained by using the method of moments. The domain D containing the unknown object is discretized in N = N x × N y small square pixels so that the electric eld and the contrast can be considered as constants in each of them (Fig. 2.1).

+ ⇝ Source(s) Receiver(s) D Ω ε (r) , µ 0 ε 0 , µ 0 ŷ x ẑ r s

Discretization of the problem

The method of moments is a method commonly used to transform functional equations into matrix equations by projection on sets of basis functions and of test functions [START_REF] Richmond | Scattering by a dielectric cylinder of abritrary cross-section shape[END_REF].

The problem being reduced under the form AX = B, where A is a matrix, X is a vector of unknown and B is a known vector, intuitively, when the matrix A is invertible, this allows us to get the desired solutions.

Let L be a linear operator L : G -→ H, and two functions g ∈ G and h ∈ H such as h = L(g). One seeks to determine g knowing L and h. The function g can be written in the form g = i=1,••• ,I g i u i , where u 1 , u 2 , • • • , u I are constant coecients which are a set of basis functions u i . The property of linearity of the operator L gives us

h = L(g) = L I i=1 g i u i = I i=1 g i Lu i (2.15)
The goal is to approach equation (2.15) by a discrete linear system. It denes a set of test functions t 1 , t 2 , • • • , t N on which projected the two members of this equation:

t n , h H = I i=1 g i t n , Lu i H , ∀n = 1, • • • , N, (2.16)
where ., . H is a scalar product on H such as t n , h H = H t n (x)h * (x)dx.

Formulation of the Direct Problem

Thus, by omitting the index H, one obtains the algebraic writing h = Lg with

L =       t 1 , Lu 1 t 1 , Lu 2 • • • t 1 , Lu I t 2 , Lu 1 t 2 , Lu 2 • • • t 2 , Lu I . . . . . . . . . . . . t N , Lu 1 t N , Lu 2 • • • t N , Lu I       , g =     g 1 . . . g I     , h =     t 1 , h . . . t N , h     ,
The choice of basic functions and test functions depends essentially upon the nature of the physical problem to be treated. Note that, if we choose test functions identical to the basis functions, we come across with the special case of the Galerkin method [START_REF] Harrington | Field Computation by Moment Methods[END_REF].

The two coupled equations (2.9) and (2.11) have a bilinearity property with respect to the two variables χ and E. Then, the method of moments can be applied to these equations. Indeed, if one considers the equations of observation and of state like two equations linear with respect to the variable E, equations (2.9) and (2.11) can be written as follows

E di (r r , r s ) = D G(r r , r )χ(r )E(r , r s )dr ⇒ E di = L 1 E (2.17) E inc (r, r s ) = E(r, r s ) - D G(r, r )χ(r )E(r , r s )dr ⇒ E inc = L 2 E, (2.18)
where L 1 and L 2 are linear operators.

Then, in order to apply the principle of the method of moments, we must choose the basis functions u 1i and u 2i and the test functions t 1i and t 1i in order to construct scalar products

D G(r r , r )χ(r)E(r, r s )dr , t 1i = L 1ij (2.19) E inc (r, r s ) = E(r, r s ) - D G(r, r )χ(r)E(r, r s )dr ⇒ E inc = L 2 E, (2.20)
Here, we will take basis functions as the characteristic functions of the pixels partitioning the domain D and test functions as distributions of Dirac δ located at the centers of each of the elementary pixels for the coupling equation and at the points of measurement for the observation equation.

Thanks to a classical method of moments using pulse basis/point matching method [START_REF] Richmond | Scattering by a dielectric cylinder of arbitrary cross section shape[END_REF] where the domain under test D is discretized into N = N x × N y pixels, a discretized version of the previous equations is obtained. In so doing χ(r) is approximated by using pulse basis functions dened over square cells:

χ(r) = N n=1 χ n C n (r), (2.21)
where C n (r) is the basis function dened as

C n (r) =    1, r ∈ D n 0, r / ∈ D n , (2.22)
in a `pixel' representation. Here, χ is a N × 1 vector storing the samples of χ(r). Additionally, the discretized version of (2.9) stands as

E di i = G or J i , i = 1, . . . , N s , (2.23)
where E di i and J i are complex vectors of size N r and N respectively and G or a complex matrix of size N r × N . The discretized version of (2.14) is then

J i = diag (χ)E i , i = 1, . . . , N s , (2.24)
where E i and χ are complex vectors of size N and diag (χ) is a diagonal matrix of size N × N obtained from χ. Finally the discretized version of (2.11) is

E i = E inc i + G oo J i , i = 1, . . . , N s (2.25)
where E inc i is a complex vector of size N and G oo a matrix of size N × N .

Thus, if we denote the center of the elementary pixel µ i (permeability) by r i and assuming that the contrast χ and the electric eld E are constant in each pixel µ i and equal to χ(r i ) and E(r i ) respectively, the elements of the coupling matrix are written as

L 2ij = δ ij -χ(r i ) µ j G(r i , r )dr , (2.26)
where δ ij is Kronecker delta. The numeric integration of the Green's function on the square cell ∆ i is analytical following [START_REF] Richmond | Scattering by a dielectric cylinder of abritrary cross-section shape[END_REF] where the square cell is approximated by a circular cell of the same surface leading to:

∆ i G(r j , r )dr =        1 k 2 B i 2 πk B RH (1) 1 (k B R) -1 , if i = j, i 2k B πRH (1) 0 (k B r ij ) J 1 (k B R) , if i = j, (2.27)
with r ij = |x j -x i | and R = ∆x∆y π ; H

(1) 1 being the Hankel function of rst order and rst kind and J 1 the Bessel function of the rst kind.

Validation of the model

Firstly our numerical code has been validated by comparison with numerical results provided by an analytical solution available for the case of a cylindrical obstacle illuminated by a line source for a single frequency and a single incidence.

The validation in Fig. 2.2 has been done with the frequency being 3 GHz and the radius of the cylinder equal to λ 0 /2 whereas λ 0 = 1. The relative permittivity of the object is 2 whereas it is 1 for the embedding medium. The source is located at (0, 2λ 0 ) while there is 36 receivers.

Generally the dierence between the model and the data is low but increases with frequency. This is explained by the fact that the discretization, i.e. pixel size, is constant whatever the frequency is. Thus, we have seen that a good convergence of the calculations carried out using the method of moments requires to index this size to the wavelength.

Comparison of the data with dierent models

In this section we studied an example where has the true value of the relative permittivity of the object is 2 whereas it is 1 for embedding medium. The scattering object which is a dielectric square sided λ (1 m) under test is contained in a l = 3 × λ-sided square investigation area D centered at the origin, and the discretization size is n × n = 36 × 36

for the forward problem. The number of transmitters and receivers located around the investigation area is 29. The frequency of the transmitters is 300 MHz. The measured eld samples are generated by adding either 10 dB or 20 dB noise with a zero mean additive Gaussian noise with unknown variance.

From Fig. 2.4 we notice that the outputs of the models coincide relatively well with the data, both in modulus and in phase. However, as we increase the level of noise, the discrepancy between the model and the data increases. On the other hand, Fig. 2.5

shows us that whenever we add noise of the same level to our model we get dierent randomization of the same noise level with respect to each other and also to the data.

This specied issue led us to be careful while choosing the model to solve the inverse problem due to the sensitivity of the data to noise.

Conclusion

In this chapter, we presented the modeling of the direct problem where incident wave and object are known. Primarily, we are interested in the calculation of the scattered eld resulting from their interaction. The resolution of the latter corresponds to a rst necessary step to solve the inverse problem in which it is a matter of reconstructing the unknown object, the incident wave and the scattered eld then being known. In the absence of experimental data, we have studied the electromagnetic eld calculations carried out using this forward model by comparison with synthetic results obtained using the forward model added with noise.

Methods to Solve Nonlinear

Inverse Problem

Motivation

In this chapter we propose two methods in order to solve the nonlinear inverse problem.

Initially, we propose a method where sparsity constraint is directly applied to the problem of reconstructing the complex internal dielectric properties of an object based on knowledge of the external scattered eld which is generated by the interaction between the object and a known incident eld. The nonlinear optimization problem is solved by the iterative algorithm of soft shrinkage in order to enforce the sparsity constraint. Sparsity is applied at each iteration by a soft thresholding function.

Since the electromagnetic inverse scattering problem which we consider is nonlinear and ill-posed, the reduction of the number of unknown parameters is important to alleviate the ill-posedness [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF][START_REF] Candes | Curvelets: A surprisingly eective nonadaptive representation for objects with edges[END_REF][START_REF] Li | A contrast source inversion method in the wavelet domain[END_REF]. As a matter of fact, wavelet expansion (as considered hereafter) should allow (this is in eect a topic of discussion in the present work, noticing that optimality of a basis is strongly linked to the way the scatterers at hand t it) to reduce the number of unknown parameters with respect to the usually exploited pixel representation so that it advances the inversion reliability [START_REF] Isernia | On the local minima in a tomographic imaging technique[END_REF].

Moreover, the adopted basis plays a key role as it has to be accurate and ecient.

Eciency requires a reduced number of coecients in the representation, while accuracy involves a low representation error, which quanties the mismatch between the actual function and its projection onto the considered basis. In this respect, this chapter argues a second approach which studies whether the wavelet basis is a good choice to accommodate for the trade o between eciency and accuracy of the representation or not in microwave imaging with simulation results provided.

Soft Shrinkage Method for Nonlinear equations

Soft shrinkage is an approach which minimizes a nonlinear Tikhonov functional with sparsity promoting penalty term. The algorithm is based on the iterated soft shrinkage approach originated for linear operators in the work [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]. A generalization to nonlinear inverse problems has been studied in [START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF][START_REF] Bonesky | A generalized conditional gradient method for nonlinear operator equations with sparsity constraints[END_REF].

The algorithm performs a gradient descent step which involves the adjoint gradient of the cost function with a step size τ , like Landweber method, and then a shrinkage step.

The latter enforces the sparsity of the reconstruction by setting the small coecients to zero. Following this, Barzilai and Borwein (BB) method is suggested in order to choose the step size to overcome the slow convergence of the iterative soft shrinkage algorithm with xed step size. The solution of the inverse problem can be obtained by minimizing the cost function which is the dierence, in other words called an error, between the measured quantity and the solution obtained by a forward problem. The cost function to be minimized is of the form

F (χ) = 1 2 ζ(χ) -E di 2 K(χ)
+α χ 1 .

(

In other words, by combining (2.23) and (2.25) and using J inc i = diag (χ)E inc i i = 1, . . . , N s the inverse problem can be rewritten as

χ = argmin χ ζ i -G or diag (χ) [I -diag (χ)G oo ] -1 J inc i 2 + α χ 1 , (3.2)
where ζ i is a vector of size N r which gathered the signal due to the source #i measured by the N r receivers. It can be seen from (3.2) that the inverse problem is nonlinear in χ and ill-posed. The dimensions of these matrices may become quickly important as soon as the domain D is more than a few wavelengths. Indeed, the convergence of the calculations carried out using the method of moments imposes an upper limit on the size of the pixels to be considered [START_REF] Richmond | Scattering by a dielectric cylinder of abritrary cross-section shape[END_REF][START_REF] Hagmann | Upper bound on cell size for momentmethod solutions[END_REF]. Pixel sides a ≤ λ/10 generally ensure a good convergence.

Iterative soft shrinkage is a standard approach for minimizing a functional involving the l 1 -penalty and it has the form as in (3.1) where ζ : X → Y is a bounded and nonlinear operator. The l 1 penalty can promote a-priori knowledge of the sparse representation. At rst, the algorithm is started by choosing an initial guess χ 1 , and the iteration continues 3. Methods to Solve Nonlinear Inverse Problem

3.2.1. Gradient K as χ k+1 = S α χ k -τ ζ * (χ k ) ζ(χ k ) -E di , (3.3
(S α (χ)) i =    (|χ i | -α)sign(χ i ), if |χ i | > α 0, otherwise. (3.4) The term ζ * (χ k ) ζ(χ k ) -E di is the gradient of the discrepancy 1 2 ζ(χ) -E di 2 . 3.2.1. Gradient K
The gradient is obtained by using the adjoint method. The main idea is to obtain the gradient (Fréchet derivative) of the whole discrepancy term and avoid calculating ζ * (χ) in the iteration rule (3.3). For the calculation of K(χ) we use the adjoint method which is presented in [54] to reduce the calculation costs. This has been achieved by solving the adjoint problem as in A.1. For simplicity, at rst it can be again stated with simpler notations that the least squared cost function describing the goodness of t of a hypothesized case to a measured case can be written as

K(χ) = 1 2 ζ (χ) -E di 2 = 1 2 ζ (χ) -E di , ζ (χ) -E di (3.5)
If χ is perturbed by a small amount δχ, then the K(χ) changes according to

K (χ + δχ) = K (χ) + δχ, ζ (χ) * ζ (χ) -E di + δχ, ζ (χ) * (ζ (χ) -E di ) + O (δχ) . (3.6)
By denition of the Fréchet derivative,

∇χK = ζ (χ) * (ζ (χ) -E di ), (3.7)
where the overlined terms are the complex conjugate ones and gradient

K (χ) = ∇ χ K (χ)| χ .
The detailed work can be found in A.1 as stated above.

Smooth Gradient K s

It is emerged in practice that the gradient K (χ) has unnatural oscillating properties which can be avoided by using the smoother gradient. This process is also called denoising [54]. Therefore, we look for a Sobolev smoothed gradient K s (χ). For instance, K (χ)ζ is the image of ζ ∈ D under application of the adjoint operator considered as an operator mapping from D into P = L 2 (Ω). Its image denoted by K s (χ) under the adjoint operator with respect to the newly dened weighted inner product mapping into the smaller space

P such as K (χ)x, ζ D = x, K (χ)ζ P = x, K s (χ)ζ P (3.8) following K s (χ)ζ = (δI -β∆) -1 K (χ)ζ, (3.9)
where I refers to the identity, and ∆ refers to the Laplacian operator. A proper choice of the weighting parameters δ and β will allow us to drive the regularization properties of our algorithm in an ecient and predictable way. We should choose β close to one in order not to lose the dierentiability for the next iteration. For theoretical justication we refer to [START_REF] Jin | A reconstruction algorithm for electrical impedance tomography based on sparsity regularization[END_REF]. Applying Green`s formula to the right hand side of equation (3.8) yields equation (3.9). Following this, we can rewrite the iteration formula with this property following as

χ k+1 = S α (χ k -τ K s (χ k )).
(3.10)

The Step Size

The step size τ can be determined in order to fasten the algorithm. The motivation for increasing the rate of convergence is the comparison with the classical Landweber iteration whose slow convergence results from using a constant step size which is very small. Therefore, we select the step size in a way to increase the convergence speed where we consider only the steepest descent operation χ kτ K s (χ k ) of the algorithm. The selection is done by the two-point rule of Barzilai and Borwein (BB) which calculates the step size as

τ k = arg min τ τ (χ k -χ k-1 ) -(K s (χ k ) -K s (χ k-1 )) 2 .
(3.11)

The choice of this functional is motivated by the secant equation

K (χ k ) = K χ k-1 + B(χ k -χ k-1
) with τ times the identity operator as the approximation of the Hessian B and describes the approximation performance of the last iteration step. This equation does not necessarily have a solution, so it is solved in a least-squared sense. In one-dimensional real case of this procedure implies the secant method. By minimizing the functional we get the following formula for the step size

τ k = χ k -χ k-1 , χ k -χ k-1 χ k -χ k-1 , K s (χ k ) -K s (χ k-1 )
.

(3.12)

Another approach, also from BB, is that τ times the identity operator imitates the inverse of the Hessian B over the last step. This results in

τ k = arg min τ (χ k -χ k-1 ) -τ (K s (χ k ) -K s (χ k-1 ) (3.13)
and therefore we get [START_REF] Jin | Sparsity regularization for parameter identication problems[END_REF][START_REF] Jin | A reconstruction algorithm for electrical impedance tomography based on sparsity regularization[END_REF]

τ k = χ k -χ k-1 , K s χ k -K s χ k-1 K s (χ k ) -K s (χ k-1 ) , K s (χ k ) -K s (χ k-1 ) (3.14) with K s χ k = (δI -β∆) -1 ζ * (χ k ) ζ(χ k ) -E di . (3.15)
In the implementation we use this step size as an initial guess and it is decreased geometrically until the Armijo condition [START_REF] Armijo | Minimization of functions having Lipschitz continuous rst partial derivatives[END_REF] is satised. In general, descent methods determine a descent direction d k with a step size α k ∈ (0, 1] by an inexact line search such as Armijo, Wolfe or Goldstein backtracking [START_REF] Armijo | Minimization of functions having Lipschitz continuous rst partial derivatives[END_REF] schemes forming x k+1 = x k + τ k d k . The scheme is repeated until a stopping criteria is reached. Firstly, choosing a suitable inexact line search maintains that the sequence of function values is monotonically decreasing (f k+1 ≤ f k ). Secondly, the sequence x k is converging globally. In other words, the method is convergent even if the initial point is far away from the minimizer. The rst case implies the minimization of the objective function and the second case provides that the method is not depending on the initial point. Particularly, Armijo`s line search satises

f (x k + τ k d k ) ≤ f (x k ) + c 1 τ k ∇f (x k )d k , (3.16)
where c 1 ∈ (0, 1 2 ) and τ k is the largest τ ∈ {s, ρs, . . .} with s > 0 and ρ ∈ (0, 1) such that (3.16) is satised. This concludes that function values satisfy the condition f (x k+1 ) ≤ f (x k ) imposing the monotonicity to the sequence of functions generated by this scheme and this scheme is globally convergent [START_REF] Armijo | Minimization of functions having Lipschitz continuous rst partial derivatives[END_REF][START_REF] Grippo | A nonmonotone line search technique for Newton's method[END_REF].

Methods to Solve Nonlinear Inverse Problem

However, in [START_REF] Grippo | A nonmonotone line search technique for Newton's method[END_REF] a modied version of Armijo condition is proposed such as

F χ 1 + S α χ k -τ K s χ k ≤ max k-M +1≤n≤k F (χ n ) -τ s S α χ k -τ K s χ k -χ k , (3.17)
where s is a small number and M is an integer. The right hand side of the new Armijo type line search is greater than the original Armijo`s rule implying that the new method can take bigger step sizes compared to the descent methods using original Armijo condition.

Thus, we can get faster convergence. This is also called weaker monotonicity. In original Armijo condition if no step size can be obtained to satisfy the condition the algorithm usually stops by rounding errors and preventing further progress. Further details can be found in [START_REF] Jin | Sparsity regularization for parameter identication problems[END_REF][START_REF] Grippo | A nonmonotone line search technique for Newton's method[END_REF]. The equation in (3.17) is motivated by the Taylor approximation of F χ k+1 as following

F χ k+1 = F χ k + F χ k + ξ χ k+1 -χ k χ k+1 -χ k , (3.18)
where ξ ∈ (0, 1). We can now restrict the gradient to a negative direction by choosing τ such that F χ k+1 ≤ F χ k . Hence,

F χ k+1 = F χ k -F χ k + ξ χ k+1 -χ k χ k+1 -χ k . (3.19)
Furthermore, we take a suciently small and xed s > 0 such as

sτ χ k+1 -χ k , χ k+1 -χ k ≤ F χ k + ξ χ k+1 -χ k χ k+1 -χ k . (3.20)
The aim is to get a negative direction, therefore we can estimate it as

F χ k+1 ≤ F χ k -sτ χ k+1 -χ k , χ k+1 -χ k . (3.21)
However, this estimation causes slow convergence because of strong monotonicity. Therefore, we change F χ k to max k-M +1≤n≤k F (χ n ).

Adding χ k+1 = S α χ kτ K s χ k leads us back to equation (3.17).

We should point out that we cannot use the rule developed by Brazilai and Borwein and the monotonicity criterion in the rst iteration of the algorithm. Therefore, we use an initial step size obtained by `fminsearch' which nds the minimum of a scalar function of several variables, starting at an initial estimate zero. This is generally referred to as order to obtain the required step size τ . We have also discovered that weak monotonicity does not improve the convergence rate of the cost function in our formulations compared with the back tracking line search with Armijo condition as can be seen in Fig. 3.

1
where iteration number which is 300 is our stopping criteria for the scatterer Fig. 3.4f. A standard monotonicity condition holds for M = 1. We have used the value M = 2 in our reconstructions.

Algorithm 1 The Barzilai-Borwein (BB) approach for choosing τ in line 6 of Algorithm 1.

1: Input The iteration counter k, τ 0 = 0,

h k = χ k -χ k-1 and g k = K (χ k ) -K (χ k-1 ) 2: if k = 0 then 3: τ = fminsearch(K (χ 1 ), τ 0 ) 4: else 5: output τ k = h k , h k / h k , g k 6: end for 3.2.

Soft Thresholding

Recently, algorithms known as iterative soft-thresholding (or forward-backward splitting) for optimization with sparse regularizers, including [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF][START_REF] Hale | A xed-point continuation method for l1regularized minimization with applications to compressed sensing[END_REF][START_REF] Nesterov | Gradient methods for minimizing composite functions[END_REF] have been used in many applications. These methods address problems of the form

min x f (x) = J 1 (x) + J 2 (x) (3.22)
wherein J 2 (x) is convex and possibly non-dierentiable while J 1 (x) is assumed to be differentiable and convex with a Lipschitz-continous gradient as stated before. These algorithms solve the non-smooth optimization problem directly with a projection-like operator instead of converting this problem into a constrained optimization problem. Particularly, these methods take the form of

x k+1 ← S J 2 (x k -α J 1 (x k ), α) , (3.23)
where S J 2 (x, α) is the solution of a `soft-threshold' problem at x with step size α and regularizer J 2 (x). Moreover, the soft-threshold operator is given by the solution to the soft-threshold problem

S J 2 (x, α) = argmin x 1 2 y -x 2 2 + αJ 2 (x). (3.24)
In our case, J 2 (x) = λ x 1 so the soft-threshold step for problem

min x f (x) = J 1 (x) + λ x 1 (3.25) would be argmin x 1 2 y -x k -α∇J 1 (x k ) 2 2 + α λ x 1 .
(3.26)

The soft-threshold rule is the nonlinear function dened as

(S α (t)) =          t -α, t > α 0, |t| ≤ α t + α t < -α (3.27)
or more compactly the shrinkage function dened as [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF][START_REF] Jin | Sparsity regularization for parameter identication problems[END_REF][START_REF] Gehre | Sparsity reconstruction in electrical impedance tomography: an experimental evaluation[END_REF][START_REF] Jin | A reconstruction algorithm for electrical impedance tomography based on sparsity regularization[END_REF] (S α (χ)

) i =    (|χ i | -α) sign(χ i ), if |χ i | > α 0,
otherwise.

(3.28)

It truncates small values to zero and shrinks large values as in Fig. 3.2. If we consider J 2 (x) = λ x 0 , S α (t) is termed hard thresholding function [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF]. However, we apply the sparsity constraint through the soft-thresholding function because of the reasons mentioned in the Introduction section recalling that the l 0 norm counts the number of non-zero elements in the solution, which is directly related to the solution's sparseness. However, this makes the minimization problem non-convex such that the problem might contain more than a single inmum [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Moreover, a global minimization of the cost function under l 0 -norm penalty is NP-hard problem (extensively burden computational cost is required to achieve the solution) [START_REF] Grasmair | Sparse regularization with lp penalty term[END_REF]. We rst take a step along the negative gradient of the function, and then compute this projection-like soft-threshold operator to take into account the eect of the regularizer.

The latter step eectively sparsies the result of the (generally dense) gradient step.

As discussed by [START_REF] Combettes | Signal recovery by proximal forward-backward splitting[END_REF], the soft-threshold operator is a generalization of the projection operator, and we recognize the iterative soft-thresholding algorithm as the classic gradientprojection algorithm but with projection replaced by soft-thresholding. Similar with the classic gradient projection algorithm, this algorithm may converge very slowly. However, analogous to the SPG algorithm which is a non-monotone projected gradient algorithm that combines the classical projected gradient method with the spectral gradient choice of step length and a non=monotone line-search strategy, [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF] proposes to use Barzilai-Borwein steps and a non-monotonic line search which we have discussed in previous section to speed up the convergence of the method.

It is important to underline that the regularization parameter determines the relative weight of the regularization term with respect to the other terms of the cost functional.

Therefore, if it is too small the sparsity regularization does not play a signicant role in the minimization process, whereas a large value entails that less importance is given to the tting of the data with respect to the sparsity enforcement. As we are dealing with a nonlinear problem, an optimal choice of this parameter is not straightforward and in any case depends on the available a-priori information.

Positivity Constraints

Imposing a priori constraints can improve the quality of solutions to the inverse problems in a great portion [54]. Non-negativity is important in applications like imaging [START_REF] Daubechies | An iterative thresholding algorithm for linear inverse problems with a sparsity constraint[END_REF]54].

We know that in order to have a physical solution with a convergence, there should not be a negative conductivity arising after a few iterations. However, the real and imaginary parts of the contrast function are not necessarily non-negative. On the other hand, there are possible ways which enforce the positivity constraints.

Projection Operator:

One possible way is to dene a closed set such as C = {χ ∈ Re n |χ ≥ 0} so that this condition can be imposed to our concerned problem. Thus, we can demonstrate the projection of χ onto a set C as

P C (χ) = arg min v∈C χ -v . (3.29)
In other words, P C (χ) is the closest point to χ in C [54]. We drop the subscript C from the projection operator for simplicity and dene the ith component of P (χ) as

[P (χ) i ] = max(χ i , 0) =    χ i if χ i ≥ 0 0 if χ i < 0. (3.30)

Reformulation of Parameters:

The other possibility can be to choose two new realvalued unknowns κ and η such that they can keep a-priori information. Let ε r (r) = 1 + κ 2 and σ(r) = η 2 with ε r (r) ≥ 1 and σ(r) ≥ 0. Thus, we can write the gradient of contrast function with respect to these unknowns such as dχ(r) dκ(r) = 2k 2 0 κ(r), and dχ(r) dη 2 (r) = 2k 2 0 η(r) ωε 0 (r) .

(3.31)

Hence, by applying a chain rule, the gradient of the cost functional is obtained as

dK(χ(η, κ)) dκ(r) = dK(χ(η, κ)) dχ(η, κ) × dχ(η, κ) dκ(r) .
(3.32)

In our simulations we have applied positivity constraint by projection rather than the reformulation of parameters. The reason for this is that positivity constraint by reformulation of parameters has a computationally burden and it slows down the convergence of cost function.

Stopping Criteria

The algorithm is terminated if the following holds

max k-ς+1 n k K(χ n+1 ) -K(χ n ) K(χ n+1 ) < s. (3.33)
This means the algorithm is terminated if the last ς steps have suciently small changes.

Another possible stopping criterion is to check τ when falling below a small positive constant. It can be interpreted as when maximum absolute value of step size times search direction goes below a small positive constant.

One other termination criteria can be the number of iterations. Having enough number of iterations can lead us to have a convergence in the cost function. In our simulations we have used 100 iterations unless it is stated otherwise.

Methods to Solve Nonlinear Inverse Problem

The algorithm which has been put forth, is as follows:

Algorithm 2 Steepest descent reconstruction algorithm (+ sparsity constraint)

1: Initialization χ 1 and α

2: for k = 1, • • • , T do 3: Solve the direct problem E di χ k 4: Compute the gradient K χ k = ∇ χ K (χ)| χ k
Adjoint method

5:

Smooth the gradient K s χ k by solving (3.9)

6:

Determine the step size τ j

7:

Update inhomogeneity by

χ k+1 = χ k -τ k K s χ k 8:
Threshold χ k+1 by S α χ k+1 (3.4) Sparsity constraint

9:

Imposing prior constraint via projection Positivity constraint 10:

check stopping criterion.

11: end for 12: output approximate the minimizer of (3.1)

We have used this algorithm throughout our simulations with addition of projection constraint when we call projection and the algorithm without projection constraint when we call pixel basis approach.

Soft Shrinkage Method in Wavelet Transform

The intrinsic multiresolution feature of the wavelet transform [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF] allows an accurate representation of the unknown function with a reduced number of coecients. In particular, the wavelet transform decomposes a given prole into two sets of coecients: coarse and detail. The coarse coecients account for the prole's low frequency content acting as a low-pass lter of the original function, while detail coecients account for high frequency content and allow representing the ner details of the function [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. As the decomposition level increases, the number of nonzero coecients is reduced. These properties motivate an approach that starts from a high level of decomposition, which is gradually reduced to retrieve ner details in the image. In practice, this is achieved by starting from a high order coarse representation and progressively moving to lower order coarse images to improve resolution.

In this section we will give a brief avor about the wavelet transform and reformulate the forward and the inverse problem of interest in terms of wavelet decomposition.

y(r) = k,l a J 0 ,k,l Φ J 0 ,k,l (r) + 3 s=1 J j=J 0 2 J -1 k,l=0 w s j,k,l Ψ s j,k,l (r), (3.34)
where r = (x, y) represents the spatial coordinates in the Cartesian system and J 0 is the initial decomposition level. Scaling functions Φ J 0 ,k,l and wavelet function Ψ s j,k,l are dened as tensor products of the scaled and translated scaling and wavelet functions in R such as

Φ j,k,l (r) = φ j,k (x)φ j,l (y), (3.35)
Ψ V j,k,l (r) = φ j,k (x)Ψ j,l (y), (3.36)
Ψ H j,k,l (r) = Ψ j,k (x)φ j,l (y), (3.37)
Ψ D j,k,l (r) = Ψ j,k (x)Ψ j,l (y), (3.38)
where Ψ j,k (x) is derived from the mother wavelet Ψ (x) as follows:

Ψ j,k (x) = 2 -j/2 Ψ (2 j x -k), (3.39)
and φ j,k (x) is the corresponding scaling function. This means that there is a unique scale function to compute the low frequency components in the previous decomposition level and three wavelet functions to compute the detail coecients along horizontal, vertical and diagonal directions (H,V,D).

The coecients a J 0 ,k,l and w s j,k,l are obtained by inner products in R 2 , a j,k,l = y, Φ j,k,l w j,k,l = y, Ψ j,k,l .

(3.40)

For the sake of simplicity, we rewrite the equation (3.34) as

y(r) = j c j C j (r), (3.41)
where c j and C j (r) are wavelet coecients and scaling basis functions for j = 1,

• • • , 2 2J 0 and for j = 2 2J 0 + 1, • • • , 2 2J
. We denote the calculations in (3.40) in an operator form following as c = Wy,

(3.42)
where W is a wavelet transform operator which maps y from R to a sparse set of coecients, c [START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF]. In our application, we are dealing with complex-valued functions. For instance, by using the linear property of wavelet transforms we obtain c = Wy re + iWy im = c re + ic im . LL, HL, LH and HH. The band HL indicates the variation of the original signal along the horizontal direction x and the subband LH for the vertical variation y.

Forward Problem in the Wavelet Domain

First, we explore the applicability of the wavelet representation for the contrast sources by solving the forward scattering problem. We start by rewriting the domain integral equation which is in the form of contrast source such as

J inc (r, r s ) = Ω [I -χ(r)G(r, r )]
A(r,r ) J(r , r s )dr .

J inc (r, r s ) = ∞ i=1 A i (r, Ω)J i (r s ) (3.45)
with A i (r, Ω) = Ω A(r, r )ψ i (r )dr , ψ i (r ) being wavelet basis and r s represents the position of sources. Applying basis decomposition along r on (3.45) leads to

J inc j (Ω, r s ) = ∞ i=1 A ij (Ω, Ω)J i (r s ), j = 1, • • • , ∞ (3.46)
with A ij (Ω, Ω) = Ω A i (r, Ω)ψ i (r)dr. Moreover, N pixels leads to N basis functions so that

J inc (r, r s ) = N i=1
A i (r, Ω)J i (r s .

(3.47)

Applying the basis decomposition along r on (3.47) leads to

J inc j (Ω, r s ) = N i=1 A ij (Ω, Ω)J i (r s ), j = 1, • • • , N (3.48)
with A ji (Ω, Ω) = Ω A i (r, Ω)ψ i (r)dr. Solving the state equation in the wavelet domain can be obtained as

      J 1 (r s ) J 2 (r s ) . . . J j=N (r s )       =       A 1,1 (Ω, Ω) A 1,2 (Ω, Ω) • • • A 1,N (Ω, Ω) A 2,1 (Ω, Ω) A 2,2 (Ω, Ω) • • • A 2,N (Ω, Ω) . . . . . . . . . . . . A N,1 (Ω, Ω) A N,2 (Ω, Ω) • • • A N,N (Ω, Ω)       -1       J inc 1 (Ω, r s ) J inc 2 (Ω, r s ) . . . J inc N (Ω, r s )       (3.49)
Combining (3.48) and (3.49) shows that the direct problem can be worked in the wavelet basis domain.

Inverse Algorithm in Wavelet Domain

Since we would like to represent the contrast function in the wavelet domain, we rst write the data equation as follows:

E di (r r , r s ) = Ω G(r r , r )W * (W(χ(r )E(r , r s )))dr , (3.50)
where W * is the inverse wavelet transform operator which maps χ back to a spatial domain. Similarly, the domain equation is written in terms of the contrast sources in the wavelet domain

W ((χ(r))E(r, r s )) = W χ(r)E inc (r) -W χ(r) D G(r, r )W * (W (χ(r )E(r , r s ))) .
(3.51) By using the above equations, we solve the nonlinear problem which minimizes the following cost functional:

F (χ) = ζ(χ) -G(r r , r )W * (W (χ(r )E(r , r s )) 2 K(χ) + W(χ) 1 . (3.52)
We update the unknown wavelet domain contrast function with the formula

W(χ k+1 ) = W(χ k ) -τ k W(K (χ k )), (3.53)
where

τ k = W(χ k ) -W(χ k-1 ), W(K χ k ) -W(K χ k-1 ) W(K (χ k )) -W(K (χ k-1 )), W(K (χ k )) -W(K (χ k-1 )) . (3.54)
In other words, to retrieve the unknown contrast function by solving (2.9), a wavelet expansion of the problem unknown in discretized form is considered as

E di (r r , r s ) = N i=1 c i D G or (r r , r )C i (r )E(r , r s )dr (3.55)
The algorithm used to solve inverse scattering problem is stated below:

Algorithm 3 Steepest descent reconstruction algorithm with sparsity constraint in wavelet domain

Initialize χ 1 and α

for k = 1, • • • , T do Solve the direct problem E di χ k Compute the gradient W(K χ k ) = ∇ χ K (χ)| χ k Determine the step size τ k Update inhomogeneity by W(χ k+1 ) = W(χ k ) -τ k W(K χ k ) Threshold W(χ k+1 ) by S α W(χ k+1 ) Back to spatial domain W * (χ k+1 )
check stopping criterion end for

Simulations

In this chapter, we present a number of numerical results to justify the eectiveness and accuracy of the proposed methods. We have tried dierent scatterers having dierent permittivities and conductivities in order to identify and compare the performance and efciency of each method with dierent constraints (i.e, addition of smoothness, projection, etc.). The description of physical characteristics scatterers and simulation conguration for each scatterer has been given in Tabs. 3.1 and 3.2 and in Fig. 3.4. The relative error norms as a comparison criterion are expressed as

ε err r = ε rec r -ε r 2 ε r 2 , σ err = σ rec -σ 2 σ 2 , (3.56)
where ε rec r and σ rec are the reconstructed permittivity and conductivity, respectively, and ε r and σ the exact ones. The minimum of the cost function is the one reached at the end of the process for each regularization parameter α. Additionally, the relative error norm on χ is dened as

χ err = Tr ( diag(χ) -diag(χ true ) 2 ) Tr ( diag(χ true ) 2 ) , (3.57)
where χ true is the true contrast to be found and Tr (A) stands for the trace of the matrix A. The discrepancy between the model and the measurement is as where E di (χ) is computed by solving the direct problem with the reconstructed χ. The laptop that has been used in order to run the simulations has a processor such as: Intel Core i7-4600U CPU@2.10 GHz×4

E di = N i i=1 ζ i -E di i (χ) 2 N i i=1 ζ i 2

Choice of Regularization Parameters

One of the key points of such an inversion either in spatial domain (pixel basis decomposition) or wavelet domain is the choice of the regularization parameters α in (3.4) and β in (3.9) (δ being kept constant and equal to 1). Dierent tests have been performed in order to evaluate the sensitivity of the choice of those parameters on the solution in spatial domain. The weighting parameters α and β should be chosen suitably so that we can drive the regularization properties of our algorithm in an ecient and predictable way.

In order to achieve a reconstruction which is close to the real case we make modications that aect the parameters. (g) Scatterer 7 From a practical point of view, supported by an extensive numerical analysis, it has been observed that a convenient choice for α can be 1/#of pixels. This choice is suggested by the fact that the penalty functions acts on the amplitude of the unknowns and it is so normalized by their number. On the other hand, this value is not the optimal one in our case. We have followed trial and error approach to be able to choose the most appropriate regularization parameter α in this application.

θ = [0, -π] θ = [0, -π]
At rst we consider the reconstruction as a function of α with a constant β as in Fig. 3.5 for the scatterer in Fig. 3.4f. Secondly, we examine the same reconstruction with a constant α as a function of β as in Fig. 3.6.

In the rst case we observe a better reconstruction of the exact case by increasing the parameter α such as in Figs. 3.8e and 3.8g or in Figs. 3.9e and 3.9g. The higher the value we choose for α the sharper the reconstructions are. However, the choice of α is not arbitrary as can be seen in Fig. 3.5. The error in permittivity gets larger for larger α. We get a minimum error when α is 2.5 × 10 -2 .

In the other case, we would expect the inclusion to get smoother when we increase β. However, in Fig. 3.6 it can be observed that there are no crucial dierences when β is 1, 1.5 and 2.5 dierently than 2. Fig. 3.7 shows that even if we do not have smoothness, small error can be obtained by a properly chosen regularization parameter. However, if we cannot choose the appropriate regularization parameter, addition of smoothness can help us to improve the quality of reconstruction. Nevertheless, we should choose β carefully and close to 1 in order not to loose the dierentiability for the next iteration. Using β as 1.5 shows that the reconstruction becomes smoother as can be seen in Figs. 3.8i, Figs. 3.8j, 3.9i, and 3.9j compared to the gures without any smoothness.

Addition of Projection

Addition of projection to the sparsity reconstruction method improves the quality of the reconstruction even more and the error range for both permittivity and conductivity decreases until a value of α such as 1.5 × 10 -2 for permittivity as in Figs. 3.10b and 3.11b in spatial domain. In here, dierent randomization of 10 dB noise levels have been tested to see how the addition of noise data can be discretized. Figs. 3.10b, 3.10d, 3.11b, and 3.11d show us that we can choose the regularization parameter in a wider range when we add projection to our algorithm. On the other hand, we cannot base our solutions on the minimum cost function as it does not give us the best solution even if we add projection as in Fig. 3.14. This shows us that our problem is ill-posed and we need a-priori information. Figs. 3.12 and 3.13 present a comparison of the maps of the permittivity and conductivity, respectively for various inversion parameters where Tab.3.3 exhibits the error in permittivity in each case. The inuence of the α parameter without any projection constraint onto the permittivity and conductivity is important as we get better reconstruction when α is 5 × 10 -3 than when α is 1 × 10 -6 .

On the other hand, a clearer image is obtained by addition of projection and choice of α has wider range. Comparing Figs. 3.12 and 3.13 shows us that when the permittivity is low, the reconstruction quality is better. Moreover, Tab.3.3 highlights that we can get the minimum error no matter the regularization parameter is, if we have addition of projection. Even if we have less sparsity enforcement (smaller regularization parameter), addition of projection to our inverse problem yields a better quality of reconstruction.

Analysis of Proposed Approaches

In this section, the eectiveness of the proposed approaches has been assessed and compared against dierent scatterers as in Fig. The SWT is a wavelet transform achieving translation invariance, a property that is missing in the DWT [START_REF] Candes | Compressed sensing with coherent and redundant dictionaries[END_REF]. Translation-invariance is achieved by removing the downsamplers and upsamplers in the DWT and by upsampling the lter coecients by a factor of 2j at the (j -1) level of the algorithm. For this reason, the output of each level of the SWT contains the same number of samples as the input. In comparison with the previous case of step functions, this means that the unknown is now expressed by means of a redundant dictionary in which there are more columns than rows [START_REF] Candes | Compressed sensing with coherent and redundant dictionaries[END_REF]. In other words, such a representation is not univocal, as the considered functions are not orthonormal.

The use of these overcomplete dictionaries is now widespread in signal processing and data analysis, as there are numerous practical examples in which a signal is not sparse in an orthonormal basis or incoherent dictionary, but it is instead sparse in terms of a truly redundant dictionary [START_REF] Candes | Compressed sensing with coherent and redundant dictionaries[END_REF].

On the other hand, `Wavelab' which is a library of Matlab routines for wavelet analysis, wavelet-packet analysis, cosinepacket analysis and matching pursuit is used in Matlab to study our problem in Discrete Wavelet Transform. DWT is any wavelet transform where the wavelets are discretely sampled. In practice, we only have a nite number of values of the contrast function in a rectangular domain D, which is discretized into 2 J × 2 J grid cells. The values of the contrast function on these cells represent the scaling coecients at the highest scale J (decomposition level). In our case, we have used 32 × 32 (J = 5)

discretization for each geometry while solving our inverse problem in the wavelet domain.

As mentioned before the choice of the wavelet family plays an important role since the priority is to choose a wavelet basis that produces as many zero coecients as possible. This choice is based on the analysis of three main criteria which are the vanishing moments, support size and the regularity. High-amplitude coecients occur when the support of the wavelet overlap with transitions such as edge and they are proportional to the width of the wavelet support, which needs preferably to be as small as possible.

However, the number of vanishing moments is proportional to the support size so that the choice of the optimal wavelet is a trade-o between the number of vanishing moment and support size. Furthermore, the wavelet regularity can reduce the visibility of artifacts.

Considering all these features and the Tab.3.4 it can be observed that the Daubechies and Symmlet families represent good options to be employed in the numerical analysis.

However, as it was pointed out before, support size with respect to our unknown scatterers plays an important role. Taking account the sparsity of our unknown scatterers within these wavelet basis, we will rstly exploit the `Haar' (db1) wavelet and do the comparison with Daubechies (db4) wavelet basis.

Secondly, we will do the optimization on detail coecients of the contrast function, which are its ner details by exploiting Haar basis for SWT (swt-db1-opt.details) and

Wavelab case (WL-db1-opt.details). Before passing to the comparisons of the methods with dierent characteristic properties of wavelet functions, we would like to provide some notations for simplicity:

• Haar wavelet: db1

• Daubechies 4 wavelet: db4

• Optimization on both detail and approximation coecients by wavelet decomposition using Wavelab: WL

• Optimization on both detail and approximation coecients by wavelet decomposition using Stationary Wavelet Transform: SWT

• Optimization on detail coecients by exploiting Haar basis through SWT: swt-db1opt.details

• Optimization on detail coecients by exploiting Haar basis through Wavelab: WL-db1-opt.details

• level two: l = 2 and level three: l = 3.

We have used 20 dB and 10 dB noise data set in order to illustrate our proposed approaches. Moreover, we exclude the smoothing the gradient constraint in soft-shrinkage algorithm both in spatial and wavelet domains for the rest of the simulations since addition of smoothness does not improve the solution accuracy with the proper choice of regularization parameter.

20 dB noise data: For instance, the plots of the real part and the imaginary one in Fig. 3.17 of the contrast retrieved from blurred scattering data with a SNR on the scattered eld amplitude equal to 20 dB prove a good reconstruction of the actual contrast despite its nonsparse nature with respect to a standard single-pixel basis. Comparing the performance Property Haar

Daubechies N Symmlet N Coiet N Bi-orthogonal Bi-orthogonal Reverse N r .N d N r .N d Support width 1 2N -1 2N -1 6N -1 2N r + 1 2N d + 1 Number of 1 N N 2N N r N d vanishing moments Regularity ≈ 0.2 ≈ 0.2N 2N Not dened N r -1 N d -1
Table 3.4.: Features of the main wavelet families of the approach in the wavelet domain with the one in pixel basis, wavelet based soft shrinkage reconstruction scheme is faster. However, as it can be noticed in Tab. 3.5, the wavelet based approach improvement is of about 0.3% which is quite low. Moreover, addition of projection to the pixel based reconstruction algorithm gives us better result.

On the other hand, choosing Daubechies (db4) basis decomposition does not give us a good quality of reconstruction while its performance is slower when compared to others.

For the case of Scatterer 9 (Fig. 3.4i), we can see from Tab. 3.5 when we increase the decomposition level from level 2 to 3 we get a slightly better solution. Even if there is not a big dierence in the reconstruction error of permittivity between Haar and Daubechies (db4) basis decompositions, the reconstruction obtained by Wavelab with Haar basis is smoother than the one obtained with db4.

10 dB noise data: When we increase the noise such as in Tab. 3.6, the reconstruction error increases directly. However, interestingly, the error dierence between WL-db1 and pixel reconstruction schemes is around 0.0015, which is not remarkable while the scatterer is not exactly sparse in the domain of the pixel basis function. On the contrary, the performance of WL-db1 is better than the one with pixel basis decomposition.

When we consider scatterer 7 (Fig. 3.4g), which counts as sparse with respect to the pixel basis, we can see in Tab. 3.5, 3.6 that the soft shrinkage algorithm gives a slightly better reconstruction in the spatial domain. On the other hand, it is worthwhile to mention that the algorithm in wavelet decomposition also achieves a good resolution of the image whereas wavelet decomposition in terms of swt is not preferable because of its low performance and lower reconstruction quality.

We also consider the applicability of the proposed methods with the scatterers which have partial views (Fig. 3.4d, 3.4j) and a case where the sources and receivers are located on the top half of the region of interest as in Fig. 3.4h and Fig. 3.4k. Compared to the full view scatterers, these congurations have worse reconstruction quality. However, scatterer can be seen in Tab. 3.7 and Tab. 3.8 especially when we have the wavelet decomposition with swt. For the scatterer 11 (Fig. 3.4k) we have better results by using wavelet basis functions than the pixel functions albeit the scatterer is sparse with respect to pixel basis.

When we have both reection and transmission congurations the eectiveness of the proposed approaches reduces as the a-priori information is limited due to the information obtained at the receivers. However, we can still obtain the best approximation to the actual proles by enforcing sparsity through projection as it can be seen in Tabs. 3.7 and 3.8. Like for all other results, application of projection gives us the better minimum error for the permittivity and conductivity for both 20 dB and 10 dB noise data sets while keeping its slow performance compared to performances of other approaches.

It is also worthwhile to mention that in [START_REF] Ambrosanio | Exploiting wavelet decomposition to enhance sparse recovery in microwave imaging[END_REF] the authors studied wavelet decomposition in breast cancer imaging and according to their results the wavelet decomposition recovers the unknown prole around 15% more than the pixel basis case based on the fact that there is no frequency hopping for wavelet decomposition. However, in studies [START_REF] Li | A contrast source inversion method in the wavelet domain[END_REF][START_REF] Bevacqua | Non-linear inverse scattering via sparsity regularized contrast source inversion[END_REF][START_REF] Anselmi | Wavelet-based compressive imaging of sparse targets[END_REF] wavelet decomposition gives better results than pixel basis decomposition taking into consideration scatterers which are sparse in wavelet domain and not in spatial one. We extend our proposed inversion algorithm based on soft shrinkage enforcement over the reconstruction of small isolated scatterers. By adopting wavelet basis functions we discretized the formulation of the 2D inverse scattering problem where the imaging problem has been solved with a soft shrinkage reconstruction algorithm as a second technique.

Selected results from numerical experiments have been presented and discussed to give some insights about the robustness, the exibility, and the accuracy of the proposed approaches as well as to illustrate their advantages and limitations with respect to dierently constrained inversion methods.

Since sparsity is a relative concept, the eectiveness of each approach depends on the actual scatterer and the adopted representation basis, i.e, on the available a-priori information on the scattering scenario to be used for choosing the basis. We have studied dierent scatterers within this respect. The results showed us that whenever the proper wavelet basis is selected for the scenario at hand (Fig. 3.4), the proposed method based on soft shrinkage enforcement gives good results in terms of computational eciency and accuracy (Tabs. 3.5, 3.6, 3.7, 3.8). The wavelet based method can be reliably applied to a wide set of scattering congurations by handling a reduced set of scattering data.

Overall, even if the soft shrinkage algorithm which has been studied both in spatial domain and wavelet domain gives us satisfactory results, an inversion based on projection constraint gives us slightly better results yet its performance is more time consuming.

Two-Step Inversion Method

Motivation

Additional to the interest of exploiting sparsity which has been mentioned in previous chapters, the other key point of the sparse estimation problem is to reveal the identication of the support, which denotes the indices of the nonzeros. If the support is known, the estimation of the sparse vectors reduces to a standard overdetermined linear inverse problem [START_REF] Lee | Subspace methods for joint sparse recovery[END_REF].

A wide variety of the studies starts to contemplate the structure information of the solutions in order to facilitate a better estimation. For instance, `group sparsity' structure [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF] has been presented in many applications, where the components of solutions are likely to be either all zero or all nonzero in a group. Thus, one aims to decrease the dispersion to enhance the solution by taking account the grouping prior. In this chapter, we focus on joint sparsity, which is a special case of the group sparsity. Specically, joint sparsity means that multiple unknown sparse vectors (x j ∈ R n , j = 1, • • • , s) share a common unknown nonzero support set [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF][START_REF] Fan | Enhanced joint sparsity via iterative support detection[END_REF].

We propose a two-step inversion approach in order to solve a nonlinear inverse problem as ours by applying joint sparsity to get the equivalent current, then the unknown contrast.

We show the eciency, accuracy and the limitations of our proposed method with the results obtained in 2D and 3D.

Two-Step Inversion

In compressive sensing, joint sparsity desires to reconstruct unknown signals from m measurement vectors based on a common measurement matrix. This is also called the multiple measurement vectors (MMV) problem [START_REF] Sarvotham | Distributed compressed sensing of jointly sparse signals[END_REF]. Given the vectors y j , ∈ R m , and a measurement matrix A ∈ R m×n , we want to recover the x j from the noisy underdetermined systems y j = Ax j + n j (j = 1, • • • , s), where n j is the noise vector.

= Ax j + n j , j = 1, • • • , s, (4.1)
where |M | is the cardinality of M [START_REF] Heckel | Joint sparsity with dierent measurement matrices[END_REF]. Since (4.1) is NP-hard [START_REF] Natarajan | Sparse approximate solutions to linear systems[END_REF][START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF], this problem is usually relaxed with a convex alternative which is computationally ecient at the expense of more being required measurements. Like l 1 -norm being the convex relaxation of l 0 -norm [START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF], the (weighted) l 2,1 -norm is widely used as the convex replacement of |M | as stated below [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF][START_REF] Heckel | Joint sparsity with dierent measurement matrices[END_REF]:

min X X w,2,1 := n i=1 w i x i 2 s.t AX = B + N, (4.2)
where A ∈ R m×n and B ∈ R m×l is an available measurement matrix consisting of l measurement vectors, and x j ∈ R n denote the i-th row and the j-th column of X, whereas X = [x 1 , • • • , x s ] ∈ R n×l denotes a collection of l jointly sparse solutions (unknown source matrix) while w i ≥ 0. N ∈ R m×l is an unknown noise matrix.

A key assumption in the MMV model is that the support of every column of X is identical. Similarly with the constraint in the single measurement vector model, the number of nonzero rows in X has to be below a threshold to maintain a unique global solution [START_REF] Cotter | Sparse solutions to linear inverse problems with multiple measurement vectors[END_REF]. This leads to the fact that X has a small number of nonzero rows.

Following all the above properties, the other study of the thesis is coming from adopting a two-step method which rstly consists of nding the equivalent current J i for i = 1, . . . , N s (N s represents the number of sources) using (2.9) and then looking for χ by combining (2.14) and (2.11) instead of solving directly the nonlinear inverse problem given by (3.2) as we studied in Chapter 3. This approach leads us to have two linear minimization problems solved by exploiting the jointly-sparse aspect of the sought equivalent currents solution of the rst step and a classical l 2 -minimization of a linear problem for the second step.

First step: Reconstruction of the equivalent currents

The following optimization problem is solved

J i = argmin J i 1 2 ζ i -G or J i 2 i = 1, • • • , N s . (4.3)
where ζ i is a vector of size N r which gathered the signal due to the source #i measured by the N r receivers.

The main idea of the approach is to take into account that, as shown by (2.24), J i and χ share the same support which means that when χ j = 0 then J j i = 0, ∀i = 1, . . . , N s where χ j and J j i are the j th element of χ and J i respectively. Taking into account this hypothesis,

G or J i = ζ i (4.4)
while w j is the corresponding weight [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF].

In [START_REF] Yang | Alternating direction algorithms for l 1 -problems in compressive sensing[END_REF], Yong et al. to the group version for solving the group sparse optimization with l 2,1 -norm regularization in [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF]. In this respect, we use YALL1 Group [START_REF] Deng | Group sparse optimization by alternating direction method[END_REF] package which encodes the joint sparsity model in order to solve (4.4).

Alternatively, we use T-MSBL which is a block sparse Bayesian learning to solve (4.4).

T-MBSL identies the MMV model in order to exploit the correlation that exists in each nonzero row of X while automatically choosing the optimal regularization value.

More details related to this algorithm can be found in [START_REF] Zhang | Sparse signal recovery with temporally correlated source vectors using sparse Bayesian learning[END_REF]. The MATLAB code for both algorithms can be obtained online as well.

In [START_REF] Yang | Alternating direction algorithms for l 1 -problems in compressive sensing[END_REF], it has been also proved that YALL1 gives better performance and solution accuracy compared to state of art algorithms such as spectral projection gradient method (SPGL1) [START_REF] Chen | Atomic decomposition by basis pursuit[END_REF], SpaRSA (a sparse reconstruction algorithm for more general regularizers) [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF], FISTA (a fast iterative shrinkage thresholding algorithm that attains an optimal convergence rate in function values) [START_REF] Beck | A fast iterative shrinkage-thresholding algorithm for linear inverse problems[END_REF]. Therefore, we have used YALL1 toolbox and T-MSBL which applies joint sparsity approach by adopting a probabilistic approach to incorporate correlation structure in each nonzero row of the solution matrix dierently than the existing algorithms to apply our proposed method.

Second step: Reconstruction of the contrast function

Once J i is known, the contrast function is obtained by solving the following minimization problem [START_REF] Abubakar | Iterative forward and inverse algorithms based on domain integral equations for three-dimensional electric and magnetic objects[END_REF] χ = argmin

Ns i=1 J i -diag (χ)E i 2 = Ns i=1 J i • Ēi Ns i=1 E i • Ēi , (4.5)
where E i has been obtained using (2.25) and Ēi being its conjugate.

Two-Step Inversion Method Results in 2D

It is worthwhile to restate that sparseness is a relative concept with respect to a basis [4].

In the case of a pixel basis a fast way to estimate the sparsity of our problem is to dene it as the ratio of the obstacle's areas to the investigation domain D area. By this way, sparseness of the rst example (Fig. 3.4b) is around 3 % whereas the one of the second example (Fig. 3.4f) is around 11 %. One of the key points of our examples is that in the rst one the scatterers are sparse with respect to their pixel expansion basis while in the second one the scatterer is less sparse and our proposed approach in this chapter performs dierently for these two cases.

We exploit the joint sparsity by using YALL1 and T-MSBL algorithms in our twostep inversion method. Thereafter, we compare these approaches with SpaRSA and Born Approximation (B.A). SpaRSA is a method based on iterative shrinkage/thresholding (IST) [START_REF] Wright | Sparse reconstruction by separable approximation[END_REF]. It computes the steepest descent direction on the l 2 norm in (4.3) and uses a very simple soft thresholding function related to the regularization term.

As for the parameters of YALL1 algorithm we choose the primal-based solver where the linear system is exactly solved. The weights have been chosen as w i = 1 (i = 1, • • • , N ) and the initialization J ≡ 0 (no prior information). We have also not added projection while keeping the number of iterations as 10 × # of measurements.

It can be observed that when we have sparse scatterers (as in Fig. 4.1) the two-step inversion approach gives a better reconstruction in terms of quality of χ and the processing time compared to the soft shrinkage algorithm as shown in Tab. 4.1. However, application of joint sparsity through T-MSBL is more favorable than the YALL1 when we have sparse scatterers with respect to their basis function. Error on contrast and permittivity is calculated as in equations 3.56 and 3.57

On the other hand, even if the two-step inversion approach is faster in the reconstruction of the scatterer in Fig. 4.2, the soft shrinkage algorithm yields a better reconstruction quality with an appropriate regularization parameter. Furthermore, one can also observe that when we decrease the number of discretization (comparing Fig. 3.4b to Fig. 3.4g), the dierence in reconstruction quality is quite small such as 1.21%. Also, Born approximation and SpaRSA approaches give us poor reconstruction even if they are less time consuming as it can be seen in Tab. 4.1.

When we transform the two-step inversion approach from spatial to wavelet domain, our minimization problem transforms into

min J J w,2,1 := N j=1 w j Jj 2 s. t. W(G or J i ) = W(ζ i ), (4.6)
where J = N j=1 J j C j implying J = W(J) as in (3.42). {C j } N j=1 represent the wavelet functions and {J j } are the decomposition coecients.

When we compare both scatterers (Fig. 3.4f and Fig. 3.4g) in the wavelet domain (WD)

and the spatial domain, soft shrinkage algorithm provides approximately the same error.

However, the soft shrinkage algorithm in wavelet domain achieves slightly better results than the two-step inversion method applied through YALL1 and T-MSBL as shown in Tab. 4.2. On the other hand, when we compare the approaches studied in the spatial and the wavelet domain, we can see in Tab. 4.1 and Tab. 4.2 for the case of Fig. 3.4g where we have 5 small scatterers, the two-step inversion method through T-MSBL gives better results in the spatial domain than is the wavelet domain.

Furthermore, as it can be seen in Fig. 4.3 when we reduce the noise level on our synthetic data, the error on the contrast function shows us that the eectiveness of the approximation gets better. For instance, when the scatterer is as in Fig. 3.4f, the error obtained through the soft shrinkage algorithm gives us better estimate than the other methods applied.

Contrarily, the two-step inversion method applied by T-MSBL in the spatial domain reduces the error more than the other methods are doing when the scatterer is small as in Fig. 3.4g. On the other hand, we can observe that the soft shrinkage algorithm provides a good approximation both in the spatial and the wavelet domains even if we have our data disturbed with a high level of noise. Let us consider a complex, multi-layer structure, where each of the n layers is assumed to be non-magnetic (µ = µ 0 ) homogeneous uniaxial. (This type of structure usually results from a large-scale (homogenized) view of bered planar laminates as in aeronautics.) In the local coordinate system (material frame), each layer is characterized by a diagonal complex permittivity tensor

¯ (n) e =    (n) 11 0 0 0 (n) 22 0 0 0 (n) 22    Ξ (θ n ) =    cos θ n sin θ n 0 -sin θ n cos θ n 0 0 0 1    . (4.7)
One is able to carry the local coordinate system to the global Cartesian one where θ n is the rotation Euler angle. The electromagnetic response of such a structure which satises the radiation condition at innity (the Sommerfeld condition) and in a complex multi-layer continuity conditions at each interface is computed as

E di j (r) = iωµ 0 V Ḡee (r, r )• χ (r )•E tot j (r ) dr , with χ (r) = -iω 0 Ξ-1 (θ n )• ¯ i -¯ (n) e • Ξ (θ n ) (4.8) χ(r) • E inc j (r) = J j (r) -iω 0 -χ(r) • iωµ 0 V Ḡee (r, r ) • J j (r ) d r (4.9)
where Ḡee (r, r ) is the electric-electric dyadic Green's function, ¯ i the background per- mittivity tensor, ¯ e the permittivity tensor of an inclusion of volume V within the back- groundand volume V . E di j (r), E tot j (r) and E inc j (r) are the scattered, total and incident elds respectively due to the j th source. For theoretical and numerical details one should refer to [START_REF] Richmond | Scattering by a dielectric cylinder of arbitrary cross section shape[END_REF], [START_REF] Zhong | Electromagnetic response of anisotropic laminates to distributed sources[END_REF] and [START_REF] Zhong | Fast calculation of scattering by 3-d inhomogeneities in uniaxial anisotropic multilayers[END_REF]. Here is the description of the conguration on which the two-step inversion is applied: • Fiber-glass composite ¯ 1 = 0 diag [5.46 + i2.29, 5.21 + i2.08, 5.21 + i2.08]

• Frequency f = 6 GHz; λ 0 = 5 cm; λ 1 = 2.
θ 1 = 0 • • Two inclusions of complex permittivity Ω = (6, 0), ¯ inclusion = Ī 0 Ω Size (0.2 × 0.2 × 0.2) λ 1 , centered at (1.4, 0.7, 0.7) × λ 1 Size (0.2 × 0.2 × 0.2) λ 1 , centered at (0.5, 1.4, 1.4) × λ 1 • Region Of Interest: Size: l x × l y × l z with l x = l y = 2λ 1 and l z = 2λ 1 Discretization: n x × n y × n z , n x = n y = n z = 10
Depth of top of ROI z = 0.25λ 1

In this section, the unknown scatterer is not too large and/or its contrast not too high.

The preliminary results obtained using two-step inversion method already show some good results. In [START_REF] Oliveri | 3-D crack detection in anisotropic layered media through a sparseness-regularized solver[END_REF] a Bayesian Compressive Sensing solver, is used whereas, in this work, an algorithm exploiting a joint sparsity regularization is proposed. Let us notice that results using the one-shot MUSIC retrieval method and the iterative Subspace Optimization Method in the 3D anisotropic case are put in perspective in [START_REF] Lesselier | On inverse scattering and imaging solutions for objects buried within uniaxially anisotropic media[END_REF], this being behind the scope of our discussion here however.

From Fig. 4.6 we can see that the two-step inversion method applied through YALL1

gives us a good localization and estimation of the contrast when there is no noise and it maintains a good convergence as shown in Fig. 4.5. However, the method does not give us good results when we have noise in our data as in Fig. 4.7.

Similarly, the two-step inversion method applied through T-MSBL gives us a good reconstruction for the real and imaginary parts of the contrast (see Fig. when there is no noise. However, it can still keep its good reconstruction quality with the addition of the noise to the data as can be seen Fig. 4.10 and Fig. 4.11. We have used 3000 iterations while running our simulations. On the other hand, this method is sensitive to initialization.

In 3D, nonlinearity and ill-posedness become even more severe than 2D applications.

This reminds us of the importance of the enforcement of the sparsity constraint. As it can be observed, the two-step inversion gives better results when it is applied through T-MSBL than the ones of YALL1. This might be due to the stronger correlation that applied through T-MSBL to the coecients of each row of the unknown contrast. This point as well as other open questions will be considered in the Conclusion chapter.

Conclusion

In this chapter, we presented a two-step inverse process which allows sparse recovery of the unknown (complex) dielectric proles of scatterers. The proposed approach is correlated with joint sparsity which gives multiple sparse solutions that share a common nonzero support. The principal interest of such a method is that it accurately reconstructs the unknown scatterers without linear approximation and presents an ecient recovery algorithm of sparse scatterers by reducing possible ambiguities on the scatterer sparsity deduced from the null values of the equivalent current [4,[START_REF] Candès | An introduction to compressive sampling[END_REF].

Thereafter, we compared the method proposed in this chapter with the approach that directly obtains the contrast through enforcement of sparsity by soft shrinkage thresholding. Both approaches produce sharp and good reconstruction of dielectric proles in sparse domains and keep their convergence during the reconstruction.

On the other hand, when we have a scatterer which is sparser with respect to pixel basis then two-step inversion method through T-MSBL gives better results in the spatial domain than the wavelet domain. This might be due to the correlation of the nonzero rows of the unknown contrast that we are interested in. When we have a scatterer which is not that sparse in spatial domain then there is not a big dierence for both two-step inversion and soft-shrinkage approaches either in spatial or wavelet domains. Moreover, a three-dimensional inversion strategy is introduced for the detection of scatterer in uniaxially anisotropic layered media having principal axes with arbitrary orientation. The two-step method is adopted for the solution of the imaging problem. A set of preliminary numerical results is reported to assess the accuracy of the proposed method even if more studies have to be done in order to evaluate the limitations of the proposed method.

Recently, several computational advances have been made in the nonconvex sparse regularization since its performance is better than the one of the convex sparse regularization.

While there exist many algorithms for solving the nonconvex sparse regularized models, it is still a challenging problem to obtain the global optimal solution eciently. In addition, the behavior of a local solution is hard to analyze and, more importantly, structural information of the solution is also hard to be incorporated into these algorithms [START_REF] Fan | Enhanced joint sparsity via iterative support detection[END_REF].

Summary, Conclusions and Future

Work

In this last chapter, the work that has been done and presented in the previous chapters is summarized. The main points, contributions, results and inferences are summarized to provide a global view and highlight the goal of the project. The conclusions are subsequently followed by a list of suggestions and ideas for future work. These perspectives aspire to complement the work carried out and presented here and, if possible, to bring it closer to real life applications.

Summary and Conclusions

In this dissertation, we have considered the non-linear microwave imaging problem by transforming it into a problem of sparse signal representation using dictionaries. This is a very attractive way of looking at retrieving proles of unknown scatterers in a region of interest because when the scatterers can be well-modeled with an appropriate basis, and the number of their unknown coecients is small, then the true underlying spatial spectrum is sparse. The problem of signal representation in bases is an ill-posed nonlinear inverse problem, and as such, it requires regularization to have unique well-behaved solutions. We are interested in sparse signal representations, so the regularization has to enforce sparsity.

To enforce sparsity we utilized l p penalties with p ≤ 1. There is an important distinction between l 1 penalties (p = 1), and penalties with p < 1. For the l 1 case, the penalty leads to convex optimization problems, whereas for p < 1, the associated optimization is nonconvex. For the sake of simplicity in terms of computational complexity we focused on optimization involving l 1 penalties for which we used an algorithm that is a special case of gradient descent. The algorithm has the important benet of allowing ecient global solutions by adapting dierent properties.

In this work, thresholded nonlinear Landweber iterations (special case of gradient descent) are used to solve the sparse minimization problem constructed directly from the nonlinear scattering equations. It is adapted from the classical iterative soft shrinkage algorithm, and its main ingredients include a Sobolev smoothing of the estimated gradients, a soft shrinkage iteration, and an adaptive step size selection based on the Barzilai and Borwein rule. The resulting scheme requires the user to set only one simulation parameter before the execution, which signicantly simplies the application of the method to different problems. Additionally, since the regularization is applied directly to the nonlinear problem without any linearization approximation, the scheme can be applied to higher contrast level of scatterers. Furthermore, a projected steepest descent algorithm, which increases convergence rate of the nonlinear Landweber iterations, is used for reducing potential high computational cost. The projection operator replaces the thresholding function and enforces the sparsity constraint. Indeed, numerical results demonstrate that the resulting projection scheme is very ecient in recovering permittivity proles. We adopted as already said the Barzilai-Borwein direction as a step size selection, and the numerical results show the signicant development in eciency of the gradient method. Being computationally ecient and needing low memory requirement makes this scheme interesting to solve large-scale optimization problems.

On the other hand, the decoding process requires nding a sparse solution of an underdetermined linear system. What makes such a scheme work is sparsity; i.e., the original signal must have a sparse or nearly sparse representation under some known basis. Following this fact, we used the wavelet transform to retrieve the unknown scatterer prole by solving the sparse Tikhonov minimization problem constructed directly from the nonlinear scattering equations. This inversion is usually more robust than simultaneously inverting all the pixels in the spatial domain because the number of unknowns is less and the non-uniqueness is thereby reduced. We explained how to extract the dierence of information between successive resolutions and thus dene a new (complete) representation called the wavelet representation. This representation is computed by decomposing the original signal using a wavelet orthonormal basis.

Let us remind that the scatterers are sparse with respect to their expansion basis.

In this work, we have studied dierent proles of scatterers which are sparse in pixel basis and/or sparse in Haar or in Daubechies basis. The integration of wavelet bases within sparseness-regularized formulations for microwave imaging has been done. Then, the introduction of a generalized regularized imaging strategy applicable to a very wide set of scattering scenarios depending on the chosen wavelet family (i.e., the available apriori information about the scattering scenario) has been worked out while keeping the advantages of pixel-basis optimization techniques.

The wavelet transform is localized in space, and we can rst invert only the lower level wavelet coecients to get a rough image. This inversion is usually more robust than simultaneously inverting all the pixels in the spatial domain because the number of unknowns is less and the non-uniqueness is thereby reduced. After this, we can improve the inversion results by adding wavelet coecients at higher levels. On the other hand, when we increase the decomposition level in order to have the optimization on sparser coecients, there is no such large dierence both in accuracy and in performance for the examples that we have studied.

The analysis in Chapter 3 shows that a soft shrinkage algorithm is suciently robust to work even when the a-priori knowledge about the scatterer is only approximate, that is, the scatterer is not exactly sparse within the domain of the considered basis functions. Moreover, the reconstruction accuracy turns out to be acceptable also when higher contrast and lower SNRs are.

It has been analysed that the discrete wavelet transform shows a better performance and eectiveness in the reconstructions compared to the stationary wavelet transform in some cases. However, the error on retrieving various proles is quite small between each transformation. One can also conclude that the accuracy, the convergence rate, and the robustness vs. noise of the proposed approaches are not dependent only on the problem setup, but also on the choice of the wavelet basis for the contrast expansion, which is of course strongly related to the available a-priori information.

The third direction of the work carried out in the thesis is exploiting the joint-sparsity benets by applying a two-step method to solve our non-linear problem of concern. The two-step inversion method uses the contrast-source formulation of the scattering equations, which allows for straightforward application of the sparsity constraint to the solution. The nonlinear system of equations is solved using an iterative scheme, which calls for the solution of a linear system at every iteration.

The two-step inversion method is applied through available YALL1 and T-MSBL algorithms both in spatial and wavelet domains. Both methods exploit joint-sparsity which reduces possible ambiguities on the scatterer sparsity deduced from the null values of the equivalent current. This approach consumes less time and ensures better imaging quality compared to an iterative method with soft thresholding in the reconstruction of sparse scatterers. However, even if the two-step method is faster than the usual soft shrinkage algorithm, the reconstruction of the scatterers is better with soft shrinkage algorithm in wavelet domain. Above of all this, addition of projection to our proposed methods always gives us better approximation for the unknown.

Future Work

A very important issue in the framework of this work is the choice of the regularization parameters, α in the l 1 formulations and β and δ in smooth gradient formulations. It is worthwhile to continue the investigation of methods for regularization parameter selection from other elds since being able to identify the most accurate parameter is challenging.

In our application we have not be able to use `L-curve'. There are other methods such cross-validation, and universal and min-max regularization parameter selection rules as well. The viability of these methods for our problem is on discussion for our specic case of application. However, much more work has to be made to get insights into how to select the regularization parameter for our problem, or to dismiss these methods as inappropriate for our problem.

In the thesis we exclusively used l 1 regularization for enforcing sparsity. However, regularization that favors sparsity is not limited to this. Many other forms exist, such as Huber regularization, entropy-based regularization [START_REF] Karl | Regularization in image restoration and reconstruction[END_REF], and the other ones stated in Introduction section. Also, an analysis of the specic features that are necessary for the regularizing term to favor sparsity would provide much insight into the selection of a particular functional. Such analysis has been previously done at some level [START_REF] Donoho | Maximum entropy and the nearly black object[END_REF][START_REF] Charbonnier | Deterministic edgepreserving regularization in computed imaging[END_REF], but deeper understanding can be gained by putting the analysis on rm theoretical grounds and considering much wider sets of regularizing functionals.

All of the schemes proposed in this work, namely soft shrinkage algorithm in spatial and wavelet domain, and two-step inversion method, can be extended to imaging of domains residing in layered media. This extension calls for replacing the free-space Green function with the planary or cylindrically layered medium Green functions. The resulting inversion tools have applications in subsurface and borehole imaging, as an example.

The eciency of the sparsity-constrained regularization can be increased by making use of the fact that scatterers are represented as block sparse vectors after discretization [START_REF] Dorn | Level set methods for inverse scattering[END_REF][START_REF] Wright | Sparse reconstruction by separable approximation[END_REF]. This a-priori information about the distribution of non-zero elements can be incorporated within the regularization scheme to increase the robustness of the solution and the quality of the recovered images.

According to some studies such as [START_REF] Li | A contrast source inversion method in the wavelet domain[END_REF][START_REF] Bevacqua | Non-linear inverse scattering via sparsity regularized contrast source inversion[END_REF][START_REF] Anselmi | Wavelet-based compressive imaging of sparse targets[END_REF], wavelet decomposition is an ecient transformation while retrieving scatterers which are sparse with respect to their expansion functions. In our specic cases, the proposed methods do not exhibit a remarkable dierence between wavelet and pixel basis functions. This might be due to the fact that the scattering matrix, which is obtained by sampling the Green function between investigation domain and receiver locations, does not satisfy the restricted isometry property [START_REF] Enyuva | Electromagnetic imaging of closely spaced objects using matching pursuit based approaches[END_REF][START_REF] Sandhu | A modied CoSaMP algorithm for electromagnetic imaging of two dimensional domains[END_REF]. Following the work done in [START_REF] Sandhu | A modied CoSaMP algorithm for electromagnetic imaging of two dimensional domains[END_REF], we have adapted our scattered eld by adding a constant parameter to its diagonal entries. However, we have not been able to get better results so that it can be a good idea to investigate the problem within this perspective in the future.

Although Haar and Daubechies wavelet basis functions beside pixel basis function are studied in this work to represent the unknown model in Sec. 3.3, the optimal choice of the basis set is an application-dependent open problem. The recovery problem solution is the sparsest one with respect to a basis function and the sparsest reconstruction (fewest non-zero coecients) has to be picked up. Therefore, the available a-priori (physical) solution information plays a fundamental role in the selection of the most suitable basis for a given application. Within this respect, more research can be carried out in order to deduce which basis function can be the most appropriate to specic applicative domains in a more practical and faster manner. Other options such as Curvelets [START_REF] Candes | Curvelets: A surprisingly eective nonadaptive representation for objects with edges[END_REF] are of interest to a future research also. Future works can be aimed at extending the validation of the proposed inversion scheme to other wavelet families more suitable for dealing with dierent applicative domains such as nondestructive testing/evaluation. Furthermore, more investigations can be performed for 3D electromagnetic imaging, where the nonlinearity and ill-posedness become even more severe. For instance, we have already started looking into this direction of study in which reconstruction of 3D defects aecting an anisotropic laminate is dealt with the two-step inversion method as stated in Sec. 4.2.4. We can obtain good localization and estimation of the contrast by applying joint sparsity through T-MSBL. However, this method is applicable for small scatterers so far and more work needs to be performed for larger scatterers. This method can also be studied for the localization and estimation of the contrast of a small defect in a complex anisotropic structure. Improvements with respect to the convergence speed and radius of the convergence of the method should be investigated as well. so that a small change in the cost function w.r.t the real part of contrast is described as δK Re S L Θ(r,r )δE di (r,r )drdr , (A.4) where overlined terms are denoted as complex conjugate ones [START_REF] Brandwood | A complex gradient operator and its application in adaptive array theory[END_REF].

Considering the terms forming the δK, δE di can be expressed as δE di (r) = D G or (r,r )δJ(r)dr . (A.9)

Therefore, the gradient of the discrepancy K w.r.t χ real (r) is given by K (χ real (r)) = 2 × Re(P (r)E tot (r s , r))

(A.10)
The same process can be followed for the gradient w.r.t χ imag (r) as well. It is known that dK(χ(r)) dχ(r) = dK(χ(r)) dχ real (r) -i dK(χ(r))

dχ imag (r) . Therefore, K (χ(r)) = 2 × P (r)E tot (r s , r)

(A.11)
It can be noted that the equation in (A.8) has the same form as equation [START_REF] Ramlau | A Tikhonov-based projection iteration for nonlinear ill-posed problems with sparsity constraints[END_REF]. Therefore, it can be solved through MoM in the same way. By discretization, we can obtain

P i = N R r=1 G ir Θ * k + N j=1
G ij χ j P j j = 1, 2, . . . , N (A.12)

where N R is the number of receivers, G ir is the Green`s function between the position of the ith point in the region of interest and the position of the kth receiver. The equation above can also be written as the following linear equations such as

      P 1 P 2 . . . P N       =       N R r=1 G 1r Θ * k N R r=2 G 2r Θ * k . . . N R r=N G N r Θ * k       +       G 11 χ 1 G 12 χ 2 • • • G 1N χ N G 21 χ 1 G 22 χ 2 • • • G 1N χ N . . . . . . . . . . . . G N 1 χ 1 G N 2 χ 2 • • • G N N χ N           
P N       =       1 -G 11 χ 1 -G 12 χ 2 • • • -G 1N χ N -G 21 χ 1 1 -G 22 χ 2 • • • -G 1N χ N . . . . . . . . . . . . -G N 1 χ 1 -G N 2 χ 2 • • • 1 -G N N χ N       -1       N R r=1 G 1r Θ * k N R r=2 G 2r Θ * k . . . N R r=N G N r Θ * k       (A.14) whereas       N R r=1 G 1r Θ * k N R r=2 G 2r Θ * k . . . N R r=N G N r Θ * k       =       G 11 G 12 • • • G 1N R G 21 G 22 • • • G 2N R . . . . . . . . . . . . G N R 1 G N R 2 • • • G N N R                   ξ 1 ξ 2 . . . ξ N R       -       E di 1 E di 2 .
. .

E di N R             * (A.15) Furthermore,       E di 1 E di 2 .
. .

E di N R       =       N i=1 G i1 χ i E tot i N i=1 G i2 χ i E tot i . . . N i=1 G iN R χ i E tot i       (A.16)       E di 1 E di 2 . . . E di N R       =       G 11 G 21 • • • G N R 1 G 12 G 22 • • • G N R 2 . . . . . . . . . . . . G 1N R G 2N R • • • G N N R             χ 1 0 • • • 0 0 χ 2 • • • 0 . . . . . . . . . . . . 0 0 • • • χ N             E tot 1 E tot 2 .
. .

E tot N R             E di 1 E di 2 .
. .

E di N R       =       G 11 G 12 • • • G 1N R G 21 G 22 • • • G 2N R . . . . . . . . . . . . G N R 1 G N R 2 • • • G N N R       †       χ 1 0 • • • 0 0 χ 2 • • • 0 . . . . . . . . . . . . 0 0 • • • χ N             E tot 1 E tot 2 .
. .

E tot N R       (A.17)
These two equations refer to the ltering operation made by lters h(k) and g(k) using the convolution products of f with their impulse responses respectively, and the factor 2k refers to the undersampling. Therefore, the discrete wavelet transform (DWT) can be summarized as follows:

f -→ (Gf, GHf, GH 2 f, • • • , GH j-1 f, H j f ) = (d (j-1) , d (j-2) , • • • , d (0) , c (0) )

(B.7)
where d (j-1) , d ( j-2), • • • , d (0) are the details and a ( 0) is the approximation coecient dened as: c (j-1) = Hc (j) and d (j-1) = Gd (j)

(B.8)
The reconstruction process is similar with the decomposition. For each decomposition level, the signal is oversampled by two and passed by the synthetic high pass and low pass lters Ḡ and H that will be added later. Therefore, we dene the operators Ḡ and H as follows:

( Hf

) k = n h(n -2k)f (n)( Ḡf ) k = n h(n -2k)f (n).
(B.9)

The recursive calculation gives: where D j and A are the details and the approximation, respectively.

f = n-1 j=0 (

1. 1 . 3 . 6 . 3 . 7 . 3 . 8 . 3 . 9 . 74 4. 2 .

 136373839742 (a) 1-D plot of x p , p = 0.5, 1, 2. (b) 2-D level sets of x p , same p. . . . . 2.1. Geometry of scattering experiment. . . . . . . . . . . . . . . . . . . . . . . 2.2. Validation of the model when N = 125 × 125. . . . . . . . . . . . . . . . . 2.3. Evolution of the normalized error of the total eld (left) and on the scattered eld (right) as a function of the discretization number N . . . . . . . 2.4. Example for 1 source and N = 29 receivers for various SNR. . . . . . . . . 2.5. Example for 1 source and N = 29 receivers for various 10 dB noise level. . . 2.6. Dierence of scattered eld obtained with and without Born approximation. 3.1. Comparison between weak monotonicity and back-tracking with Armijo condition, error on cost function (vertical axis) w.r.t iteration number (horizontal axis). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Performance of soft and hard thresholding. . . . . . . . . . . . . . . . . . . 3.3. Filter bank of the 2D wavelet transform [48]. . . . . . . . . . . . . . . . . .3.4. Measured conguration of actual permittivity proles for dierent scatterers and source-receiver locations. . . . . . . . . . . . . . . . . . . . . . . . 3.5. Error of permittivity w.r.t dierent values of α when β is equal to 0 for the cases with (w) sparsity, without (w/o) sparsity and equal to 1.5 for sparsity and smoothness (smoothing the gradient (Sec. 3.2.2)) . . . . . . . . . . . . Error of permittivity w.r.t dierent values of β when α is equal to 2.5 × 10 -2 with (w) sparsity and α is equal to 1 × 10 -6 without (w/o) sparsity. . . . . Error of permittivity w.r.t dierent values ofβ and α with 10 dB noise data. Scatterer 6 (Fig. 3.4f): Retrieval of permittivity ε r (left) and conductivity σ (right) using sparsity and smoothness without noise. . . . . . . . . . . . Scatterer 6 (Fig. 3.4f): Retrieval of permittivity ε r (left) and conductivity σ (right) using sparsity and smoothness with 10 dB noise data. . . . . . . . xi List of Figures 3.10. Scatterer 1 (Fig. 3.4a): Error in ε r and error in σ as a function of α by using sparsity with (w) and without (w/o) projection with 20 noise levels of 10 dB (vertical axis corresponds to an error and horizontal axis corresponds to a regularization parameter α). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.11. Scatterer 2 (Fig. 3.4b): Error in ε r and error in σ as a function of α by using sparsity with (w) and without (w/o) projection with 5 noise levels of 10 dB (vertical axis corresponds to an error and horizontal axis corresponds to a regularization parameter α). . . . . . . . . . . . . . . . . . . . . . . . 3.12. Scatterer 1 (Fig. 3.4a): Retrieval of permittivity ε r (top) and retrieval of conductivity σ (bottom) by using sparsity and projection with 10 dB noise data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.13. Scatterer 2 (Fig. 3.4b): Retrieval of permittivity ε r (top) and retrieval of conductivity σ (bottom) by using sparsity and projection with 10 dB noise data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.14. Scatterer 4 (Fig. 3.4d): Error in ε r and minimum cost function obtained by using sparsity and positivity with 20 noise levels of 10 dB data. . . . . . . . 3.15. Scatterer 3 (Fig. 3.4c): Retrieval of permittivity (ε r , top) and conductivity (σ, bottom) by using sparsity and projection with 10 dB noise data. . . . . 3.16. Scatterer 4 (Fig. 3.4d): Retrieval of permittivity (ε r , top) and conductivity (σ, bottom) by using sparsity and projection with 10 dB noise data. . . . . 3.17. Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.18. Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.19. Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.20. Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 5 × 10 -3 without projection and smoothness). . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1. Scatterer 2 (Fig. 3.4b): Retrieval of permittivity (ε r , left) and conductivity (σ, right) both in the spatial domain and the wavelet domain (WD) by using sparsity with 10 dB noise data. . . . . . . . . . . . . . . . . . . . . . Scatterer 6 (Fig.3.4f): Retrieval of permittivity (ε r , left) and conductivity (σ, right) both in the spatial domain and the wavelet domain (WD) by using sparsity with 10 dB noise data. . . . . . . . . . . . . . . . . . . . . . 75 4.3. Error in χ (χ err ) against varying SNR both in the spatial and the wavelet domain . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77 4.4. Conguration for inverse imaging testing. . . . . . . . . . . . . . . . . . . . 78 4.5. Convergence of the cost function, vertical axis corresponds to number of iterations and horizontal axis corresponds to error on cost function (twostep method applied via YALL1). . . . . . . . . . . . . . . . . . . . . . . . 79 4.6. Reconstruction of contrast ( χ) without noise through YALL1 . . . . . . . . 81 4.7. Reconstruction of contrast ( χ) with 30 dB noise through YALL1. . . . . . . 82 4.8. Reconstruction of real part of contrast (Re{ χ}) without noise through T-MSBL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.9. Reconstruction of imaginary part of contrast (Im{ χ}) without noise through T-MSBL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84 4.10. Reconstruction of real part of contrast (Re{ χ}) with 10 dB noise through T-MSBL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 4.11. Reconstruction of imaginary part of contrast (Im{ χ}) with 10 dB noise through T-MSBL. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86 B.1. Filter bank of the 2D wavelet transform. . . . . . . . . . . . . . . . . . . . 111 xiii 1. Introduction

  Figure 1.1.: (a) 1-D plot of x p , p = 0.5, 1, 2. (b) 2-D level sets of x p , same p.

Chapter 3 :

 3 Methods to solve Nonlinear Inverse ProblemWe start by giving a brief overview of discrete ill-posed inverse problems, and motivate 1. Introduction 1.3. Outline of thesis the need for regularization. Then, we describe an important nonlinear inverse problem, sparse representation of signals using suitable bases. This problem serves a central role in the thesis: the basis of our work is the transformation of the microwave imaging problem into the problem of sparse signal representation. In this chapter we describe numerical optimization of the objective functions corresponding to l 1 regularization.

Chapter 5 :

 5 Summary, Conclusion and Future WorkThis chapter summarizes the main ideas of the thesis and gives suggestions for further research in the area.

Figure 2 .

 2 Figure 2.1.: Geometry of scattering experiment.

Fig. 2 .

 2 Fig. 2.2 validates that with a high discretization (N = 125) our model is well-tted with the data. On the other hand, Fig. 2.3 shows the evolution of the normalized error on to the total eld and the scattered eld as a function of the discretization number N .It is also worthwhile to mention that we are not in Born approximation case as it can be seen in Fig.2.6. The scattered elds are calculated for the case of without Born approximation throughout this manuscript.

Figure 2 . 2 .

 22 Figure 2.2.: Validation of the model when N = 125 × 125.

Figure 2 . 3 .

 23 Figure 2.3.: Evolution of the normalized error of the total eld (left) and on the scattered eld (right) as a function of the discretization number N .

3 .

 3 Figure 3.1.: Comparison between weak monotonicity and back-tracking with Armijo condition, error on cost function (vertical axis) w.r.t iteration number (horizontal axis).

Figure 3

 3 Figure 3.2.: Performance of soft and hard thresholding.

(3. 43 )

 43 Figure 3.3.: Filter bank of the 2D wavelet transform [48].

(3. 44 ) 3 .

 443 Methods to Solve Nonlinear Inverse Problem 3.3.2. Forward Problem in the Wavelet Domain Using the same basis decomposition as in (2.21) along r we can write

Figure 3 . 4 .

 34 Figure 3.4.: Measured conguration of actual permittivity proles for dierent scatterers and source-receiver locations.

,Figure 3 . 5 .

 35 Figure 3.5.: Error of permittivity w.r.t dierent values of α when β is equal to 0 for the cases with (w) sparsity, without (w/o) sparsity and equal to 1.5 for sparsity and smoothness (smoothing the gradient (Sec. 3.2.2))

Figure 3 . 6 .

 36 Figure 3.6.: Error of permittivity w.r.t dierent values of β when α is equal to 2.5 × 10 -2with (w) sparsity and α is equal to 1 × 10 -6 without (w/o) sparsity.

Figure 3 . 7 .

 37 Figure 3.7.: Error of permittivity w.r.t dierent values of β and α with 10 dB noise data.

Figure 3 . 8 .

 38 Figure 3.8.: Scatterer 6 (Fig. 3.4f): Retrieval of permittivity ε r (left) and conductivity σ (right) using sparsity and smoothness without noise.

Figure 3 . 9 .

 39 Figure 3.9.: Scatterer 6 (Fig. 3.4f): Retrieval of permittivity ε r (left) and conductivity σ (right) using sparsity and smoothness with 10 dB noise data.

Figure 3 .

 3 Figure 3.10.: Scatterer 1 (Fig. 3.4a): Error in ε r and error in σ as a function of α by using sparsity with (w) and without (w/o) projection with 20 noise levels of 10 dB (vertical axis corresponds to an error and horizontal axis corresponds to a regularization parameter α).

Figure 3 .

 3 Figure 3.11.: Scatterer 2 (Fig. 3.4b): Error in ε r and error in σ as a function of α by using sparsity with (w) and without (w/o) projection with 5 noise levels of 10 dB (vertical axis corresponds to an error and horizontal axis corresponds to a regularization parameter α).

Figure 3 .

 3 Figure 3.12.: Scatterer 1 (Fig. 3.4a): Retrieval of permittivity ε r (top) and retrieval of conductivity σ (bottom) by using sparsity and projection with 10 dB noise data.

Figure 3 .

 3 Figure 3.13.: Scatterer 2 (Fig. 3.4b): Retrieval of permittivity ε r (top) and retrieval of conductivity σ (bottom) by using sparsity and projection with 10 dB noise data.

Figure 3 .

 3 Figure 3.15.: Scatterer 3 (Fig. 3.4c): Retrieval of permittivity (ε r , top) and conductivity (σ, bottom) by using sparsity and projection with 10 dB noise data.

Figure 3 .

 3 Figure 3.16.: Scatterer 4 (Fig. 3.4d): Retrieval of permittivity (ε r , top) and conductivity (σ, bottom) by using sparsity and projection with 10 dB noise data.

3 . 4 .

 34 In particular, we study the Stationary Wavelet Transform (SWT) in terms of Haar basis of level two and three, and the Discrete Wavelet transform (DWT).

Figure 3 .

 3 Figure 3.17.: Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness).

Figure 3 .

 3 Figure 3.18.: Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness).

Figure 3 .

 3 Figure 3.19.: Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 2.5 × 10 -2 without projection and smoothness).

  Figure 3.20.: Retrieval of permittivity (ε r , 1st and 3rd column) and conductivity (σ, 2nd and 4th column) with 20 dB and 10 dB noise data (α = 5 × 10 -3 without projection and smoothness).

3 .

 3 Methods to Solve Nonlinear Inverse Problem 3.5. Conclusion Two imaging techniques have been introduced in this Chapter in order to solve the nonlinear inverse scattering problem. The rst technique is based on the iterative algorithm of soft shrinkage type which enforces the sparsity constraint at each nonlinear iteration. We have introduced constraints such as Barzilai and Borwein (BB) step size selection criteria and projection constraint on contrast in order to favor the performance and quality of the proposed algorithm.

( 4 . 3 )

 43 based on the original joint sparsity model (4.2) can be recast as a minimization problem under the constraints that J i has the same sparse support for each source i leading to the use of a (weighted) l 2,1 -regularization to enforce joint sparsity such as min J

  applied the alternating direction method (ADM) technique to solve the l 1 problem in compressed sensing and developed the corresponding MATLAB package termed Your ALgorithms for L1 (YALL1). Furthermore, Deng et al. extended the YALL1

Figure 4 . 1 .

 41 Figure 4.1.: Scatterer 2 (Fig. 3.4b): Retrieval of permittivity (ε r , left) and conductivity (σ, right) both in the spatial domain and the wavelet domain (WD) by using sparsity with 10 dB noise data.

Figure 4 . 2 .

 42 Figure 4.2.: Scatterer 6 (Fig.3.4f): Retrieval of permittivity (ε r , left) and conductivity (σ, right) both in the spatial domain and the wavelet domain (WD) by using sparsity with 10 dB noise data.

Figure 4 . 3 .

 43 Figure 4.3.: Error in χ (χ err ) against varying SNR both in the spatial and the wavelet domain

  Figure 4.4.: Conguration for inverse imaging testing.

1

 1 cm • Acquisition: Surface area of L x × L y at height L acquisition L x = L y = 6λ 1 , L acquisition = 0.3λ 0 Number of antennas: N x × N y ; N x = N y = 13Each antenna: 3 orthogonal unit dipoles along x, y and z With 30 dB noise data

Figure 4 . 5 .

 45 Figure 4.5.: Convergence of the cost function, vertical axis corresponds to number of iterations and horizontal axis corresponds to error on cost function (two-step method applied via YALL1).

  4.8 and Fig. 4.9)

Figure 4 . 6 . 30 dBFigure 4 . 7 .

 463047 Figure 4.6.: Reconstruction of contrast ( χ) without noise through YALL1

Figure 4 . 8 .

 48 Figure 4.8.: Reconstruction of real part of contrast (Re{ χ}) without noise through T-MSBL.

Figure 4 . 9 .

 49 Figure 4.9.: Reconstruction of imaginary part of contrast (Im{ χ}) without noise through T-MSBL.

Figure 4 .

 4 Figure 4.10.: Reconstruction of real part of contrast (Re{ χ}) with 10 dB noise through T-MSBL.

Figure 4 .

 4 Figure 4.11.: Reconstruction of imaginary part of contrast (Im{ χ}) with 10 dB noise through T-MSBL.

(A. 5 ) 7 )

 57 where δJ(r) can be also dened asδJ(r) = δχ(r)E inc + δχ(r) D G oo (r, r)χ(r)E tot (r)dr + χ(r) D G oo (r, r)δJ(r)dr = δχ(r)E tot (r) + χ(r) D G oo (r, r)δ(χ(r)E tot (r))dr = δχ(r)E tot + χ(r) D G oo (r, r)δJ(r)dr (A.6) Considering G oo (r,r ) = G oo (r , r) T G ro (r,r ) = G or (r ,r) T and substituting A.5 into A.4 we can write δK as follows By denition as in [108], the adjoint is P (r) = L G ro (r, r )Θ(r)dr + D G oo (r, r )χ(r)P (r )dr (A.8) so that substituting A.8 into A.7, δK becomes δK = Re S D δJ(r). P (r) -D G oo (r,r )χ(r )P (r )dr dr = Re S D P (r).δJ(r) -δJ(r) D G oo (r,r )χ(r)P (r )dr dr = Re S D P (r) δJ(r) -χ(r) D G oo (r,r )δJ(r )dr dr = Re S D P (r, r s )δχ(r)E tot (r, r s )dr .

  H) j Ḡd j + ( Hf ) n c (0) (B.10) or in the temporal domain: D j = ( H) j Ḡd j and A = ( H) n a (0) (B.11)

  ) where τ is the step size, ζ (χ) is the gradient of the nonlinear function ζ(χ) with respect to χ, and ζ * (χ) is the adjoint of the operator ζ (χ). S α is the soft shrinkage operator dened componentwise by

  3. Methods to Solve Nonlinear Inverse Problem 3.3.1. Review of wavelet representation and the wavelet transform 3.3.1. Review of wavelet representation and the wavelet transform Any continuous function in R 2 can be approximated in terms of scaling and wavelet functions up to some decomposition level J according to the following formula[START_REF] Li | A contrast source inversion method in the wavelet domain[END_REF]

Table 3 .

 3

1.: Description of the scatterers, x, y being the coordinate of the center of the obstacle (in m), L x , L y its lengths (in m) and ε r and σ its relative permittivity and conductivity (the latter in S m -1 )

Table 3 .

 3

	2.: Properties of the geometry of each system: Description of the region of in-
	terest (ROI), frequency, discretization, and number of sources (Nsources) and
	receivers (Nreceivers).

Table 3 .

 3

	3.: Error on permittivity with or without addition of projection

  The vectors x 1 , • • • , x s share the sparsity pattern M , i.e., the nonzero entries of x 1 , • • • , x s appear at the same positions. A common recovery model is

	min X	|M | s.t y j

Table 4 .

 4

2.: Error and average simulation time in seconds with 10 dB-Wavelet Domain

n, des méthodes de régularisation linéaires et non linéaires, qui minimisent une fonction de coût pondérée entre mesure/inadéquation des données et un terme de pénalité, ont été développées. Le choix le plus populaire du terme de pénalité est la norme l 2 de la solution[START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. Le problème de minimisation qui en résulte peut être résolu en utilisant le schéma bien connu de Tikhonov[START_REF] Tveito | Parameter estimation and inverse problems[END_REF][START_REF] Fornasier | Theoretical Foundations and Numerical Methods for Sparse Recovery[END_REF]. En outre, l'utilisation d'itérations tronquées de Landweber ou de gradient conjugué conduit à un type similaire de régularisation[START_REF] Pastorino | Microwave Imaging[END_REF][START_REF] Tveito | Parameter estimation and inverse problems[END_REF].En programmation linéaire, toutes ces méthodes ltrent ecacement les valeurs singulières du problème de minimisation discrétisé, qui sont inférieures à un niveau de seuil spécié, ou font appel à un paramètre de régularisation pour atténuer le caractère mal-posé.La régularisation est utilisée pour résoudre des problèmes mal-posés en incorporant des connaissances a-priori sur l'objet de façon à stabiliser le problème et fournir des solutions raisonnables et utiles. Par exemple, si l'on sait que la solution doit être une discrétisation d'une fonction continue, cette connaissance nous permet de rejeter les candidats les plus variables et de réduire considérablement l'ensemble des solutions possibles. Le choix d'une régularisation appropriée dépend étroitement de la propriété que l'on souhaite mettre en avant, et cela dépend de l'application particulière. Dans de nombreux problèmes inverses mathématiques, les choix a-priori sont des formes diérentes de régularités, et les pénalités correspondantes sont les normes l 2 de x ou de ses dérivése. Des a priori de parcimonie sont utiles lorsque les signaux que nous recherchons doivent être parcimonieux. Nous dénissons la parcimonie d'un vecteur par la présence d'un petit nombre d'éléments non nuls comme dans le travail de Mallat[START_REF] Mallat | A Wavelet Tour of Signal Processing[END_REF], celui de v

(Fig.3.4d) has lower error in permittivity compared to its full view case (Fig.3.4f) as

Université Paris-Saclay Espace Technologique / Immeuble Discovery Route de l'Orme aux Merisiers RD 128 / 91190 Saint-Aubin, France

I would like to thank my supervisors Dr. Marc Lambert, chargé de recherche at Centre and Prof. Andrea Massa, director of the years. Additionally, I would like to thank DIGITEO France for nancing of this thesis through the SIRENA project (20142017) under the Call for Chairs 2014. I would like also to thank Dominique Lesselier, Directeur de Recherche at CNRS, and his

A. Adjoint Gradient of Cost Function

A.1. Adjoint Gradient of Cost Function

The discrepancy D under consideration is the standard least squared cost function such as

where ζ(r) is the measured scattered eld at the receiver position L while E di is the scattered eld obtained analytically and S is the source position. By multiplying both sides of the state equation by χ(r), the following equation J(r) = χ(r)E inc + χ(r) Ω G(r,r )χ(r)E tot (r,r )dr

is obtained where J(r) = χ(r)E tot (r) and χ(r) = χ real (r) + χ imag (r). Taking the rst-order development of the kind K(χ(r) + δχ real (r):

B.1. The Wavelet Transform

The wavelet transform has been proposed to overcome the limits of the Fourier transform.

However, if we apply the Fourier transform on a given signal, we cannot determine the moments corresponding to high frequencies. The Short Fourier Transform (SFT) uses a windowing technique to calculate the spectrogram, which provides information in the time-frequency domain.

The problem of the SFT is the Heisenberg uncertainty principle [START_REF] Busch | Heisenberg's uncertainty principle[END_REF]. This principle states that it is impossible to know the exact time frequency correspondence of a signal, i.e., we cannot know which spectral component corresponds to a signal at given instant. Although the temporal and frequency resolution problems come from physical phenomenon and they are independent of the used transform, it is possible to analyze any signal using an alternative approach called wavelet transform (WT). This transform allows to analyze the signals at dierent frequencies resolutions unlike the Fourier transform that perform a unique frequency resolution analysis. The wavelet transform is designed to give a poor temporal resolution and good frequency resolution at high frequencies and vice-versa at low frequencies. This compromise makes sense especially when the analyzed signal has short duration high frequency components and long duration low frequencies components. This is the case of most biological and natural signals [START_REF] Daubechies | Ten Lectures on Wavelets[END_REF]. So the major dierence between the Fourier and wavelet transforms is that the window length varies with the transform computed for each spectral component.

B.1.1. Approximation theory and multiresolution analysis

The theory of multiresolution analysis in wavelet domain is founded by Mallat in 1989 [START_REF] Mallat | A theory for multiresolution signal decomposition: the wavelet representation[END_REF]. Its principle consists to assume that there is a series of closed subsets V j j∈Z that approximate the integrable L 2 (R) functions. The subsets V j are designed in such a way that

. Then, for a given function f ∈ L 2 (R), there are f N ∈ V N approximating f . If we set g i ∈ W i and f i ∈ V i then we can write:

This is the wavelet decomposition of the function f . In the signal processing context, this decomposition is established under a pyramidal coding algorithm. The time frequency representation of a signal in wavelet domain is obtained using ltering techniques. The series h(k), and g(k), k ∈ Z are quadrature mirror lters and the relationship between h and g is given by:

The series h(k) is a family of low pass lters and g(k) are high pass lters. These two lters families belong to the Finite Impulse Response (FIR) lters. By using the Fourier transform and the orthogonality property between h and g, we can prove that [START_REF] Graps | An introduction to wavelets[END_REF]:

The wavelet decomposition process begins by passing the signal through the low frequency band h(n). The ltering operation of the signal consists to calculate its convolution with the impulse response of the lter. Note that a half band low pass lter removes all frequencies that are higher than half of the maximum frequency in the signal. So after applying the low pass lter, half of the samples will be eliminated and the ltered signal will have only half of samples in the next decomposition level. On other hand, the low pass ltering divides the resolution of the signal by two since this resolution depends only on the quantity of information in the signal. Then, the signal is undersampled by two since redundancy of the half of the samples and this operation doubles the scale.

B.1.2. Continuous Wavelet Transform

The continuous wavelet transform coecients applied to a given signal f are dened as follows [START_REF] Graps | An introduction to wavelets[END_REF]:

is the window function of the WT, a is the scale parameter and b is the translation parameter.

The translation coecient b is related to the window location, it corresponds to the temporal information in the wavelet domain. On the other hand, we have not the same parameter for the frequency information, like in the Fourier transform, but the scale parameter a. This latter is used to delay or compress the transformed signal. The small scale values correspond to a dilation and the large values correspond to a compression.

Note that the relationship between the scale and the frequency is reverse, i.e., low scale values correspond to high frequencies and vice-versa. To express the multiband coding process in the wavelet domain, let us set the series f = f n that represents the discrete signal to be decomposed and the operators H and G dened as follows: Title : Sparsity and electromagnetic imaging in non-linear situations Keywords : electromagnetic imaging, non-linear, sparsity, wavelet decomposition Abstract : So-called quantitative electromagnetic imaging focused onto here is the problem of determining material properties from scattered elds measured away from the domain under investigation. Solving this inverse problem is a challenging task because it is ill-posed due to the presence of (smoothing) integral operators used in the representation of scattered elds in terms of material properties, and scattered elds are obtained at anite set of points through noisy measurements. Moreover, the inverse problem is nonlinear simply due the fact that scattered elds are nonlinear functions of the material properties. The work described in this thesis deals with the ill-posedness of the electromagnetic imaging problem using sparsity-based regularization techniques, which assume that the scatterer(s) capture only a small fraction of the investigation domain and/or can be described in sparse fashion on a certain basis.

B.1.3. Multiband coding in wavelet domain

The primary aim of the thesis is to intensively investigate sparsity regularization for nonlinear inverse problems. Therefore, we focus on sparsityregularized nonlinear Tikhonov method which directly solves the nonlinear minimization problem using Landweber iterations, where a thresholding function is applied at every iteration step to promote the sparsity constraint. This scheme is accelerated using a projected steepest descent method and replaces the thresholding operation to enforce the sparsity constraint. This approach has also been implemented in wavelet domain which allows an accurate representation of the unknown function with a reduced number of coecients. Additionally, we investigate a method correlated with the joint sparsity which gives multiple sparse solutions that share a common nonzero support in order to solve concerned nonlinear problem.