Modélisation physique de la diffraction d’ondes électromagnétiques

Télécharger le fichier pdf d’un mémoire de fin d’études

Comparaison aux données réelles de l’Institut Fresnel

Nous avons également souhaité valider le modèle en s’appuyant sur des données réelles. L’Institut Fresnel (Marseille, France) propose librement un ensemble de fichiers de mesures issues d’expérimentations réalisées en chambre anéchoïque. Ces données ont été présentées au monde scientifique dans l’objectif de valider les méthodes d’inversion et de reconstruire des objets de diﬀérentes formes, tailles et propriétés diélectriques. Toutes les informations nécessaires à la reproductibilité des expérimentations sont données dans Geﬀrin et Sabouroux (2009). Le chapitre 6 se focalise sur la reconstruction de ces objets et les détails du montage d’acquisition y sont donnés.
Pour valider notre modèle, nous avons reproduit les expériences de l’Institut Fresnel et com-paré les mesures expérimentales aux valeurs obtenues par le modèle numérique. Ici, nous donnons uniquement les résultats pour l’objet TwoCubes. La figure 2.5 est tirée de Geﬀrin et Sabouroux (2009) et présente l’objet. L’objet est composé de deux cubes purement diélec-triques de 25 mm de côté et d’une permittivité relative de 2.35. Le volume d’étude englobant les deux cubes est discrétisé avec une résolution de 2.5 mm. Les figures 2.6 et 2.7 présentent les parties réelle et imaginaire des données aux 243 récepteurs pour les 36 sources aux fré-quences 5 GHz et 8 GHz. On observe que les valeurs du champ diﬀracté obtenues par le modèle numérique sont très proches des données expérimentales.
Pour mesurer quantitativement ces erreurs, nous avons calculé l’erreur quadratique. Contrai-rement aux comparaisons avec la solution analytique de Mie, nous comparons ici à des données expérimentales. Les données étant bruitées, il est important de comprendre que nous ne cherchons pas à obtenir une erreur nulle mais plutôt à obtenir une variance faible entre le modèle et les mesures expérimentales. D’après les chercheurs de l’Institut Fresnel, le modèle numérique développé au sein de leur équipe aﬃchait une erreur entre 5 et 10% par rapport aux mesures expérimentales pour les diﬀérentes fréquences utilisées. Le tableau 2.4 donne les erreurs pour les diﬀérentes mesures. Malgré quelques mesures où les erreurs atteignent 17-18%, la plupart des erreurs sont comprises entre 5 et 10%. Les erreurs les plus impor-tantes apparaissent pour les plus faibles fréquences et pour les plus hautes fréquences. Une explication possible de ces erreurs pour les faibles fréquences est que les valeurs du champ diﬀracté sont faibles alors que le bruit de mesure est constant pour toutes les fréquences. Il semble alors normal que ces forts pourcentages correspondent à la variance du bruit de me-sure. Pour les hautes fréquences, une explication éventuelle est qu’à cette fréquence, l’objet est de grande taille électrique et la résolution est grande par rapport à la longueur d’onde. Tout comme pour la solution analytique, le modèle numérique pourrait montrer ses limites, ce qui justifierait ces erreurs plus élevées que pour les données aux autres fréquences.

MÉTHODES D’INVERSION EN IMAGERIE MICRO-ONDES

Le chapitre précédent a présenté le modèle direct permettant de calculer le champ électrique diﬀracté par un objet diélectrique. Ce chapitre se concentre sur le problème inverse et sur les méthodes de reconstruction. Le problème inverse en TMO consiste en l’estimation des propriétés diélectriques d’un milieu inconnu à partir de mesures de champs diﬀractés. Cette estimation repose généralement sur la minimisation d’un critère d’adéquation aux données : on cherche les propriétés électriques pour lesquelles l’erreur entre les observations et le mo-dèle numérique est la plus faible. Les problèmes inverses sont souvent mal posés au sens de Hadamard et c’est le cas du problème de TMO (Pastorino, 2010). Il est alors fréquent d’introduire un terme de régularisation permettant de stabiliser la solution. Enfin, un algo-rithme d’optimisation doit être choisi afin de trouver un optimum au critère à minimiser. Finalement, trois éléments sont primordiaux dans la définition de notre problème inverse en TMO : le critère, la régularisation et l’algorithme d’optimisation. Ce chapitre propose de comparer les diﬀérentes approches fréquentes dans la littérature pour ces trois éléments et nous présenterons notre méthode d’inversion en justifiant nos choix.

Ce chapitre s’articule autour de cinq sections :

– la première section présente le problème général d’imagerie micro-ondes ;

– la deuxième section traite de l’inversion vue comme un problème d’optimisation. Il s’agit de présenter les fonctions coût fréquemment utilisées dans la littérature et de comparer les avantages et inconvénients de chacune ;

– les diﬀérents types de régularisation sont discutés dans la troisième section et nous décri-vons en détail le choix adopté pour notre problème ;

– la section 4 dépeint les diﬀérentes familles d’algorithmes d’optimisation utilisées pour la reconstruction d’images en TMO et compare les propriétés de chacune ;

– enfin, la dernière section présente le calcul du gradient de la fonction coût choisie. Ce calcul met en évidence les diﬃcultés de la résolution du problème inverse.

L’imagerie micro-ondes

L’imagerie micro-ondes est un ensemble très vaste de techniques de détection, de localisation et d’estimation d’objets dissimulés à l’aide d’ondes électromagnétiques aux fréquences micro-ondes (300 MHz – 300 GHz). Cette science regroupe une gamme variée de problèmes et d’applications. Certaines méthodes d’imagerie ne se focalisent que sur la détection d’objets particuliers cachés, c’est-à-dire de déceler la présence ou non d’un corps ; d’autres méthodes recherchent la position et/ou la forme de ces objets. Enfin, des méthodes plus complexes cherchent à estimer aussi bien la position, la forme et les propriétés diélectriques des objets, ce qui revient à réaliser une cartographie des propriétés diélectriques dans un volume d’intérêt. Dans ce manuscrit, nous nous intéressons uniquement à ces dernières, les méthodes dites quantitatives. Dans la littérature, l’expression inverse scattering methods est souvent utilisée pour nommer ce type d’approches. Un grand nombre de méthodes de reconstruction a été développé depuis quelques décennies mais il n’est pas envisageable de mentionner toutes ces approches. Ici nous nous concentrerons sur les méthodes basées sur la minimisation d’un critère pénalisé par des algorithmes itéra-tifs d’optimisation locale. Ces méthodes sont les plus courantes dans la littérature car elles donnent des résultats convenables pour un coût de calcul raisonnable. Plus d’informations sont données dans les sections suivantes. Nous pouvons citer malgré tout quelques méthodes “diﬀérentes” de celles que nous verrons par la suite, s’appuyant par exemple sur :

– diﬀérents modèles numériques. Nous avons parlé de la discrétisation des équations régissant la propagation d’onde et des méthodes numériques comme la méthode des moments, la mé-thode des éléments finis. Mais il existe des modèles bien diﬀérents comme celui basé sur la déformation d’un volume où un nombre fini de points de la surface fait oﬃce d’inconnues (Li et al., 2012), ou celui basé sur une représentation binaire du volume d’intérêt (Duchêne, 2001) ;

– diﬀérentes approximations. Certaines méthodes utilisent des approximations linéaires afin de contourner la non-linéarité du modèle direct et donc d’accélérer la reconstruction. C’est le cas de l’approximation de Born ou de Rytov, ou encore l’approche Linear Sampling Method (LSM) ;

– diﬀérents outils mathématiques, par exemple ceux permettant de résoudre un problème linéaire comme la troncature de la décomposition en valeurs singulières (TSVD).

Pour plus de références, Pastorino (2010) fournit de nombreuses descriptions d’algorithmes de reconstruction.

Pour la suite, on pose les notations suivantes : on suppose un volume homogène V dans lequel un ou plusieurs objets inconnus sont placés. Le volume est discrétisé en N voxels cubiques. Ce volume est illuminé successivement par NS ondes incidentes notés einci (i = 1, . . . , NS). NR récepteurs sont situés autour de ce volume et mesurent le champ électrique diﬀracté noté escati (i = 1, . . . , NS). Le vecteur escati de longueur NR contient les mesures aux NR récepteurs pour l’illumination i. On cherche alors à reconstruire la fonction de contraste sur les N voxels, représentée par le vecteur x de longueur N. On rappelle les équations matricielles d’observation et du domaine issues du chapitre précédent : escat = GoXetot (équation d’observation) etot = einc + GcXetot (équation du domaine), ou de façon équivalente, l’équation non-linéaire : escat = GoXL−x1einc avec Lx = I − GcX.

L’inversion vue comme un problème d’optimisation

Puisqu’à partir des mesures de champs diﬀractés, on ne peut pas directement obtenir une estimation du contraste, on considère le problème inverse comme la minimisation d’un critère permettant de réduire l’erreur entre les données mesurées et le modèle direct. En problème inverse, un critère des moindres carrés, c’est-à-dire l’erreur quadratique entre les mesures et le modèle, est souvent considéré. L’estimation du contraste correspond alors au minimiseur de ce critère qui est généralement obtenu par des algorithmes d’optimisation itératifs. S’ap-puyant sur le modèle numérique vu précédemment, deux formulations du critère ressortent fréquemment dans la littérature : une formulation bilinéaire et une formulation non-linéaire. Notons que ces formulations du critère ne sont pas réservées à une forme intégrale du modèle direct, elles peuvent être appliquées à une approche par diﬀérences finies ou éléments finis par exemple.

Formulation bilinéaire

On peut voir le problème inverse comme un problème d’optimisation de l’équation d’observa-tion sous la contrainte de l’équation du domaine. Il est en général plus simple de se ramener à un problème sans contrainte où la contrainte d’égalité est introduite d’une certaine façon dans le critère. La formulation bilinéaire s’appuie justement sur l’équation d’observation pour le terme d’attache aux données pénalisée par l’erreur commise sur la contrainte. Dans cette équation, aussi bien le contraste x que les champs totaux etoti sont inconnus : on cherche alors à les estimer conjointement. Dans la formulation bilinéaire, on ne cherche pas à résoudre exac-tement les champs totaux dans l’équation du domaine ; on relâche la contrainte d’égalité sur cette équation. Dans la littérature, ces méthodes prennent le nom MGM (Modified Gradient Method) (Kleinman et van den Berg, 1992) ou CSI (Contrast Source Inversion) (van de Berg et Kleinman, 1995). Le critère à minimiser prend la forme générique suivante : CSI tot tot 1 NS X F (x, e1 , . . . , eNS ) = 2 i=1 γ NS kescati − GoXetotik2 + 2 X ketoti − GDXetoti − eincik2. (3.4) i=1

Le paramètre γ permet de contrôler la pénalisation que l’on souhaite attribuer à la contrainte sur les champs totaux : plus γ est élevé, plus l’erreur sur la contrainte dans le second terme va tendre vers zéro. Il serait nécessaire d’avoir un poids γ très grand afin de s’assurer que le champ total estimé correspond bien à l’équation etoti = L−x1einci. De nombreux auteurs ont proposé de faire varier γ au cours des itérations de l’algorithme d’optimisation (van den Berg et Kleinman, 1995). En revanche, Barrière et al. (2011) a montré qu’il était judicieux de fixer empiriquement la valeur de γ et de garder cette valeur au cours de la reconstruction afin d’assurer la convergence de l’algorithme d’inversion.

Ces approches bilinéaires sont fréquemment utilisées car elles ont l’avantage majeur que le critère est linéaire par rapport à x et par rapport aux champs totaux. De plus, il n’est pas nécessaire de résoudre le problème direct (inversion par Lx) ; le coût de calcul par itération est considérablement réduit. En revanche, le nombre d’inconnues est très conséquent :(NS +1)×N inconnues. A cause du nombre important d’inconnues, il est souvent nécessaire d’eﬀectuer un grand nombre d’itérations pour converger.

Les algorithmes de reconstruction utilisant cette formulation bilinéaire du critère sont très fréquents notamment en 2-D car le nombre d’inconnues reste raisonnable (Abubakar et al., 2005; Barrière et al., 2011; van den Berg et Abubakar, 2001; Bozza et Pastorino, 2009; Gilmore et al., 2009b,a; Li et al., 2013; Mojabi et LoVetri, 2010; Ozdemir, 2014; Rubaek et al., 2011; Zakaria et LoVetri, 2012). En revanche, les méthodes d’inversion basées sur cette formulation pour des problèmes tridimensionnels sont bien plus rares. On peut tout de même en citer quelques-unes (Abubakar et van den Berg, 2004; Zhang et Liu, 2004; Catapano et al., 2006; Asefi et al., 2015).

Formulation non-linéaire

Contrairement à la formulation bilinéaire où la contrainte d’égalité est exprimée sous la forme d’un terme de pénalisation, la formulation non-linéaire s’appuie sur l’inclusion de cette contrainte dans le modèle direct. Cela revient donc à considérer l’équation non-linéaire (3.3). Le critère de cette formulation exprime directement l’erreur entre les données et le modèle non-linéaire. Le critère d’adéquation aux données est alors défini par : FD(x) = 1 NS (3.5) 2keiscat − GoXLx−1eiinck2 X i=1.

Les champs totaux ne sont plus des inconnues du problème, seule la fonction de contraste x en est une : le nombre d’inconnues est réduit à N. En revanche, le calcul du critère nécessite la résolution du problème direct, c’est-à-dire de la résolution des systèmes linéaires impliquant Lx. Cela entraîne un coût de calcul conséquent à chaque itération. De plus, la forte non-linéarité du modèle direct peut engendrer la présence de multiples minima locaux et rendre la résolution du problème inverse diﬃcile. En conclusion, le coût d’une itération du problème inverse est nettement plus important que celui pour le critère bilinéaire. En revanche, le nombre d’itérations est en général bien plus faible que pour l’approche bilinéaire puisque le nombre d’inconnues est bien plus faible.

En TMO 2-D, l’utilisation d’un critère non-linéaire a été utilisée (Roger, 1981; Chew et Wang, 1990; Remis et van den Berg, 2000; Franchois et Tijhuis, 2003; Fang et al., 2004; Hu et al., 2009). Pour des problèmes 3-D de grande taille, le fait que le nombre d’inconnues soit plus faible que pour les approches bilinéaires a incité les auteurs à se tourner davantage vers la formulation non-linéaire. Il s’agit de l’approche la plus fréquente dans la littérature ces dernières années (Harada et al., 2001; Soleimani et al., 2006; Yu et al., 2009; De Zaeytijd et Franchois, 2009; Chaumet et Belkebir, 2009; Winters et al., 2010; Grzegorczyk et al., 2012; Abubakar et al., 2012b,a; Estatico et al., 2013; Scapaticci et al., 2015). Notre choix s’est donc porté sur la formulation non-linéaire.

De nombreuses autres méthodes s’appuient sur diﬀérents critères. Ce manuscrit ne prétend pas comparer toutes ces méthodes. Nous pouvons citer tout de même les méthodes itératives basées sur des approximations du modèle direct, comme la méthode de Born itérative (Ali et Moghaddam, 2010; Zhang et Liu, 2015). Mudry et al. (2012) a proposé une méthode hybride mêlant les critères non-linéaire et bilinéaire. Des informations supplémentaires peuvent être obtenues dans Pastorino (2010).

La régularisation

Les problèmes inverses sont en général mal posés et c’est le cas en tomographie micro-ondes (Pastorino, 2010). Un problème bien posé respecte les trois conditions suivantes (Idier, 2001) :

– existence d’une solution,

– unicité de la solution,

– robustesse de la solution au bruit.

Obtenir ces trois conditions est très rare en pratique. Par exemple, la résolution d’un problème sous-déterminé donne une infinité de solutions possibles et l’unicité de la solution n’est donc pas vérifiée. C’est le cas en TMO, notamment en 3-D, où le nombre de voxels N est bien supérieur au nombre de données.

Pour éviter le caractère mal posé du problème en TMO, il est nécessaire d’ajouter des contraintes ou des informations a priori sur la solution. L’utilisation d’un terme de régu-larisation est une méthode très fréquente pour améliorer la stabilité de l’inversion et donc obtenir une solution correcte (Idier, 2001). Cette régularisation intervient comme un terme supplémentaire dans le critère à minimiser. Dans la suite, le terme de régularisation est noté R(x).

Régularisation “`2`1”

Diﬀérentes questions se posent quant au choix de la régularisation. Quelle connaissance a priori a-t-on sur l’objet ? Que cherche-t-on à pénaliser ? Quelle fonction de régularisation utiliser ? Comment intégrer le terme de régularisation dans le critère à minimiser ? Le choix de la régularisation reste ouvert et la quantité impressionnante d’approches de régularisation dans la littérature en est la preuve. Malgré tout, une régularisation conservant les discon-tinuités (dite « edge-preserving ») semble mettre d’accord la plupart des auteurs. Ce type de régularisation fait l’hypothèse que l’objet à reconstruire est un objet homogène par morceaux, ce qui en eﬀet représente la majorité des objets. Par exemple, pour une application biomédi-cale, on peut faire l’hypothèse que les propriétés diélectriques d’un tissu humain (os, muscle, etc.) sont constantes dans un même tissu.

L’idée est alors de pénaliser les variations de contraste, c’est-à-dire les diﬀérences premières : on pénalise fortement les faibles variations afin de favoriser des milieux constants et on pénalise dans une moindre mesure les fortes variations de contraste afin de permettre des discontinuités entre les milieux homogènes. De nombreuses fonctions de régularisation per-mettent ce type de pénalisation. Une liste est fournie dans Idier (2001). Nous avons fait le choix d’utiliser ce type de pénalisation dans cette thèse à l’aide des fonctions “`2`1” sur les diﬀérences premières. Ces fonctions ont l’avantage d’être continûment diﬀérentiables et donc d’avoir un gradient défini pour tout x. De plus, ces fonctions sont convexes, ce qui permet d’assurer la robustesse à la solution et d’améliorer la convergence des algorithmes d’optimi-sation. Le nom “`2`1” provient du fait qu’elles ont un comportement quadratique à l’origine et asymptotiquement linéaire à l’infini. Elles conviennent donc pour pénaliser fortement les faibles variations et préserver les grandes discontinuités. Il existe d’autres types de régu-larisation pour conserver les discontinuités comme les fonctions “`2`0” (asymptotiquement constantes à l’infini) mais celles-ci ne sont pas convexes, ce qui peut entraîner des instabilités lors de l’optimisation.

Le terme de régularisation utilisé dans nos travaux est le suivant : Nx−1 Ny−1 Nz−1 X X X R(x) = ϕ(χn+1,m,p −χn,m,p)+ϕ(χn,m+1,p −χn,m,p)+ϕ(χn,m,p+1 −χn,m,p). (3.6) n=1 m=1 p=1

Rappelons que le terme χn,m,p correspond à la valeur du contraste χ sur le voxel d’indice (n, m, p) selon la représentation volumique des voxels et les entiers Nx, Ny et Nz représentent les dimensions du volume V discrétisé. On retrouve ici la pénalisation sur les diﬀérences premières à travers l’approximation au premier ordre : χn+1,m,p − χn,m,p pour la direction x par exemple. La fonction ϕ(u) = δ2 + |u|2 est la fonction de pénalisation “`2`1” qui permet de préserver les fortes discontinuités. Dans cette fonction, le paramètre δ > 0 doit être réglé. Dans ce manuscrit, nous avons fait le choix de fixer δ de façon empirique.

Ce type de régularisation semble faire l’unanimité en TMO bien que l’appellation ne soit pas partagée. Certains auteurs introduisent la variation totale comme terme de pénalisa-tion (van den Berg et Kleinman, 1995; Abubakar et van den Berg, 2001) mais ce terme est équivalent à notre expression dans (3.6).

Régularisation additive ou multiplicative ?

Dans la majorité des méthodes d’inversion, les auteurs proposent d’intégrer le terme de régularisation dans le critère soit :

– sous forme additive. Le critère global prend la forme : F = FD + λR

– sous forme multiplicative. Le critère global prend la forme F = FD × R

La régularisation additive est très fréquente en problèmes inverses. Le paramètre de régula-risation λ permet d’ajuster l’équilibre entre l’adéquation aux données FD et l’information a priori de la régularisation R . Le réglage de ce paramètre est important car il influence le ré-sultat de la reconstruction. Choisir λ faible signifie que l’on fait plutôt confiance aux données et au modèle (par exemple, faible bruit sur les données). À l’inverse, choisir λ élevé signifie que l’on attribue plus de crédit à l’information a priori. Le choix de λ est donc fortement lié au niveau de bruit sur les données (Idier, 2001) : plus le niveau de bruit est élevé, plus λ sera élevé. Des techniques ont été développées pour ce réglage (Idier, 2001), on peut citer la méthode de la “courbe en L” par exemple. Ces techniques peuvent être coûteuses et le réglage de λ est alors réalisé empiriquement. Par exemple, Barrière et al. (2011) proposent une procédure expérimentale pour ce réglage : une première étape consiste à choisir un objet de référence suﬃsamment grand et contrasté pour représenter le type d’objets à reconstruire pour l’application donnée, la seconde étape consiste à eﬀectuer un certain nombre de recons-tructions de l’objet référence en faisant varier les paramètres de la régularisation (ici, λ et δ) sur une grille logarithmique. Les valeurs des paramètres sont alors choisies en fonction de la qualité de la solution pour les diﬀérentes reconstructions de l’objet référence.

La régularisation multiplicative est apparue dans les algorithmes d’inversion en TMO depuis quelques années et semble avoir convaincu un grand nombre d’auteurs. L’avantage avancé pour une telle régularisation est le fait qu’il n’y ait pas besoin de paramètre de régularisation, ce qui évite la phase délicate du réglage de celui-ci. Barrière et al. (2011) montrent en revanche que ce type de régularisation présente des limites notamment pour la reconstruction d’objets fortement contrastés. D’autres approches de régularisation ont été étudiées dans la littérature, comme par exemple De Zaeytijd et al. (2007) qui mélange les régularisations additives et multiplicatives.

Dans la suite, nous avons choisi une régularisation additive. Pour le réglage du paramètre de régularisation λ, nous avons choisi de le fixer de façon empirique en suivant la procédure décrite dans Barrière et al. (2011).

Les algorithmes d’optimisation utilisés

Après avoir fait le choix d’un critère à minimiser et d’une régularisation qui convient au pro-blème étudié, il est nécessaire de choisir un algorithme d’optimisation approprié. Comme le critère n’est pas convexe, des minima locaux peuvent être présents. La minimisation d’un tel critère nécessite généralement des méthodes d’optimisation globale, qu’elles soient détermi-nistes ou stochastiques (Pastorino, 2007). Cependant, à cause du grand nombre d’inconnues et de la diﬃculté du problème, ces approches sont trop coûteuses. Malgré la non-linéarité du problème et le risque de tomber dans un minimum local non global, les approches d’optimisa-tion locale sont tout de même robustes et permettent d’obtenir des résultats de reconstruction acceptables. Ce sont les algorithmes les plus fréquents pour minimiser le critère et nous énu-mérerons les principales méthodes ci-dessous. Pour assurer la convergence vers un minimum acceptable, il est nécessaire d’initialiser correctement ces algorithmes afin que l’initialisation appartienne au bassin d’attraction de ce minimum acceptable. Certaines techniques per-mettent d’obtenir une initialisation convenable. On peut citer deux d’entre elles :

– initialisation du contraste avec une solution approchée. Par exemple, la méthode de Born (linéarisation du modèle) permet d’obtenir une première idée de l’objet inconnu. Cette technique est surtout valable dans les conditions de validité de l’approximation de Born. On peut également citer la méthode Linear Sampling Method (LSM) qui peut être utilisée pour déterminer une première approximation de la forme et de la position de l’objet (Catapano et al., 2009) ;

– la technique par saut de fréquence (frequency hopping). L’idée est de reconstruire l’objet à partir de mesures à diﬀérentes fréquences. Les basses fréquences permettent de résoudre un problème faiblement non-linéaire puisque l’objet est de taille électrique faible à ces fréquences. On obtient alors une solution lissée de l’objet inconnu. En reconstruisant l’objet à des fréquences de plus en plus élevées, il est possible d’obtenir une solution de plus en plus fidèle de l’objet. Cette technique est très utilisée, c’est par exemple le cas dans la plupart des méthodes reconstruisant les objets de la base Fresnel. Cette technique nécessite d’avoir fait l’acquisition des mesures à diﬀérentes fréquences en amont de la reconstruction.

Nous allons nous focaliser sur les méthodes d’optimisation locale qui couvrent la majeure partie des méthodes d’inversion en TMO. L’idée générale est d’estimer une direction de descente d et un pas de descente α à partir d’une estimation courante du contraste. Notons à l’itération k, l’itéré courant x(k). L’itéré suivant x(k+1) est défini sous la forme : x(k+1) = x(k) + α(k)d(k), (3.7) de façon à ce que F(x(k+1)) < F(x(k)). L’algorithme est arrêté quand l’itéré est proche d’un minimum local : en pratique, quand la norme du gradient rxF(x(k)) est inférieure à un seuil défini par l’utilisateur. La diﬀérence principale entre les algorithmes itératifs est la façon dont la direction d est calculée. On peut classer les algorithmes en deux catégories : les algorithmes du premier ordre, c’est-à-dire n’utilisant que l’information du gradient pour calculer les directions de descente ; et les algorithmes du second ordre, utilisant le gradient et la matrice Hessienne afin de calculer une direction de descente plus eﬃcace. Les méthodes du premier ordre sont généralement moins coûteuses pour calculer la direction de descente à chaque itération puisqu’il n’est pas nécessaire d’estimer la matrice Hessienne. En revanche, elles requièrent souvent un plus grand nombre d’itérations avant convergence. À l’inverse, le coût de calcul du Hessien pouvant être important, le coût d’une itération pour les méthodes du second ordre est souvent très élevé. Mais la convergence est atteinte en moins d’itérations, notamment quand le critère est proche d’une forme quadratique. Concernant le pas de descente α(k), trouver le pas optimal minimisant le critère selon la direction de descente est trop coûteux pour un critère non-linéaire. En général, le pas de descente est choisi de façon à respecter certaines conditions. Les conditions fortes de Wolfe peuvent par exemple être appliquées (Nocedal et Wright, 1999). Une condition moins stricte, la condition d’Armijo, est souvent utilisée. Plusieurs procédures de recherche de pas ont été développées afin de trouver un pas respectant une de ces conditions. On peut citer la procédure de backtracking fréquemment utilisée.

Revenons aux algorithmes d’optimisation où nous présentons ceux qui sont largement appli-qués à la reconstruction en imagerie micro-ondes. Les méthodes ne seront pas détaillées ici ; pour plus d’informations, le lecteur pourra se référer à Bertsekas (1999).

La méthode du gradient conjugué non-linéaire

La méthode du gradient conjugué non-linéaire est l’extension du réputé gradient conjugué pour des fonctions non-linéaires. On peut exprimer la direction de descente à l’itération k sous la forme : d(k) = −rxF(x(k)) + βkd(k−1), (3.8) où βk est calculé selon une certaine formule. Les formules les plus connues sont celles de Fletcher-Reeves et celle de Polak-Ribière (Bertsekas, 1999).

La méthode de Newton

La méthode de Newton est une méthode du second degré basée sur la minimisation successive d’approximations quadratiques du critère. La direction de descente est donnée par : d(k) = −Hk−1 rxF(x(k)), (3.9) où la matrice Hk, de taille N × N, est la matrice Hessienne. Pour des problèmes de grande taille, l’inversion directe de Hk est impossible ; il est nécessaire d’utiliser des méthodes ité-ratives pour résoudre le système linéaire afin de déterminer la direction de descente d(k). De plus, de par sa grande taille, la matrice Hk ne peut pas être stockée en mémoire directe-ment. Tout comme la matrice Lx dans le modèle direct, seules quelques quantités matricielles exprimant une formulation implicite de Hk doivent être stockées.

La méthode originale de Newton consiste à choisir un pas de descente α(k) = 1. Si la matrice Hessienne n’est pas définie positive, ce choix peut entraîner non pas une direction de descente mais une direction de montée. Pour contourner ce problème, la variante relâchée de la méthode de Newton propose de choisir un pas de descente selon la même démarche que les autres algorithmes (par exemple par backtracking) qui peut être négatif pour assurer la descente de l’algorithme.

La méthode de Gauss-Newton

La méthode de Gauss-Newton (GN) est une méthode utilisable uniquement pour la minimi-sation de critères des moindres carrés. On peut voir cette approche comme une variante de la méthode de Newton dans laquelle la matrice Hessienne est approchée par Hk ≈ J†J, où J est le jacobien du modèle en x(k). On obtient donc la direction de descente suivante : d(k) = −(J†J)−1 rxF(x(k)). (3.10)

La matrice J†J étant semi-définie positive, il est assuré que la direction est une direction de descente. Tout comme la méthode originale de Newton, la méthode de Gauss-Newton propose un pas de descente α(k) = 1. Cela n’assure pas nécessairement que F(x(k+1)) < F(x(k)). On peut alors ajouter une recherche de pas eﬃcace afin de satisfaire aux conditions d’Armijo par exemple. Cet algorithme est très utilisé pour le problème inverse en TMO, on le retrouve notamment dans la littérature récente (Mojabi et LoVetri, 2009; Abubakar et al., 2012a; De Zaeytijd et Franchois, 2009; Ostadrahimi et al., 2013; Estatico et al., 2013).

Tout comme l’algorithme de Newton, il est impossible de stocker la matrice J†J directement pour des problèmes de grande taille. Il est nécessaire d’utiliser d’autres solutions pour calculer la direction de descente. Une solution possible est la formulation implicite du jacobien décrite dans Abubakar et al. (2012a).

La méthode de Levenberg-Marquardt

La méthode de Levenberg-Marquardt (LM) est une amélioration de la méthode de Gauss- Newton qui permet de stabiliser la solution en améliorant le conditionnement de la matrice J†J. La direction de descente est : d(k) = −(J†J + λM)−1 rxF(x(k)). (3.11)

On peut voir la matrice M comme un terme de régularisation quadratique, autrement appelé régularisation de Tikhonov. Le paramètre λ peut aussi être considéré comme un facteur de régularisation. Le choix le plus courant pour M est de prendre la matrice identité.

Il a été prouvé que cette méthode est équivalente à la méthode Distorted Born Iterative Method (DBIM) avec régularisation quadratique (Chew et Wang, 1990; Joachimowicz et al., 1991; Franchois et Pichot, 1997). De nombreux auteurs utilisent le dénommé DBIM (Had-dadin et Ebbini, 1998; Remis et van den Berg, 2000; Yu et al., 2009; Gilmore et al., 2009b;

Scapaticci et al., 2015). De plus, l’algorithme de Gauss-Newton est équivalent à la méthode DBIM pour laquelle aucune régularisation n’est utilisée (Carfantan, 1996; Azghani et al., 2015).

La méthode BFGS

La méthode Broyden-Fletcher-Goldfarb-Shanno (BFGS) est un algorithme de quasi-Newton, c’est-à-dire que l’inverse de la matrice Hessienne est approchée selon une formule donnée, uniquement à partir des gradients aux itérations précédentes. La direction de descente est définie par : d(k) = −Bk rxF(x(k)), (3.12)

où la matrice Bk est une approximation de l’inverse de la matrice Hessienne. Cette matrice est calculée à partir de Bk−1, de la direction précédente et du gradient courant rxF(x(k)). Les diﬀérents algorithmes de quasi-Newton se distinguent par la façon dont est calculée la matrice Bk. L’intérêt de ces méthodes est qu’il n’y a pas de système linéaire à résoudre pour calculer la direction de descente, contrairement à la méthode de Newton, GN ou LM. Quelques auteurs ont utilisé cet algorithme pour la reconstruction en TMO (Franchois et Tijhuis, 2003; Hu et al., 2006; De Zaeytijd et al., 2007).

En revanche, comme les méthodes de Newton, GN et LM, la matrice Bk est de taille N × N et il n’est donc pas envisageable de la stocker directement en mémoire. Pour les problèmes de grande taille, une version permettant de limiter la mémoire nécessaire a été développée : L-BFGS (Limited-Memory BFGS) (Nocedal et Wright, 1999; Byrd et al., 1995). Cet algo-rithme permet de ne stocker qu’un certain nombre de vecteurs afin de représenter de façon implicite la matrice Bk. Dans la suite de nos travaux, nous avons choisi d’utiliser l’algorithme L-BFGS.

Le rapport de stage ou le pfe est un document d’analyse, de synthèse et d’évaluation de votre apprentissage, c’est pour cela clepfe.com propose le téléchargement des modèles complet de projet de fin d’étude, rapport de stage, mémoire, pfe, thèse, pour connaître la méthodologie à avoir et savoir comment construire les parties d’un projet de fin d’étude.

Table des matières

LISTE DES ANNEXES
CHAPITRE 1 INTRODUCTION GÉNÉRALE
1.1 La tomographie micro-ondes (TMO)
1.2 Méthodes de reconstruction en TMO
1.3 Difficultés et objectifs de recherche
1.3.1 Optimiser le nombre de systèmes linéaires
1.3.2 Accélérer la résolution des systèmes linéaires
1.3.3 Étudier l’effet de la troncature et améliorer la convergence
1.4 Plan de la thèse
1.5 Cadre du travail
CHAPITRE 2 LE MODÈLE DIRECT
2.1 Modélisation physique de la diffraction d’ondes électromagnétiques
2.1.1 Équation différentielle régissant la propagation d’une onde électromagnétique
2.1.2 Solution intégrale de l’équation d’onde dans un milieu homogène
2.1.3 Équations intégrales en TMO
2.2 Discrétisation du modèle physique et formulation matricielle
2.2.1 État de l’art des différentes approches de discrétisation des équations électromagnétiques
2.2.2 Discrétisation des équations intégrales par la méthode des moments
2.3 Implémentation du problème direct
2.3.1 Résolution du système linéaire
2.3.2 Multiplication par Lx : comment réduire l’espace mémoire
2.4 Validation du modèle direct
2.4.1 Comparaison à la solution analytique de Mie
2.4.2 Comparaison aux données réelles de l’Institut Fresnel
2.5 Conclusion
CHAPITRE 3 MÉTHODES D’INVERSION EN IMAGERIE MICRO-ONDES
3.1 L’imagerie micro-ondes
3.2 L’inversion vue comme un problème d’optimisation
3.2.1 Formulation bilinéaire
3.2.2 Formulation non-linéaire
3.3 La régularisation
3.3.1 Régularisation “`2`1”
3.3.2 Régularisation additive ou multiplicative ?
3.4 Les algorithmes d’optimisation utilisés
3.5 Calcul du gradient pour le critère non-linéaire
3.6 Conclusion
CHAPITRE 4 ÉTUDE ALGÉBRIQUE DU GRADIENT ET RÉDUCTION DU NOMBRE DE SYSTÈMES LINÉAIRES
4.1 Mise en œuvre du calcul du jacobien
4.1.1 L’état adjoint pour le calcul du jacobien
4.1.2 Équivalence entre l’état adjoint et la formulation matricielle du jacobien
4.2 Mises en œuvre du calcul du gradient
4.2.1 Mise en œuvre directe
4.2.2 Mise en œuvre adjointe du jacobien
4.2.3 Choix de la formulation en fonction du montage d’imagerie
4.3 Théorème de réciprocité de Lorentz
4.4 Procédure finale d’optimisation du nombre de systèmes
4.5 Exemple de configurations présentes dans la littérature
4.5.1 Montage générique d’acquisition
4.5.2 Montage simulé de Zhang et Liu (2015)
4.5.3 Configuration d’acquisition des données Fresnel
4.6 Conclusion
CHAPITRE 5 APPROCHES ITÉRATIVES PAR BLOCS POUR LA RÉSOLUTION DE SYSTÈMES MULTIPLES
5.1 Les approches de résolution de systèmes par blocs
5.2 Adaptation de l’approche Block-BiCGStab pour une application en TMO
5.2.1 L’algorithme Block-BiCGStab
5.2.2 Réglage des paramètres du Block-BiCGStab pour une application en TMO
5.2.3 Résultats justifiant le réglage de l’initialisation
5.3 Amélioration du Block-BiCGStab pour une architecture multi-cœurs
5.3.1 L’approche Partial-Block BiCGStab
5.3.2 Réglages des paramètres de l’algorithme Partial-Block BiCGStab
5.3.3 Effet des réglages du Partial-Block BiCGStab sur la convergence
5.4 Exploitation des approches par blocs dans le processus d’inversion
5.5 Évaluation des coûts de calcul des différents algorithmes pour la résolution de problèmes directs
5.5.1 Influence du contraste
5.5.2 Influence de la fréquence
5.5.3 Influence du nombre de systèmes à résoudre
5.5.4 Effet de l’architecture multi-cœurs
5.6 Évaluation des approches pour la reconstruction
5.6.1 Procédure de reconstruction
5.6.2 Résultats de reconstruction
5.6.3 Comparaison des coûts de reconstruction selon les approches
5.7 Conclusion
CHAPITRE 6 VALIDATION DE NOTRE ALGORITHME D’INVERSION SUR DONNÉES RÉELLES
6.1 Présentation de la base de données de l’Institut Fresnel
6.1.1 La configuration du montage d’acquisition
6.1.2 Les objets de la base
6.2 Mise en œuvre de notre algorithme de reconstruction
6.2.1 Réciprocité de Lorentz et polarisations utilisées
6.2.2 Volume d’étude et discrétisation
6.2.3 Procédure par saut de fréquence
6.3 Résultats de reconstruction
6.4 Comparaison des coûts de calcul des différentes approches
6.5 Conclusion
CHAPITRE 7 VERS UNE GESTION EFFICACE DES EFFETS DE LA TRONCATURE SUR LA COHÉRENCE ENTRE LE CRITÈRE ET LE GRADIENT
7.1 Mise en évidence de l’effet de la troncature sur la convergence de l’algorithme d’optimisation
7.2 Expression du champ total dans un sous-espace
7.2.1 Expression du champ
7.2.2 Choix d’un sous-espace
7.3 Définition du nouveau critère et calcul du gradient
7.3.1 Critère approché
7.3.2 Gradient correspondant au critère approché
7.3.3 Vérification du gradient pour B inversible
7.3.4 Étude de la variation du critère selon la direction du gradient
7.4 Conclusion
CHAPITRE 8 CONCLUSION ET PERSPECTIVES
8.1 Synthèse des travaux
8.1.1 Réduction du nombre de systèmes linéaires
8.1.2 Accélération des résolutions de systèmes
8.1.3 Validation de l’algorithme de reconstruction sur des données expérimentales
8.1.4 Cohérence entre le critère et le gradient
8.2 Limitations des solutions proposées
8.3 Améliorations futures et pistes de recherche
8.3.1 Étude de la formulation cohérente du gradient
8.3.2 Études complémentaires de l’approche Partial-Block BiCGStab
8.3.3 Approches de parallélisation de l’algorithme Block-BiCGStab
8.3.4 Comparer notre algorithme de reconstruction à d’autres algorithmes fréquents en TMO
8.3.5 Élargir le champ des applications
RÉFÉRENCES
ANNEXES