## II.2 Méthode des champs virtuels

La méthode des champs virtuels a été développée pour résoudre des problèmes de mécanique du solide, notamment l'identification des paramètres pilotant les lois de comportement des milieux [10]. Afin d'introduire ces paramètres, le principe des travaux virtuels (PTV) est réécrit en remplaçant les termes du tenseur des contraintes $\sigma$ dans (3.5) par leur expression suivant une loi de comportement correspondant au milieu étudié. Le PTV étant valable pour tous milieux, n'importe quelle loi de comportement peut être utilisée. Les paramètres pilotant cette loi seront ensuite estimés grâce à une sélection judicieuse de champs virtuels.

Les champs virtuels n'ont pas de réel sens physique, mais doivent être vus comme un ou- til mathématique permettant d'estimer les paramètres recherchés. L'intérêt de la méthode des champs virtuels repose donc sur la flexibilité de ces champs, qui peuvent être déterminer judi- cieusement afin d'être le mieux adaptés au problème posé. De même, il est possible de paramétrer les champs virtuels dans le but d'éliminer des termes inconnus présents dans l'équation (3.5). Par exemple, en ERM dynamique, où la force de surface est inconnue, Connesson et al. [12] ont proposé de définir $u^* = 0$ aux frontières dans le but de supprimer ce terme. Par conséquent, le choix des champs virtuels peut permettre d'obtenir un problème relativement simple où les seules inconnues sont les paramètres contrôlant la loi de comportement choisie.

La MCV sera décrite plus en détail dans les sections suivantes. Avant cela, certaines remarques importantes sur l'utilisation des champs virtuels doivent être faites.

1. Premièrement, il est primordial de bien différencier les champs de déplacement et de dé- formation réels/estimés de ceux dits virtuels, malgré leur appellation similaire. Les champs virtuels sont équivalents à des fonctions de tests. Ce sont des fonctions arbitraires, pouvant être totalement différentes des champs réels. De fait, les champs virtuels utilisés par la MCV n'ont pas de lien physique avec les champs réels. Néanmoins, ces champs étant assimilés à des déplacements, la relation entre les déplacements virtuels et les déformations virtuelles est définie par l'équation (3.4). On parlera souvent par la suite de paire de champs virtuels $(u^*, \varepsilon^*)$, abrégée CVF.

2. Deuxièmement, le PTV est valable pour tous déplacements virtuels vérifiant la condition d'admissibilité cinématique [10]. Cette condition signifie que les champs virtuels doivent être continus et respecter les conditions aux limites du déplacement réel. En pratique, il existe une infinité de possibilités.

3. Troisièmement, il est nécessaire d'utiliser autant de CVF qu'il y a d'inconnues dans le pro- blème traité. Par exemple, pour un milieu homogène, linéaire, élastique et isotrope, le pro- blème sera constitué de deux inconnues (module de Young et coefficient de Poisson). Il faudra donc au minimum utiliser deux CVF distincts.

## III Caractérisation de milieux soumis à un chargement quasi-statique

La méthode des champs virtuels décrite ci-dessus est très générale et peut être adaptée à un grand nombre de problèmes, complexifiant plus ou moins son application. En effet, cette méthode offre une grande latitude à l'utilisateur quant au choix de la loi de comportement et des champs virtuels utilisés. Généralement, certains termes du problème, comme la force de volume b, peuvent être négligés. De plus, dans un problème quasi-statique, l'accélération est considérée comme étant suffisamment faible pour ne pas intervenir dans l'équation (3.5). Ainsi, l'écriture du principe des travaux virtuels peut être simplifiée, telle que :
$$ \int_V \sigma : \varepsilon^* dv = \int_S T \cdot u^* ds \qquad (3.7) $$

Pour finir, et comme cela a été indiqué précédemment, lors de mesures expérimentales réalisées en élastographie ultrasonore quasi-statique, les données ultrasonores radiofréquences (RF) acquises au cours du temps sont généralement des images (et non des volumes). Seul le déplacement (et donc la déformation) dans le plan d'imagerie peut alors être estimé, et la composante en élévation, appelée "hors plan", reste inaccessible. Afin de pallier le manque de données 3D, le principe des travaux virtuels peut être écrit sous les hypothèses de contraintes planes ou de déformations planes. Les deux cas sont détaillés ci-dessous.

### III.1 Réduction de dimension du problème

Dans ce travail, nous nous intéressons à l'analyse des propriétés mécaniques de milieux élastiques, linéaires et isotropes. La suite de ce chapitre sera donc présentée en considérant uniquement ce type de milieu. Néanmoins, le développement de la méthode des champs virtuels reste le même quelle que soit la loi de comportement choisie. Pour rappel, la loi de comportement d'un matériau élastique, linéaire et isotrope, écrite selon le module de Young E et le coefficient de Poisson $\nu$ est :

$$ \sigma = \frac{E}{1 + \nu} \left( \varepsilon + \frac{\nu}{1 - 2\nu} \mathrm{Tr}(\varepsilon I) \right) \qquad (3.8) $$

#### III.1.1 Hypothèse des contraintes planes

L'état de contraintes planes se produit dans un milieu de faible épaisseur (par exemple une plaque mince, fig. 3.1) assujetti uniquement à un chargement dans son plan [13]. Introduisons un repère orthonormé ($O, \vec{e}_1, \vec{e}_2, \vec{e}_3$). Avec l'illustration de la figure 3.1, la contrainte ne dépendra pas de $x_3$. En d'autres termes, le tenseur des contraintes sera de la forme :

$$ \sigma = \begin{bmatrix} \sigma_{11}(x_1, x_2) & \sigma_{12}(x_1, x_2) & 0 \\ \sigma_{21}(x_1, x_2) & \sigma_{22}(x_1, x_2) & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad \text{ou} \quad \sigma = \begin{pmatrix} \sigma_1(x_1, x_2) \\ \sigma_2(x_1, x_2) \\ 0 \\ 0 \\ 0 \\ \sigma_6(x_1, x_2) \end{pmatrix} \qquad (3.9) $$

<figure><img src="image_29.png" /><figcaption>FIGURE 3.1 – Milieu quelconque de faible épaisseur dans le cas de l'hypothèse des contraintes planes.</figcaption></figure>
Le tenseur des déformations se déduit du tenseur des contraintes (3.9) et de la loi de comportement (3.8), soit

$$ \varepsilon = \begin{pmatrix} \varepsilon_1(x_1, x_2) \\ \varepsilon_2(x_1, x_2) \\ \varepsilon_3(x_1, x_2) \\ 0 \\ 0 \\ \varepsilon_6(x_1, x_2) \end{pmatrix}, \text{ avec } \varepsilon_3 = \frac{-\nu}{1-\nu}(\varepsilon_1 + \varepsilon_2). \qquad (3.10) $$

Dans le plan considéré, il est alors possible d'écrire la relation contrainte-déformation suivante :

$$ \underbrace{\begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_6 \end{pmatrix}}_{\sigma} = \underbrace{\frac{E}{1-\nu^2} \begin{bmatrix} 1 & \nu & 0 \\ \nu & 1 & 0 \\ 0 & 0 & 1-\nu \end{bmatrix}}_{Q} \underbrace{\begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_6 \end{pmatrix}}_{\varepsilon} \qquad (3.11) $$

### III.1.2 Hypothèse des déformations planes

L'état de déformations planes se produit dans un milieu où l'une des dimensions est grande face aux deux autres (fig. 3.2) et le long de laquelle il est possible de négliger la déformation. L'hypothèse des déformations planes suppose que [13] :

1. l'une des trois composantes du déplacement, celle le long du grand axe du milieu, est constante, soit dans notre exemple $u_3 = \text{cst}$.

2. les deux autres composantes du déplacement ne dépendent pas de cette troisième dimension, soit ici $u_1(x_1, x_2)$ et $u_2(x_1, x_2)$.

<figure><img src="image_30.png" /><figcaption>FIGURE 3.2 - Milieu quelconque, dont l'épaisseur (le long de l'axe ($Ox_3$)) est grande face aux autres dimensions, dans le cas de l'hypothèse des déformations planes.</figcaption></figure>

En d'autres termes, le tenseur des déformations sera de la forme :

$$ \varepsilon = \begin{bmatrix} \varepsilon_{11}(x_1, x_2) & \varepsilon_{12}(x_1, x_2) & 0 \\ \varepsilon_{21}(x_1, x_2) & \varepsilon_{22}(x_1, x_2) & 0 \\ 0 & 0 & 0 \end{bmatrix} \quad \text{ou} \quad \varepsilon = \begin{pmatrix} \varepsilon_1(x_1, x_2) \\ \varepsilon_2(x_1, x_2) \\ 0 \\ 0 \\ 0 \\ \varepsilon_6(x_1, x_2) \end{pmatrix} \qquad (3.12) $$
Le tenseur des contraintes se déduit du tenseur des déformations (3.12) et de la loi de comportement (3.8), soit

$$ \boldsymbol{\sigma} = \begin{pmatrix} \sigma_1(x_1, x_2) \\ \sigma_2(x_1, x_2) \\ \sigma_3(x_1, x_2) \\ 0 \\ 0 \\ \sigma_6(x_1, x_2) \end{pmatrix}, \quad \text{avec} \quad \sigma_3 = \nu(\sigma_1 + \sigma_2). \tag{3.13} $$

Dans le plan considéré, il est alors possible d'écrire la relation contrainte-déformation suivante :

$$ \underbrace{\begin{pmatrix} \sigma_1 \\ \sigma_2 \\ \sigma_6 \end{pmatrix}}_{\sigma} = \underbrace{\frac{E}{(1+\nu)(1-2\nu)}}_{Q} \underbrace{\begin{bmatrix} 1-\nu & \nu & 0 \\ \nu & 1-\nu & 0 \\ 0 & 0 & 1-2\nu \end{bmatrix}}_{\varepsilon} \underbrace{\begin{pmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_6 \end{pmatrix}}_{\varepsilon} \tag{3.14} $$

## III.2 Écriture de la MCV en 2D

Dans ce paragraphe, nous allons plus particulièrement nous intéresser à l'écriture de la méthode des champs virtuels dans le cas des déformations et des contraintes planes. Après simplification par l'épaisseur, le principe des travaux virtuels en 2D s'écrit :

$$ \int_S \boldsymbol{\sigma} : \boldsymbol{\varepsilon}^* ds = \int_L \boldsymbol{T} \cdot \boldsymbol{u}^* dl \tag{3.15} $$

Par la suite, la forme générale de la relation contrainte-déformation ($\boldsymbol{\sigma} = \boldsymbol{Q}\boldsymbol{\varepsilon}$ avec $\boldsymbol{Q}$ contenant l'ensemble des paramètres contrôlant la loi de comportement, soit ici les inconnues du problème) sera utilisée. Il conviendra que :

$$ \boldsymbol{Q} = \begin{bmatrix} Q_1 & Q_2 & 0 \\ Q_2 & Q_1 & 0 \\ 0 & 0 & Q_6 \end{bmatrix} \quad \text{avec} \quad Q_6 = Q_1 - Q_2 \tag{3.16} $$

$Q_1$, $Q_2$ et $Q_6$ s'obtiennent directement par identification avec (3.11) pour le cas de l'hypothèse des contraintes planes et (3.14) pour celui des déformations planes. Les valeurs des composantes de la matrice $\boldsymbol{Q}$ dépendent donc de la loi de comportement choisie et de l'hypothèse utilisée. Cela permet de simplifier la manipulation de la MCV, car ces choix pourront être faits sans modifier l'écriture du problème.

Le PTV (3.15) peut alors être développé tel que :

$$ \int_S (\sigma_1 \varepsilon_1^* + \sigma_2 \varepsilon_2^* + 2\sigma_6 \varepsilon_6^*) ds = \int_L T_i \cdot u_i^* dl \tag{3.17} $$

$$ \int_S (Q_1 (\varepsilon_1 \varepsilon_1^* + \varepsilon_2 \varepsilon_2^* + 2\varepsilon_6 \varepsilon_6^*) + Q_2 (\varepsilon_2 \varepsilon_1^* + \varepsilon_1 \varepsilon_2^* - 2\varepsilon_6 \varepsilon_6^*)) ds = \int_L T_i \cdot u_i^* dl \tag{3.18} $$

A partir de ces équations, plusieurs stratégies pour construire les champs virtuels peuvent être mises en place.
# IV Construction des champs virtuels

La construction des champs virtuels est une étape clé de la MCV. Théoriquement, peu de contraintes existent sur ces champs, bien qu'ils doivent être continus et cinématiquement admis- sibles, laissant une infinité de possibilités. En pratique, la capacité de la méthode à estimer les paramètres recherchés va dépendre des champs virtuels sélectionnés. En effet, les champs virtuels utilisés doivent mener à l'écriture d'un système d'équations indépendantes afin de pouvoir esti- mer chacune des inconnues du problème. Plusieurs approches ont été développées pour construire ces champs virtuels. Elles sont présentées ci-dessous.

## IV.1 Champs virtuels particuliers

Une des premières approches proposées consiste à construire manuellement les champs virt- uels, alors appelés "champs virtuels particuliers" [9, 14]. Ces champs sont généralement définis par des expressions simples conduisant à mettre certains termes à zéro. Cela permet d'assurer l'indépendance des équations.

Un exemple est proposé ci-dessous, illustrant la MCV et l'utilisation de champs virtuels particuliers pour un cas simple.

**Exemple pour un milieu homogène :**

Considérons un milieu homogène, élastique, linéaire et isotrope de dimensions $L_x \times L_y \times e$ m$^3$. Le milieu est en appui en $x_2 = 0$ et soumis à une compression uniaxiale $\mathbf{F} = (0, F)^t$, appliquée en $x_2 = L_y$ (fig. 3.3). Ici, les inconnues sont le module de Young $E$ et le coefficient de Poisson $\nu$.

<figure><img src="image_31.png" /><figcaption>FIGURE 3.3 – Illustration du problème</figcaption></figure>

Le milieu étant homogène, les propriétés mécaniques sont identiques en tout point. Ainsi, seule une valeur de $E$ et de $\nu$ est recherchée, ce qui nécessitera l'utilisation de deux CVF indépendants : $(\mathbf{u}^{*1}, \boldsymbol{\varepsilon}^{*1})$ et $(\mathbf{u}^{*2}, \boldsymbol{\varepsilon}^{*2})$. Dans ce cas, l'équation (3.18) s'écrit :

$$Q_1 \int_S (\varepsilon_1 \varepsilon_1^* + \varepsilon_2 \varepsilon_2^* + 2\varepsilon_6 \varepsilon_6^*) ds + Q_2 \int_S (\varepsilon_2 \varepsilon_1^* + \varepsilon_1 \varepsilon_2^* - 2\varepsilon_6 \varepsilon_6^*) ds = \int_L T_i \cdot u_i^* dl \quad (3.19)$$

En utilisant les deux CVF, deux équations sont obtenues, réécrites sous la forme d'un système linéaire $\mathbf{A}\mathbf{Q}=\mathbf{B}$, avec :

$$\mathbf{A} = \begin{bmatrix} \int_S (\varepsilon_1 \varepsilon_1^{*1} + \varepsilon_2 \varepsilon_2^{*1} + 2\varepsilon_6 \varepsilon_6^{*1}) \, ds & \int_S (\varepsilon_2 \varepsilon_1^{*1} + \varepsilon_1 \varepsilon_2^{*1} - 2\varepsilon_6 \varepsilon_6^{*1}) \, ds \\ \int_S (\varepsilon_1 \varepsilon_1^{*2} + \varepsilon_2 \varepsilon_2^{*2} + 2\varepsilon_6 \varepsilon_6^{*2}) \, ds & \int_S (\varepsilon_2 \varepsilon_1^{*2} + \varepsilon_1 \varepsilon_2^{*2} - 2\varepsilon_6 \varepsilon_6^{*2}) \, ds \end{bmatrix} \quad (3.20a)$$

$$\mathbf{Q} = \begin{pmatrix} Q_1 \\ Q_2 \end{pmatrix} \qquad (3.20b)$$

$$\mathbf{B} = \begin{pmatrix} \int_L \frac{F}{eL_x} u_2^{*1}(x_1, L_y) dl \\ \int_L \frac{F}{eL_x} u_2^{*2}(x_1, L_y) dl \end{pmatrix} \quad (3.20c)$$
où chaque ligne correspond à un CVF indépendant et cinématiquement admissible. Notez que dans **B**, seule la composante axiale du déplacement virtuel apparaît, puisque la force est verticale.
Un choix simple pour définir ces champs peut être le suivant :

$$
\text{CVF 1 :} \quad \left\{ \begin{array}{l} u_1^{*1}(x_1, x_2) = x_1 \\ u_2^{*1}(x_1, x_2) = 0 \end{array} \right. \qquad \left\{ \begin{array}{l} \varepsilon_1^{*1}(x_1, x_2) = 1 \\ \varepsilon_2^{*1}(x_1, x_2) = 0 \\ \varepsilon_6^{*1}(x_1, x_2) = 0 \end{array} \right. \qquad (3.21)
$$

$$
\text{CVF 2 :} \quad \left\{ \begin{array}{l} u_1^{*2}(x_1, x_2) = 0 \\ u_2^{*2}(x_1, x_2) = x_2 \end{array} \right. \quad \left\{ \begin{array}{l} \varepsilon_1^{*2}(x_1, x_2) = 0 \\ \varepsilon_2^{*2}(x_1, x_2) = 1 \\ \varepsilon_6^{*2}(x_1, x_2) = 0 \end{array} \right. \qquad (3.22)
$$

En injectant ces CVF dans (3.20), cela donne :

$$
\underbrace{\left[ \begin{array}{cc} \int_S \varepsilon_1 ds & \int_S \varepsilon_2 ds \\ \int_S \varepsilon_2 ds & \int_S \varepsilon_1 ds \end{array} \right]}_{\mathbf{A}} \underbrace{\left( \begin{array}{c} Q_1 \\ Q_2 \end{array} \right)}_{\mathbf{Q}} = \underbrace{\left( \begin{array}{c} 0 \\ \frac{FL_y}{e} \end{array} \right)}_{\mathbf{B}} \qquad (3.23)
$$

La matrice **Q** est alors calculée par inversion du système (**Q** = **A**⁻¹**B**). Les propriétés mécaniques *E* et *ν* sont ensuite directement déduites à partir de leur relation avec *Q*₁ et *Q*₂ (cf. III.1).

L'utilisation de champs virtuels particuliers permet d'estimer directement et rapidement les paramètres de la loi de comportement choisie. Cependant, face à l'infinité de champs virtuels possibles, ceux sélectionnés manuellement ne permettent pas toujours une estimation optimale de ces paramètres. Plus important encore, l'exemple présenté ci-dessus est très simple et nécessite de définir seulement deux champs virtuels, puisque deux valeurs sont à déterminer. Lorsque le nombre d'inconnues augmente, définir manuellement plusieurs champs virtuels cinématiquement admissibles est complexe. De plus, ils peuvent mener à des systèmes faiblement indépendants rendant le problème plus difficile à résoudre. Ce type d'approche n'est donc pas adaptée aux problèmes constitués d'un grand nombre d'inconnues, comme cela est par exemple le cas pour les milieux hétérogènes.

IV.2 Champs virtuels spécifiques

Face aux limites des champs virtuels particuliers, une seconde approche a été développée afin de construire automatiquement des champs virtuels adaptés au problème posé [10]. L'idée est de déterminer des champs virtuels qui permettent d'extraire directement les paramètres recherchés en rendant la matrice **A** égale à la matrice identité (3.24). Ces champs sont appelés "champs virtuels spécifiques".

$$
\underbrace{\begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix}}_{\mathbf{A}} \underbrace{\begin{pmatrix} Q_1 \\ Q_2 \end{pmatrix}}_{\mathbf{Q}} = \underbrace{\begin{pmatrix} \int_L T_i u_i^{*1} dl \\ \int_L T_i u_i^{*2} dl \end{pmatrix}}_{\mathbf{B}} \qquad (3.24)
$$

Construire les CVF tel que **A** = **I** (matrice identité) permet de s'affranchir de l'inversion de la matrice **A** en obtenant directement **Q**=**B**. Cela permet également d'assurer directement
l'indépendance des équations.

Les champs virtuels spécifiques sont construits à partir de fonctions mathématiques définissant le déplacement virtuel sur l'ensemble du milieu. Différentes fonctions peuvent être utilisées (fonctions polynomiales [15], fonctions par morceaux [16] ...) permettant des écritures différentes de CVF. Néanmoins, la démarche reste la même quelle que soit la forme des CVF. Pour illustrer la MCV avec des champs virtuels spécifiques, un exemple est proposé ci-après pour un cas simple de milieu homogène.

### Exemple pour un milieu homogène :

Reprenons l'exemple précédent du milieu homogène, élastique, linéaire et isotrope, en appui en $x_2 = 0$ et soumis à une compression uniaxiale $F$ appliquée en $x_2 = L_y$ (fig. 3.3). Ici encore, l'utilisation de deux CVF indépendants est nécessaire afin d'estimer le module de Young $E$ et le coefficient de Poisson $\nu$. Pour cet exemple, les champs virtuels sont construits grâce à des fonctions polynomiales. Chaque CVF est donc de la forme suivante :

$$ egin{cases} u_1^*(x_1, x_2) = \sum_{i=0}^{m} \sum_{j=0}^{n} \alpha_{ij} \left(\frac{x_1}{L_x}\right)^i \left(\frac{x_2}{L_y}\right)^j \\ u_2^*(x_1, x_2) = \sum_{i=0}^{p} \sum_{j=0}^{q} \beta_{ij} \left(\frac{x_1}{L_x}\right)^i \left(\frac{x_2}{L_y}\right)^j \\ \varepsilon_1^*(x_1, x_2) = \sum_{i=0}^{m} \sum_{j=0}^{n} \frac{i}{L_x} \alpha_{ij} \left(\frac{x_1}{L_x}\right)^{i-1} \left(\frac{x_2}{L_y}\right)^j \\ \varepsilon_2^*(x_1, x_2) = \sum_{i=0}^{p} \sum_{j=0}^{q} \frac{j}{L_y} \beta_{ij} \left(\frac{x_1}{L_x}\right)^i \left(\frac{x_2}{L_y}\right)^{j-1} \\ \varepsilon_6^*(x_1, x_2) = \frac{1}{2} \left[ \sum_{i=0}^{m} \sum_{j=0}^{n} \frac{j}{L_y} \alpha_{ij} \left(\frac{x_1}{L_x}\right)^i \left(\frac{x_2}{L_y}\right)^{j-1} + \sum_{i=0}^{p} \sum_{j=0}^{q} \frac{i}{L_x} \beta_{ij} \left(\frac{x_1}{L_x}\right)^{i-1} \left(\frac{x_2}{L_y}\right)^j \right] \end{cases} \quad (3.25) $$

où $m, n, p, q$ définissent le nombre maximal de monômes et déterminent la taille du problème à résoudre. Ici, les fonctions polynomiales permettent de satisfaire facilement la condition d'admissibilité cinématique. Cependant, dans le cas contraire, une fonction $G$ venant multiplier les déplacements peut être introduite pour satisfaire cette condition.

L'utilisation de champs virtuels spécifiques modifie donc l'enjeu du problème, où il n'est plus tant de trouver $Q = A^{-1}B$, mais de construire les CVF de manière à obtenir $A = I$ (assurant de fait $Q = B$), soit :

$$ \begin{bmatrix} \int_S (\varepsilon_1 \varepsilon_1^{*1} + \varepsilon_2 \varepsilon_2^{*1} + 2\varepsilon_6 \varepsilon_6^{*1}) \, ds & \int_S (\varepsilon_2 \varepsilon_1^{*1} + \varepsilon_1 \varepsilon_2^{*1} - 2\varepsilon_6 \varepsilon_6^{*1}) \, ds \\ \int_S (\varepsilon_1 \varepsilon_1^{*2} + \varepsilon_2 \varepsilon_2^{*2} + 2\varepsilon_6 \varepsilon_6^{*2}) \, ds & \int_S (\varepsilon_2 \varepsilon_1^{*2} + \varepsilon_1 \varepsilon_2^{*2} - 2\varepsilon_6 \varepsilon_6^{*2}) \, ds \end{bmatrix} = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix} \quad (3.26) $$

Le problème peut donc être réécrit en posant $A = DY$, où $Y$ est un vecteur contenant les coefficients $\alpha_{ij}$ et $\beta_{ij}$ déterminés de manière à satisfaire la condition (3.26) [15, 17, 18]. De même que pour les champs virtuels particuliers, une valeur de $\nu$ et $E$ est finalement estimée pour l'ensemble du milieu, suivant leur relation avec $Q_1$ et $Q_2$ (cf. III.1).

Les champs virtuels spécifiques permettent donc de construire automatiquement des CVF capables d'estimer les propriétés mécaniques du milieu étudié. Cependant, un grand nombre de
champs virtuels reste possible. Ainsi, des conditions supplémentaires peuvent donc être ajoutées par exemple, pour sélectionner des champs virtuels permettant de minimiser l'effet du bruit présent dans les données [6, 18, 19]. Pour ce faire, les champs virtuels montrant une meilleure stabilité face à l'ajout d'un bruit connu sont recherchés. On parle alors de champs virtuels optimisés.

## IV.3 Application aux milieux hétérogènes

Jusqu'à présent, nous avons présenté la méthode des champs virtuels pour des milieux homogènes, notamment parce que ces milieux constituent l'une des principales applications de la MCV en mécanique du solide. De plus, illustrer cette méthode avec des milieux homogènes permet de l'appréhender plus facilement. Cependant, la MCV peut aussi être utilisée pour reconstruire les propriétés mécaniques de milieux hétérogènes, même si les développements réalisés restent encore peu nombreux. Les approches proposées reprennent des éléments présentés précédemment et sont plus particulièrement détaillées ci-après [10]. Afin d'entrevoir la manière dont sont étudiés ces milieux, une illustration est proposée sur un cas particulier.

Considérons un milieu hétérogène, élastique, linéaire et isotrope, en appui en $x_2 = 0$ et soumis à une compression uniaxiale $F$ appliquée en $x_2 = L_y$ (fig. 3.4). Pour simplifier nos propos, le milieu est constitué de deux régions $S$ et $S'$ aux propriétés mécaniques différentes.

<figure><img src="image_32.png" /><figcaption>FIGURE 3.4 – Schéma d'un milieu hétérogène composé de deux régions homogènes</figcaption></figure>

Dès lors que ces régions sont exactement localisées, l'équation (3.18) peut être réécrite sous la forme de plusieurs intégrales correspondant aux différentes régions homogènes. Dans le cas présent, le milieu contient deux régions homogènes, aux inconnues $Q_1$ et $Q_2$ pour $S$ et $Q'_1$ et $Q'_2$ pour $S'$. On obtient ainsi :

$$ \int_S (Q_1 (\varepsilon_1 \varepsilon_1^* + \varepsilon_2 \varepsilon_2^* + 2\varepsilon_6 \varepsilon_6^*) + Q_2 (\varepsilon_2 \varepsilon_1^* + \varepsilon_1 \varepsilon_2^* - 2\varepsilon_6 \varepsilon_6^*)) ds + \\ \int_{S'} (Q'_1 (\varepsilon_1 \varepsilon_1^* + \varepsilon_2 \varepsilon_2^* + 2\varepsilon_6 \varepsilon_6^*) + Q'_2 (\varepsilon_2 \varepsilon_1^* + \varepsilon_1 \varepsilon_2^* - 2\varepsilon_6 \varepsilon_6^*)) ds = \int_L T_i u_i^* dl \quad (3.27) $$

Cette équation peut être résolue comme pour les milieux homogènes précédents, aussi bien en utilisant des champs virtuels particuliers que spécifiques. Dans l'exemple présent, quatre CVF indépendants seront nécessaires pour résoudre l'équation (3.27). Cependant, le nombre de CVF/d'inconnues sera d'autant plus grand que le nombre de régions homogènes sera élevé, rendant les champs virtuels particuliers moins adaptés. Enfin, en pratique, cette méthode trouve rapidement sa limite, car elle nécessite de connaitre au préalable les régions homogènes.

Une seconde approche consiste à définir la distribution spatiale des inconnues contenues dans la matrice $Q$ par une fonction dépendant de $x_1$ et $x_2$. Ainsi, les propriétés mécaniques peuvent être estimées en chaque point du milieu, et non plus en des régions homogènes. A titre d'exemple, $Q$, tout comme les champs virtuels, peut être défini à l'aide de fonctions polynomiales, soit pour la première composante $Q_1(x_1, x_2) = \sum_{i=0}^{m} \sum_{j=0}^{n} \alpha_{ij} \left(\frac{x_1}{L_x}\right)^i \left(\frac{x_2}{L_y}\right)^j$. Le nombre de CVF
devant être utilisés dépend alors du nombre maximal de monômes. De plus, les valeurs choisies pour m et n vont conditionner la capacité de la méthode à détecter les variations spatiales des paramètres estimés. Ainsi, Kim et al. [20] ont développé ce type d'approche pour détecter des dommages de délaminage dans des milieux composites. Il est en effet attendu que ces zones de dommages soient traduites par une variation des propriétés mécaniques, conduisant donc à un milieu mécaniquement hétérogène. Dans leur étude, l'ensemble des composantes de $Q$ sont définies par la même fonction polynomiale $P$, soit $Q = Q_0(1 - P)$, limitant le nombre de CVF utilisés [10, 20, 21]. Cette fonction permet alors de détecter des variations par rapport à une valeur de référence $Q_0$ correspondant aux propriétés mécaniques des zones intactes dans le milieu. Plus précisément, ils se sont intéressés à l'impact du degré de ces fonctions servant à caractériser le milieu et ont montré que l'utilisation de fonctions polynomiales de faibles degrés ne permettaient pas de détecter correctement les variations mécaniques et que la zone endommagée pouvait même être manquée lorsqu'elle est de petite taille. Cependant, un degré élevé implique un nombre d'inconnues plus important et peut mener à des instabilités (fig. 3.5).

<figure><img src="image_33.png" /><figcaption>FIGURE 3.5 – Illustration des résultats obtenus pas Kim et al. [20] en fonction du degré de $P$. Un degré 4 ne permet pas de localiser la zone endommagée avec précision, et à l'inverse, le résultat devient instable pour le degré 12. Les dommages sont localisés de manière satisfaisante pour les degrés 6 à 10.</figcaption></figure>

Au vu des résultats présentés ci-dessus, nous pouvons remarquer que si ce type d'approches permet de détecter des zones ayant des propriétés mécaniques différentes, elles sont cependant peu sensibles aux variations rapides au sein du milieu. Plus récemment, Nguyen et al. [22] ont proposé une approche basée sur l'utilisation de séries de Fourier, et dont les résultats en simulation ont montré une nette sensibilité aux variations rapides.

Enfin, une dernière approche consiste à définir les déplacements virtuels par morceaux, et non par une même fonction sur l'ensemble du milieu, ce qui autorise une plus grande flexibilité [10]. Dans ce cas, une subdivision du domaine en autant de sous-domaines que souhaités peut être opérée, et l'estimation faite pour chaque sous-domaine permettra de reconstruire des cartographies des propriétés mécaniques du milieu. Cette approche permet d'imposer facilement des conditions/contraintes à des nœuds spécifiques, simplifiant la construction de champs virtuels cinématiquement admissibles. Par exemple, si une condition aux limites impose le déplacement virtuel nul à une frontière, alors le déplacement des nœuds situés sur cette frontière pourra être imposé comme étant nul, tout en laissant ceux alentour libres. Cependant, ce type d'approche a aussi l'inconvénient d'augmenter considérablement le nombre de CVF utilisés, car elle nécessite d'estimer autant d'inconnues qu'il y a de sous-domaines.
# V Méthode développée pour l'élastographie ultrasonore quasi-statique

Comme nous l'avons vu précédemment, la méthode des champs virtuels a principalement été utilisée pour résoudre des problèmes de mécanique du solide. Parmi ces applications, beaucoup de travaux s'attachent à estimer les propriétés mécaniques de milieux homogènes [10, 14, 23-27], mais peu de développements ont été réalisés pour étudier les milieux hétérogènes, malgré les quelques approches présentées ci-dessus [19, 20, 22]. À notre connaissance, de premiers travaux impliquant l'utilisation de la MCV en élastographie sont présents dans la littérature [7, 8, 12, 28-30], avec là encore, essentiellement des mesures réalisées sur des milieux homogènes et quelques résultats sur des milieux hétérogènes. Dans ces travaux de thèse, nous nous intéressons à la reconstruction des propriétés mécaniques des tissus biologiques en élastographie ultrasonore quasi-statique. En raison de la complexité du problème, certaines hypothèses et certains choix, ont dû être faits. Ces hypothèses/choix seront d'abord discutés, puis la méthode développée sera présentée.

## V.1 Présentation du problème : choix et hypothèses

L'élastographie est une technique d'imagerie, destinée à l'examen des tissus biologiques. Ces derniers sont extrêmement complexes et des hypothèses simplificatrices vont donc être utilisées. Pour ces premiers travaux, les milieux étudiés seront considérés comme étant élastiques, linéaires et isotropes (cf. Chap. 1). Ce modèle est certes très simple, mais il permet de proposer et d'évaluer les premiers développements d'une méthode de reconstruction des propriétés mécaniques d'un milieu. Il est de plus important de souligner que ce modèle se retrouve dans les méthodes d'élastographie par ondes de cisaillement implémentées sur des systèmes cliniques et actuellement utilisées pour le diagnostic. Commencer nos travaux avec ce type de milieux n'est donc pas dénué d'intérêt.

Un milieu élastique, linéaire et isotrope est localement décrit par deux constantes, le module d'Young E et le coefficient de Poisson ν. Les tissus biologiques auxquels nous nous intéressons sont généralement considérés comme quasi-incompressibles, car composés en majorité d'eau. Dans ces travaux, nous ferons donc l'hypothèse supplémentaire que le coefficient de Poisson est connu, constant sur tout le milieu, et égal à 0.49. Ainsi, dans le problème présent, seule la reconstruction de la distribution spatiale du module d'Young d'un milieu à partir de la connaissance des champs de déplacement/déformation et de la force appliquée est considérée.

Enfin, l'échographie est une modalité d'imagerie qui fournit majoritairement des images (et non des volumes) des milieux au cours du temps, conduisant en élastographie à une estimation 2D du déplacement et de la déformation (cf. chap. 2). Les déplacements dits "hors plan" (ou en élévation), c.-à-d. perpendiculaires au plan d'imagerie, sont inaccessibles. Dans ces travaux, la méthode sera donc développée en 2D. Pour cela, comme nous l'avons vu dans la section III, le PTV sera écrit sous les hypothèses de contraintes planes et de déformations planes afin de pallier le manque de données 3D. Ces hypothèses sont couramment utilisées en élastographie, bien qu'en pratique, aucune d'entre elles ne permettent de modéliser parfaitement le comportement des tissus biologiques [31]. La loi de comportement sous les hypothèses de contraintes planes (3.11) ou de déformations planes (3.14), écrite précédemment $\sigma = Q\epsilon$, est réécrite :

$$ \sigma = EC\epsilon, \quad C = \begin{bmatrix} C_1 & C_2 & 0 \\ C_2 & C_1 & 0 \\ 0 & 0 & C_6 \end{bmatrix} \qquad (3.28) $$
avec *C* une matrice contenant les termes liés à *v*. Nous précisons dès à présent que les deux hypothèses de contraintes planes et de déformations planes seront utilisées dans ce travail, modifiant simplement l'écriture des coefficients de la matrice *C* (3.28).

## V.2 Description de la méthode développée

### V.2.1 Type de champs virtuels utilisé

Comme nous l'avons montré précédemment, diverses approches ont été décrites pour choisir les champs virtuels. En élastographie ultrasonore quasi-statique, les milieux imagés sont des milieux hétérogènes et les valeurs locales du module d'Young sont autant d'inconnues à déterminer dans le problème posé. Utiliser des champs virtuels particuliers n'est ici pas considéré, car définir manuellement un grand nombre de CVF indépendants est une tâche très peu réalisable. Les images dont nous disposons avant d'effectuer la reconstruction des propriétés mécaniques sont des images échographiques, des images de déplacement et des image de déformation. Mais dans le cas des tissus biologiques, aucune ne permet d'effectuer avec exactitude une segmentation pour définir des régions qui seraient mécaniquement homogènes. Ainsi, l'utilisation d'approches reposant sur la connaissance préalable de régions homogènes apparaît peu pertinente dans notre contexte. Pour ces raisons, une approche utilisant des champs virtuels définis par morceaux est ici privilégiée. En effet, cette dernière a l'avantage de permettre de reconstruire des cartes des propriétés mécaniques sans connaissances sur leurs distributions spatiales.

### V.2.2 Définition des champs virtuels

Considérons donc un milieu élastique, linéaire, isotrope et de coefficient de Poisson $v = 0.49$, soumis à un chargement quasi-statique. Afin de reconstruire une carte de module d'Young, le domaine est subdivisé en de multiples sous-domaines pour lesquels une estimation du module d'Young va être effectuée. En termes pratiques, le domaine est maillé linéairement par des éléments rectangles composés de 4 nœuds. Le maillage est choisi plus large que l'échantillonnage des données, considérant que des points spatialement proches possèdent des propriétés mécaniques similaires (fig. 3.6).

<figure><img src="image_34.png" /><figcaption>FIGURE 3.6 – Illustration de la division du domaine en de multiples sous-domaines pour l'obtention d'une cartographie du module d'Young. a) numérotation des éléments du maillage en bleu correspondant à plusieurs pixels du milieu, b) numérotation des nœuds de chaque élément.</figcaption></figure>

Dans ces travaux, le choix de champs virtuels par morceaux a été fait. De manière similaire à la méthode des éléments finis, cela implique l'utilisation de fonctions, appelées fonctions de forme, permettant d'interpoler le déplacement virtuel sur l'élément. Considérons un élément quelconque *k*. Pour chaque nœud de cet élément, une numérotation locale lui est attribuée, conduisant à $n_1, n_2, n_3$ et $n_4$ (fig 3.7), et quatre fonctions de forme $N_1, N_2, N_3$ et $N_4$ sont définies telles que :
$$
\begin{align}
N_i(x_{ni}, y_{ni}) &= 1 \nonumber \\
N_i(x_{nj}, y_{nj}) &= 0 \quad \text{pour } j \neq i \tag{3.29}
\end{align}
$$

<figure><img src="image_35.png" /><figcaption>FIGURE 3.7 – Illustration d’un élément quelconque $k$ avec la numérotation locale de chacun de ces nœuds.</figcaption></figure>

Les fonctions de forme $N_1$, $N_2$, $N_3$ et $N_4$ s'écrivent :

$$
\left\{
\begin{aligned}
N_1(x, y) &= \left(1 + i - \frac{x}{\Delta_{lat}}\right) \left(1 + j - \frac{y}{\Delta_{ax}}\right) \\
N_2(x, y) &= \left(\frac{x}{\Delta_{lat}} - i\right) \left(1 + j - \frac{y}{\Delta_{ax}}\right) \\
N_3(x, y) &= \left(\frac{x}{\Delta_{lat}} - i\right) \left(\frac{y}{\Delta_{ax}} - j\right) \\
N_4(x, y) &= \left(1 + i - \frac{x}{\Delta_{lat}}\right) \left(\frac{y}{\Delta_{ax}} - j\right)
\end{aligned}
\right.
\qquad (3.30)
$$

avec $i\Delta_{lat} \le x \le (i+1)\Delta_{lat}$ et $j\Delta_{ax} \le y \le (j+1)\Delta_{ax}$.

Notons maintenant $\tilde{u}_i^*$ le déplacement virtuel au nœud $n_i$, de composantes latérale et axiale $\tilde{u}_{i1}^*$ et $\tilde{u}_{i2}^*$. Le déplacement virtuel sur l'élément, $u_e^*$, est obtenu en utilisant le déplacement virtuel de chacun de ses nœuds $\tilde{u}_i^*$ et les fonctions de forme $N_i$, tel que :

$$
\left\{
\begin{aligned}
u_{e1}^* &= \sum_{i=1}^{4} N_i \tilde{u}_{i1}^* \\
u_{e2}^* &= \sum_{i=1}^{4} N_i \tilde{u}_{i2}^*
\end{aligned}
\right.
\qquad (3.31)
$$

Autrement dit :

$$
u_e^* = N \tilde{u}_e^* \tag{3.32}
$$

$$
\text{avec } \mathbf{N} = \begin{bmatrix} N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 & 0 \\ 0 & N_1 & 0 & N_2 & 0 & N_3 & 0 & N_4 \end{bmatrix} \quad (3.33)
$$

La déformation virtuelle $\boldsymbol{\epsilon}_e^*$ est ensuite définie comme la dérivée du déplacement virtuel $\boldsymbol{u}_e^*$,
soit :

$$
\epsilon_e^* = L u_e^* \tag{3.34}
$$
avec $L$ contenant les termes de dérivés (3.35).

$$L = \begin{bmatrix} \frac{\partial N_1}{\partial x_1} & 0 & \frac{\partial N_2}{\partial x_1} & 0 & \frac{\partial N_3}{\partial x_1} & 0 & \frac{\partial N_4}{\partial x_1} & 0 \\ 0 & \frac{\partial N_1}{\partial x_2} & 0 & \frac{\partial N_2}{\partial x_2} & 0 & \frac{\partial N_3}{\partial x_2} & 0 & \frac{\partial N_4}{\partial x_2} \\ \frac{\partial N_1}{2\partial x_2} & \frac{\partial N_1}{2\partial x_1} & \frac{\partial N_2}{2\partial x_2} & \frac{\partial N_2}{2\partial x_1} & \frac{\partial N_3}{2\partial x_2} & \frac{\partial N_3}{2\partial x_1} & \frac{\partial N_4}{2\partial x_2} & \frac{\partial N_4}{2\partial x_1} \end{bmatrix} \quad (3.35)$$

### V.2.3 Application du principe des travaux virtuels

Reprenons le principe des travaux virtuels écrit pour le problème quasi-statique, dans le cas des contraintes ou déformations planes (cf. section III.2). Pour rappel :

$$\int_S (\sigma_1 \varepsilon_1^* + \sigma_2 \varepsilon_2^* + 2\sigma_6 \varepsilon_6^*) ds = \int_L T_i \cdot u_i^* dl \quad (3.36)$$

Après injection de la loi de comportement (3.28), l'équation (3.36) devient :

$$\int_S E \left( (C_1\varepsilon_1 + C_2\varepsilon_2)\varepsilon_1^* + (C_1\varepsilon_2 + C_2\varepsilon_1)\varepsilon_2^* + 2C_6\varepsilon_6\varepsilon_6^* \right) ds = \int_L T_i \cdot u_i^* dl \quad (3.37)$$

Dans notre cas particulier, les valeurs de module d'Young sont à déterminer pour chaque élément du maillage. L'équation (3.37) est donc réécrite telle que :

$$\sum_j E_j \int_{S_j} \left( (C_1\varepsilon_1 + C_2\varepsilon_2)\varepsilon_1^* + (C_1\varepsilon_2 + C_2\varepsilon_1)\varepsilon_2^* + 2C_6\varepsilon_6\varepsilon_6^* \right) ds = \int_L T_i \cdot u_i^* dl \quad (3.38)$$

Considérons un élément quelconque $k$ pour lequel la valeur du module est notée $E_k$. Comme nous avons pu le voir précédemment (cf. section IV.2), l'idée de base des champs virtuels spécifiques est de déterminer le déplacement virtuel tel que le membre de gauche de l'équation (3.38) soit égal au seul module $E_k$ à déterminer. Et la valeur de $E_k$ est alors directement égale au travail des forces extérieures (membre de droite de (3.38)), calculé avec ce déplacement virtuel spécifique. En d'autres termes, avec cette approche, estimer $E_k$ revient à trouver un déplacement virtuel, qui sélectionne la contribution de l'élément $k$ et annule les contributions de tous les autres éléments. Nous verrons par la suite qu'une certaine flexibilité peut être permise conduisant à une estimation moyenne du module d'Young sur un certain nombre d'éléments.

Pour obtenir des champs virtuels satisfaisant cette condition, les intégrales du membre de gauche de l'équation (3.38) sont réécrites sous la forme d'un système linéaire. La procédure utilisée ici est similaire à celle décrite en section IV.2. Les déplacements nodaux virtuels deviennent les inconnues à déterminer en utilisant comme condition de sélectionner uniquement la contribution de l'élément $k$. Soit :

$$DY = A \quad (3.39)$$

où $Y$ est un vecteur contenant le déplacement virtuel axial et latéral en chaque nœud du maillage (3.40). $D$ est une matrice rectangle contenant tous les termes connus, c.-à-d. autres que les déplacements virtuels nodaux (qui sont contenus dans $Y$). La matrice $D$ est donc de dimension $[n_{element} \times 2n_{noeud}]$, avec $n_{element}$ le nombre d'éléments du maillage et $n_{noeud}$ le nombre de nœuds total. $n_{noeud}$ est multiplié par 2 du fait du déplacement virtuel axial et latéral en chaque nœud. De même, $Y$ est un vecteur de dimension $[2n_{noeud} \times 1]$. Enfin, $A$ est de dimension $[n_{element} \times 1]$ et correspond au résultat souhaité (3.42) (condition imposée à ces champs).
$$
Y = \begin{pmatrix}
\tilde{u}_{11}^* \\
\tilde{u}_{12}^* \\
\tilde{u}_{21}^* \\
\tilde{u}_{22}^* \\
\vdots
\end{pmatrix}
\qquad (3.40)
$$

$$
D = \begin{bmatrix}
d_{1,1} & d_{1,2} & 0 & 0 & \dots & d_{1,p+2} & d_{1,p+3} & 0 & 0 & \dots \\
0 & d_{2,2} & d_{2,3} & 0 & \dots & 0 & d_{2,p+3} & d_{2,p+4} & 0 & \dots \\
\vdots & \ddots & & & & \vdots & & & & \\
\vdots & & & & & \vdots & & & &
\end{bmatrix}
\quad (3.41)
$$

$$
\mathbf{A} = \begin{pmatrix} 0 \\ \vdots \\ 1 \\ \vdots \\ 0 \end{pmatrix} \tag{3.42}
$$

Il est important de remarquer que l'ensemble des déplacements virtuels nodaux du milieu interviennent pour l'estimation de $E_k$.

Dans ces travaux, l'équation (3.39) est reformulée comme un problème d'optimisation avec contraintes. Plus précisément, de manière à s'affranchir des conditions aux limites non connues et à satisfaire la contrainte de champs cinématiquement admissibles, le déplacement virtuel des nœuds sur la frontière est imposé comme étant nul, à l'exception de la région où la force est appliquée. Pour rappel, dans le problème présent la force est supposée connue. De plus, la force en élastographie étant en général appliquée verticalement, le déplacement virtuel latéral est également imposé comme étant nul. Ainsi le problème de minimisation à résoudre est le suivant :

$$
\text{Minimiser} \quad \| D\mathbf{Y} - \mathbf{A} \| \tag{3.43}
$$

$$
\text{avec } \tilde{u}_{i1}^* = 0 \quad \forall i
$$

$$
\tilde{u}_{j2}^* = 0 \quad \forall j \in J
$$

avec $\| \cdot \|$ la norme euclidienne et $J$ l'ensemble des indices des nœuds appartenant à la frontière à l'exception de la région où est appliquée la force.

Cette approche a l'avantage de permettre une sélection directe des champs virtuels, sans nécessiter aucune étape supplémentaire pour optimiser ces champs. La procédure de minimisation (3.43) est réalisée grâce à CVX. Il s'agit d'un outil MATLAB conçu pour résoudre des problèmes d'optimisation convexes [32, 33]. Le déplacement virtuel ainsi déterminé est ensuite utilisé pour calculer le membre de droite de l'équation (3.37) permettant d'estimer la valeur du module d'Young $E_k$, sur l'élément $k$. Finalement, cette procédure est répétée pour chaque élément du milieu afin d'obtenir une carte du module d'Young.
V.3 Premiers résultats et modifications apportées

Lors de nos premières utilisations de la méthode, nous avons pu constater des difficultés à reconstruire correctement des cartes de module d'Young, même lorsque toutes les hypothèses étaient satisfaites, se traduisant par des variations locales du module aussi bien dans l'inclusion que dans le milieu environnant. Une illustration en est donnée en figure 3.8-a pour une simulation en contraintes planes. Ici, toutes les composantes du vecteur **A** sont imposées égales à zéro, sauf la $k^{eme}$ composante correspondant à la valeur du module d'Young $E_k$ recherchée, soit $a_k = 1$ (fig. 3.9-a). Autrement dit, le déplacement virtuel est cherché tel qu'il permette de sélectionner la contribution de l'élément **k** et d'annuler les contributions de tous les autres éléments. Afin d'améliorer l'estimation, et notamment de lisser les variations observées, une condition moins stricte sur **A** a été utilisée, c.-à-d. que l'estimation du module $E_k$ prend en compte non seulement l'élément **k** mais aussi les éléments adjacents. Dans ce cas, le module d'Young estimé correspond à une valeur moyenne sur l'ensemble des éléments utilisés. En pratique, une contribution pondérée des éléments adjacents est prise en compte grâce à l'attribution de poids aux éléments concernés, ces poids correspondant à un filtre gaussien (fig. 3.9-b). En effet, le filtre gaussien permet de privilégier la valeur de l'élément recherché, qui aura alors un poids plus important dans ce moyennage. De plus, la taille du filtre est gardée relativement petite ([:3 × 3]) afin de conserver le caractère local de l'estimation. La figure 3.8-b montre le résultat obtenu avec cette méthode sur le même cas que celui présenté en figure 3.8-a. Les variations locales ont disparu et le milieu environnant et l'inclusion apparaissent maintenant homogènes, en parfait accord avec la simulation réalisée.

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

FIGURE 3.8 – Cartes de module de Young reconstruites pour un milieu hétérogène composé de deux régions homogènes de 30 kPa (inclusion) et 10 kPa (milieu environnant), simulé en contraintes planes avec Comsol Multiphysics® [34]. La minimisation est réalisée en prenant en compte : a) uniquement l'élément où une valeur du module de Young est recherchée, b) une contribution pondérée des éléments adjacents grâce à l'utilisation de poids correspondant à ceux d'un filtre gaussien [3 × 3].
<figure><img src="image_37.png" /><figcaption>FIGURE 3.9 – Illustration des valeurs introduites dans **A** afin d'estimer des cartes de modules d'Young pour un milieu quelconque. Dans cet exemple, la valeur du module d'Young sur l'élément en gris foncé est recherchée. (a) Version initiale de la méthode : le déplacement virtuel est cherché tel qu'il permette de sélectionner la contribution d'un élément et d'annuler les contributions de tous les autres éléments, donnant ainsi une seule composante du vecteur **A** non nulle (égale à 1, en gris foncé). (b) Version de la méthode prenant en compte une contribution pondérée des éléments adjacents. La contribution des éléments est d'autant plus forte que la couleur qui lui est attribuée est foncée. Ici, les poids, dont la somme doit être égale à 1, correspondent à ceux d'un filtre gaussien. **A** comporte ainsi 9 composantes non nulles.</figcaption></figure>

## V.4 Discussion

Plusieurs précisions peuvent être apportées concernant l'implémentation de cette méthode. Comme nous l'avons vu, une première étape consiste à définir le maillage à partir duquel les calculs seront réalisés. Dans la méthode proposée, des éléments rectangles composés de 4 nœuds ont été utilisés (fig. 3.6). Il peut être noté que n'importe quel type de maillage peut être employé. Cependant, ce type de maillage comporte plusieurs avantages. En effet, il présente une grande facilité de mise en œuvre et les fonctions de forme associées sont relativement simples (3.30) et de nombre limité du fait du faible nombre de nœuds. Par conséquent, la construction de champs de déplacement virtuel ainsi définis par morceaux reste de complexité limitée.

La méthode de reconstruction proposée ici nécessite d'estimer un nombre important d'inconnues (une valeur de module d'Young par élément), ce qui suppose de déterminer autant de champs virtuels, et donc de réaliser un nombre important de calculs. Le maillage utilisé, plus large que l'échantillonnage des images de déformation (fig. 3.6-a), permet entre autre de limiter le nombre de valeurs à estimer. Cependant, un maillage trop grossier peut aussi entrainer une perte d'informations, en masquant les régions du milieu de petite taille. Ici, le maillage a été choisi de telle sorte que la taille des éléments correspondent à 4 × 4 échantillons des images de déformation.

Pour ces premiers travaux, la méthode de reconstruction a été programmée en MATLAB, et les aspects temps de calcul n'ont pas été véritablement pris en considération, l'objectif premier étant d'aboutir à une méthode capable de fournir des images où les régions avec des propriétés mécaniques différentes sont clairement identifiées. Malgré tout, certains choix ont été faits afin de limiter ce temps de calcul. Premièrement, la méthode proposée a été programmée de manière à ce que l'estimation des différentes valeurs de $E_k$ soit exécutée en parallèle. Il doit être noté que cette méthode se prête parfaitement à ce type d'approche. Deuxièmement, il peut aussi être noté que la matrice **D** et le vecteur **A** sont creux (la plupart de leurs composantes sont nulles). Cette propriété a été exploitée afin de limiter l'utilisation de la mémoire et donc le temps de calcul. Enfin, les calculs ont été réalisés sur une station de travail, comportant 24 cœurs de processeurs (Intel(R) Xeon(R) CPU E5-2697 v2 2.70GHz). A titre d'exemple, les cartes de module d'Young
reconstruites et présentées en figure 3.8 ont été obtenues en approximativement 4.3 minutes.

En reposant sur le principe des travaux virtuels, la méthode développée présente certains avantages, comparée à d'autres méthodes de la littérature. Par exemple, par rapport à Nitta and Shiina [4], il n'est pas nécessaire de dériver les termes de contraintes (liés à la déformation, au coefficient de Poisson et au module d'Young). Ainsi, aucune dérivée de la déformation et donc aucune dérivée seconde du déplacement n'est à calculer, ce qui présente un intérêt puisque l'opération de dérivation est connue pour rehausser le niveau de bruit. De plus, aucune hypothèse d'homogénéité locale du milieu, hypothèse utilisée dans un certain nombre d'approches directes pour simplifier le problème, n'est à considérer. En effet, l'hypothèse d'homogénéité locale n'est pas adaptée aux milieux biologiques, hétérogènes, d'autant plus lorsque des lésions se développent. Cette approche pourrait donc aussi bénéficier à d'autres approches d'élastographie. Pour finir, il peut être noté ici que les développements présentés dans ces travaux pourront par la suite être relativement simplement modifiés afin d'enrichir le modèle décrivant le milieu (loi de comportement) et se rapprocher ainsi des propriétés des tissus.

Les résultats concernant la méthode que nous avons développée dans ce travail de thèse seront présentés dans le chapitre suivant.
Bibliographie

[1] M. M. Doyley, P. M. Meaney, and J. C. Bamber. Evaluation of an iterative reconstruction method for quantitative elastography. Physics in Medicine and Biology, 45(6):1521–1540, 2000. I

[2] S. Le Floc'h, J. Ohayon, P. Tracqui, G. Finet, A. M. Gharib, R. L. Maurice, G. Cloutier, and R. I. Pettigrew. Vulnerable atherosclerotic plaque elasticity reconstruction based on a segmentation-driven optimization procedure using strain measurements : theoretical framework. IEEE transactions on medical imaging, 28(7):1126–1137, 2009.

[3] S. R. Mousavi, H. Rivaz, A. Sadeghi-Naini, G. J. Czarnota, and A. Samani. Breast Ultrasound Elastography Using Full Inversion-Based Elastic Modulus Reconstruction. IEEE Transactions on Computational Imaging, 3(4):774–782, 2017. I

[4] N. Nitta and T. Shiina. A Method of Tissue Elasticity Estimation Based on Three-Dimensional Displacement Vector. Japanese Journal of Applied Physics, 39:3225–3229, 2000. I, V.4

[5] C. Hoerig, J. Ghaboussi, and M. F. Insana. Data-Driven Elasticity Imaging Using Cartesian Neural Network Constitutive Models and the Autoprogressive Method. IEEE transactions on medical imaging, 38(5):1150–1160, 2019. I

[6] S. Avril, M. Grédiac, and F. Pierron. Sensitivity of the virtual fields method to noisy data. Computational Mechanics, 34(6):439–452, 2004. I, IV.2

[7] S. Avril, J. M. Huntley, F. Pierron, and D. D. Steele. 3D Heterogeneous Stiffness Reconstruction Using MRI and the Virtual Fields Method. Experimental Mechanics, 48(4):479–494, 2008. I, V

[8] Y. Mei, J. Deng, X. Guo, S. Goenezen, and S. Avril. Introducing regularization into the virtual fields method (VFM) to identify nonhomogeneous elastic property distributions. Computational Mechanics, 67(6):1581–1599, 2021. I, V

[9] M. Grédiac. Principe des travaux virtuels et identification. Comptes rendus de l'Académie des sciences, Série 2, Mécanique, physique, chimie, sciences de l'univers, sciences de la terre : 1–5, 1989. II, IV.1

[10] F. Pierron and M. Grédiac. *The Virtual Fields Method*. Springer New York, New York, NY, 2012. II, II.1, II.2, 2, IV.2, IV.3, IV.3, IV.3, V

[11] D. Notta-Cuvier. *Méthode des Champs Virtuels pour la caractérisation de comportements viscoplastiques et d'endommagement, à partir de mesures de champs mécaniques hétérogènes*. PhD thesis, Université de Valenciennes et du Hainaut-Cambrésis, 2011. II.1

[12] N. Connesson, E. H. Clayton, P. V. Bayly, and F. Pierron. Extension of the Optimised Virtual Fields Method to Estimate Viscoelastic Material Parameters from 3D Dynamic Displacement Fields. Strain, 51(2):110–134, 2015. II.2, V

[13] J. T. Katsikadelis. The BEM for Two-Dimensional Elastostatic Problems. In *The Boundary Element Method for Engineers and Scientists*, pages 209–287. Elsevier, 2016. III.1.1, III.1.2

[14] F. Pierron, S. Zhavoronok, and M. Grédiac. Identification of the through-thickness properties of thick laminated tubes using the virtual fields method. International Journal of Solids and Structures, pages 4437–4453, 2000. IV.1, V
[15] M. Grédiac, E. Toussaint, and F. Pierron. Special virtual fields for the direct determination of material parameters with the virtual fields method. 1—Principle and definition. International Journal of Solids and Structures, 39(10):2691-2705, 2002. IV.2, IV.2

[16] E. Toussaint, M. Grédiac, and F. Pierron. The virtual fields method with piecewise virtual fields. International Journal of Mechanical Sciences, 48(3):256-264, 2006. IV.2

[17] M. Grédiac, E. Toussaint, and F. Pierron. Special virtual fields for the direct determination of material parameters with the virtual fields method. 2—Application to in-plane properties. International Journal of Solids and Structures, 39(10):2707-2730, 2002. IV.2

[18] M. Grédiac and F. Pierron. Numerical issues in the virtual fields method. International Journal for Numerical Methods in Engineering, 59(10):1287-1312, 2004. IV.2

[19] X. Ma, Y. Wang, and J. Zhao. Optimized Polynomial Virtual Fields Method for Constitutive Parameters Identification of Orthotropic Bimaterials. Advances in Materials Science and Engineering, 2020:1-27, 2020. IV.2, V

[20] J.-H. Kim, F. Pierron, M. Wisnom, and K. Syed-Muhamad. Identification of the local stiffness reduction of a damaged composite plate using the virtual fields method. Composites Part A: Applied Science and Manufacturing, 38(9):2065-2075, 2007. IV.3, 3.5, V

[21] J.-H. Kim, F. Pierron, M. Wisnom, and S. Avril. Local stiffness reduction in impacted composite plates from full-field measurements. Composites Part A: Applied Science and Manufacturing, 40(12):1961-1974, 2009. IV.3

[22] T. Nguyen, J. Huntley, I. Ashcroft, P. Ruiz, and F. Pierron. A Fourier-series-based virtual fields method for the identification of 2-D stiffness distributions. International Journal for Numerical Methods in Engineering, 98(12):917-936, 2014. IV.3, V

[23] M. Grédiac and A. Vautrin. A New Method for Determination of Bending Rigidities of Thin Anisotropic Plates. Journal of Applied Mechanics, 57(4):964-968, 1990. V

[24] Y. Pannier, S. Avril, R. Rotinat, and F. Pierron. Identification of Elasto-Plastic Constitutive Parameters from Statically Undetermined Tests Using the Virtual Fields Method. Experimental Mechanics, 46(6):735-755, 2006.

[25] M. Grédiac, E. Toussaint, and F. Pierron. Special virtual fields for the direct determination of material parameters with the virtual fields method. 3. –Application to the bending rigidities of anisotropic plates. International Journal of Solids and Structures, 40(10):2401-2419, 2003.

[26] A. Giraudeau, F. Pierron, and B. Guo. An alternative to modal analysis for material stiffness and damping identification from vibrating plates. Journal of Sound and Vibration, 329(10):1653-1672, 2010.

[27] S.-h. Yoon, I. Giannakopoulos, and C. R. Siviour. Application of the Virtual Fields Method to the uniaxial behavior of rubbers at medium strain rates. International Journal of Solids and Structures, 69-70:553-568, 2015. V

[28] F. Pierron, P. V. Bayly, and D. R. Namani. Application of the Virtual Fields Method to Magnetic Resonance Elastography data. In *Society for Experimental Mechanics*, page 8, 2010. V

[29] R. Miller, A. Kolipaka, M. P. Nash, and A. A. Young. Estimation of transversely isotropic material properties from magnetic resonance elastography using the optimised virtual fields method. International Journal for Numerical Methods in Biomedical Engineering, 34(6):e2979, 2018.
[30] F. Meng, X. Zhang, J. Wang, C. Li, J. Chen, and C. Sun. 3D Strain and Elasticity Measurement of Layered Biomaterials by Optical Coherence Elastography based on Digital Volume Correlation and Virtual Fields Method. Applied Sciences, 9(7):1349, 2019. V

[31] P. E. Barbone and J. C. Bamber. Quantitative elasticity imaging: what can and cannot be inferred from strain images. Physics in Medicine and Biology, 47(12):2147-2164, 2002. V.1

[32] M. Grant and S. Boyd. CVX: Matlab software for disciplined convex programming, version 2.1. http://cvxr.com/cvx, 2014. V.2.3

[33] M. Grant and S. Boyd. Graph implementations for nonsmooth convex programs. In V.BLONDEL, S. Boyd, and H. Kimura, editors, *Recent Advances in Learning and Control*, Lecture Notes in Control and Information Sciences, pages 95–110. Springer-Verlag Limited, 2008. http://stanford.edu/~boyd/graph_dcp.html. V.2.3

[34] *Comsol Multiphysics®*. Comsol AB, Stockholm, Sweden, 2019. URL http://www.comsol.com. 3.8
<figure><img src="image_38.png" /></figure>
# Chapitre 4

## Méthode des champs virtuels : Application et résultats

### Sommaire

<table><tr><td>I</td><td>Introduction</td><td>79</td></tr><tr><td>II</td><td>Critères d'évaluation</td><td>79</td></tr><tr><td>III</td><td>Simulations numériques</td><td>80</td></tr><tr><td>III.1</td><td>Présentation des données simulées</td><td>80</td></tr><tr><td>III.2</td><td>Résultats de simulations obtenus dans les cas des contraintes et déformations planes</td><td>81</td></tr><tr><td>III.2.1</td><td>Résultats en contraintes planes</td><td>81</td></tr><tr><td>III.2.2</td><td>Résultats en déformations planes</td><td>82</td></tr><tr><td>III.3</td><td>Résultats de simulations 3D</td><td>83</td></tr><tr><td>III.3.1</td><td>Étude de l'influence des propriétés mécaniques</td><td>83</td></tr><tr><td>III.3.2</td><td>Étude de l'influence de la géométrie interne</td><td>87</td></tr><tr><td>III.3.3</td><td>Multi-hétérogénéités</td><td>89</td></tr><tr><td>III.3.4</td><td>Étude de l'influence des conditions aux frontières</td><td>90</td></tr><tr><td>IV</td><td>Données expérimentales</td><td>93</td></tr><tr><td>IV.1</td><td>Milieux utilisés</td><td>93</td></tr><tr><td>IV.2</td><td>Résultats obtenus sur fantômes</td><td>93</td></tr><tr><td>IV.3</td><td>Résultats sur tissus mammaires in vivo</td><td>97</td></tr><tr><td>V</td><td>Discussion</td><td>98</td></tr><tr><td>VI</td><td>Conclusion</td><td>101</td></tr></table>
# I Introduction

Ce dernier chapitre est consacré à la mise en œuvre de la méthode de reconstruction introduite précédemment. Cette méthode a été développée dans le but d'estimer la distribution spatiale du module d'Young au sein de tissus à partir du tenseur de déformation obtenu lors d'expériences d'élastographie ultrasonore quasi-statique. Une évaluation de la méthode sera réalisée à l'aide de données de complexité croissante, allant de simulations numériques à des images de tissus mammaires in vivo. Les simulations ont l'avantage de permettre d'étudier la précision des cartes reconstruites car la distribution spatiale réelle du module d'Young est connue. De plus, une étude de l'influence de différents paramètres, tels que les valeurs des propriétés mécaniques ou les conditions aux frontières, pourra être réalisée sur les milieux simulés. Les données expérimentales seront ensuite utilisées pour déterminer si la méthode est toujours capable de produire des images du module d'Young dans des conditions d'acquisitions réelles. Des tests sur fantômes seront d'abord présentés, constituant une étape intermédiaire pour l'évaluation de la méthode du fait de leur moins grande complexité. Pour ces milieux, des informations sur leurs propriétés mécaniques sont fournies par le fabricant, ce qui aidera à analyser les résultats. Enfin, trois reconstructions effectuées sur des lésions mammaires seront présentées.

Le chapitre est organisé de la manière suivante : les critères d'évaluation seront d'abord introduits. Les différentes simulations réalisées seront ensuite présentées ainsi que les résultats obtenus. Puis viendront les résultats issus de données expérimentales sur fantômes et tissus in vivo avec les cartographies reconstruites. Une dernière partie sera dédiée à la discussion des résultats, suivie de la conclusion.

# II Critères d'évaluation

Dans ces travaux, différents critères d'évaluation seront utilisés afin d'évaluer la méthode développée. Un premier critère est visuel et concerne la capacité de la méthode à détecter des hétérogénéités au sein d'un milieu ainsi que l'aspect de ces différentes régions. Des critères plus quantitatifs sont également introduits, dépendant du type de données analysées.

Pour les simulations numériques, la géométrie des milieux ainsi que leurs propriétés mécaniques sont choisies par l'utilisateur et donc connues avec exactitude. Les valeurs de module d'Young estimées pourront donc être directement comparées à celles introduites dans les modèles. Ainsi, sur les cartographies reconstruites, le module d'Young moyen sera calculé dans le milieu environnant ($E_b$) et dans les inclusions ($E_i$) en excluant la zone de transition. Ces valeurs seront accompagnées de leur écart type, permettant de quantifier la dispersion des valeurs du module au sein des différentes régions. Enfin, un rapport des modules élastiques moyens entre l'inclusion et le milieu environnant $R$ sera également calculé, tel que :

$$R = \frac{E_i}{E_b} \qquad (4.1)$$

Les inclusions utilisées dans nos simulations étant circulaires, une ROI de même forme sera employée pour la sélection de ces régions spécifiques. Pour le milieu environnant, deux approches seront considérées (fig. 4.1). La première consistera à utiliser une ROI circulaire, comme pour l'inclusion (configuration 1). Cette configuration est intéressante car elle correspond à la pratique actuelle des radiologues lors d'examens cliniques. Cependant, le fait de choisir manuellement une ROI a inévitablement un impact sur le résultat. C'est pourquoi, une deuxième configuration sera considérée, qui consistera à prendre en compte l'ensemble du milieu environnant (configuration 2).
<figure><img src="image_39.png" /><figcaption>FIGURE 4.1 – Illustration de la sélection de régions dans le milieu environnant (rouge) et l'inclusion (bleu) dans la configuration 1 (a) et 2 (b).</figcaption></figure>

De plus, la racine de l'erreur quadratique moyenne (RMSE, (4.2)) pourra être estimée afin d'évaluer l'erreur commise entre les valeurs de module d'Young reconstruites et les valeurs théoriques.

$$RMSE = \sqrt{\frac{1}{N} \sum_{i=1}^{N} (E_{mes,i} - E_{sim,i})^2} \quad (4.2)$$

avec $E_{mes,i}$ la ième valeur du module d'Young estimée dans la région considérée (composée de N points) et $E_{sim,i}$ le module d'Young réel.

Pour les données expérimentales, la valeur de la force appliquée sur le milieu est inconnue. Par conséquent, seules des cartes relatives pourront être reconstruites, en attribuant une valeur arbitraire à la force. Afin d'analyser les résultats obtenus, le rapport des modules élastiques entre l'inclusion et le milieu environnant R (4.1) sera estimé. Les régions d'analyse seront sélectionnées de manière similaire aux simulations (fig. 4.1). Dans le cas où l'ensemble du milieu environnant est sélectionné, les bords de l'image pourront être exclus en présence d'artéfacts.

# III Simulations numériques

Une première évaluation de la méthode a été effectuée à l'aide de simulations numériques réalisées avec Comsol Multiphysics® [1]. Ces simulations représentent différents niveaux de complexité, allant de données respectant les conditions des contraintes ou déformations planes, à des données 3D plus proches des conditions d'acquisition réelles. En effet, les simulations 3D ne satisfont pas les hypothèses de déformations et contraintes planes, et constituent ainsi, pour l'évaluation de la méthode, une étape intermédiaire vers les données expérimentales. Les simulations réalisées sont présentées ci-dessous ainsi que les résultats obtenus.

## III.1 Présentation des données simulées

L'ensemble des milieux simulés sont élastiques, linéaires, et isotropes.

Les premières simulations sont construites afin de satisfaire les conditions de contraintes planes ou de déformations planes (Chap. 3). Ainsi des milieux de dimensions $60 \times 40 \times e$ mm³ ont été simulés, où l'épaisseur e est respectivement égale à 1 mm ou $10^3$ mm. Pour chaque cas, deux situations sont étudiées : celle d'un milieu homogène possédant un module d'Young de 20 kPa, et celle d'un milieu hétérogène constitué d'un milieu environnant homogène de 10 kPa dans lequel se trouve une inclusion de 30 kPa. Dans tous les cas, le coefficient de Poisson est fixé à 0.49. La simulation de la déformation des milieux a été effectuée de la façon suivante : leur surface supérieure est soumise à une force verticale F, tandis que la surface inférieure est fixe dans la direction verticale et libre de s'étendre horizontalement. Les surfaces verticales sont libres de tout mouvement. Un récapitulatif des différentes simulations est proposé dans la table 4.1.
<table><caption>TABLE 4.1 – Récapitulatif des milieux simulés en 2D</caption><thead><tr><th rowspan="2">cas</th><th colspan="2">Hypothèses</th><th colspan="2">Milieu</th><th rowspan="2">Force [N]</th><th rowspan="2">Module d'Young [kPa]</th><th rowspan="2">Coefficient de Poisson</th></tr><tr><th>contraintes planes</th><th>déformations planes</th><th>homogène</th><th>hétérogène</th></tr></thead><tbody><tr><td>#1.1</td><td>x</td><td></td><td>x</td><td></td><td>-0.05</td><td>20</td><td>0.49</td></tr><tr><td>#1.2</td><td></td><td>x</td><td>x</td><td></td><td>-50</td><td>20</td><td>0.49</td></tr><tr><td>#1.3</td><td>x</td><td></td><td></td><td>x</td><td>-0.05</td><td>10-30</td><td>0.49</td></tr><tr><td>#1.4</td><td></td><td>x</td><td></td><td>x</td><td>-50</td><td>10-30</td><td>0.49</td></tr></tbody></table>

Dans un deuxième temps, des simulations 3D ont été réalisées. Ces simulations reproduisent la configuration d'un milieu environnant homogène de dimensions 60 × 40 × 50 mm³, dans lequel se trouve(nt) une ou plusieurs inclusions. Différentes simulations sont réalisées en faisant varier les propriétés mécaniques des différentes régions, la taille et la position des inclusions, ainsi que la force appliquée. Par souci de lisibilité, ces paramètres seront indiqués au début de chaque section concernée. La simulation de la déformation des milieux a été effectuée de la même manière que pour les simulations précédentes (fig. 4.2). Pour finir, la méthode de reconstruction développée dans ce travail étant 2D, seuls les champs de déformation du plan médian ($x_3 = L_z/2$) ont été retenus.

<figure><img src="image_40.png" /><figcaption>FIGURE 4.2 – Schéma d'un milieu simulé en 3D, soumis à un chargement vertical et uniforme appliqué en $x_2 = L_y$ (bleu) sur l'ensemble de sa surface supérieure.</figcaption></figure>

### III.2 Résultats de simulations obtenus dans les cas des contraintes et déformations planes

Les premiers tests ont été réalisés à l'aide de simulations satisfaisant les conditions de contraintes planes ou de déformations planes. Les résultats obtenus sont présentés ci-dessous.

#### III.2.1 Résultats en contraintes planes

La figure 4.3 présente les cartes de module d'Young reconstruites pour des milieux simulés en contraintes planes (cf. cas #1.1 et cas #1.3, tab. 4.1). Sur cette figure, l'image de gauche est la cartographie du module d'Young reconstruite dans le cas du milieu homogène et l'image de droite la cartographie obtenue pour le milieu hétérogène. Dans les deux cas, les cartes reconstruites montrent une distribution spatiale du module d'Young en parfait accord avec les simulations réalisées, et la valeur du module sur les différentes régions apparaît relativement constante. De plus, que le milieu environnant soit sélectionné sur une ROI ou sur son ensemble, les mêmes résultats sont obtenus. Pour le milieu homogène (fig. 4.3a), le module d'Young moyen est très précisément estimé, à 20.00 ± 0.03 kPa, pour une valeur théorique de 20 kPa. Pour le milieu hétérogène (fig. 4.3b), l'inclusion est clairement visible et apparaît plus rigide que le milieu environnant, avec un module d'Young moyen estimé à 30.11 ± 0.57 kPa vs 10.00 ± 0.10 kPa, les valeurs étant là encore extrêmement proches des valeurs introduites dans la modélisation, respectivement 30 kPa et 10 kPa. Les faibles valeurs d'écart type observées dans les deux cas
confirment le caractère homogène des différentes régions reconstruites. De plus, la précision des estimations conduit à de très faibles valeurs d'erreur (RMSE, (4.2)) : 0.03 kPa pour le milieu homogène, et pour le milieu hétérogène, 0.10 kPa pour le milieu environnant et 0.58 kPa pour l'inclusion.

<figure><img src="image_41.png" /><figcaption>FIGURE 4.3 – Cartes de module d'Young reconstruites pour des simulations en contraintes planes d'un milieu a) homogène de 20 kPa et b) hétérogène constitué d'une inclusion de 30 kPa localisée au centre d'un milieu environnant de 10 kPa.</figcaption></figure>

### III.2.2 Résultats en déformations planes

Les résultats obtenus ensuite pour les mêmes milieux, mais dans le cas des déformations planes (cf. cas #1.2 et cas #1.4, tab. 4.1) sont présentés en figure 4.4. Ici, les mêmes observations que précédemment peuvent être faites. Pour le milieu homogène (a), la carte reconstruite est uniforme avec un module d'Young estimé à 20.00±0.03 kPa (RMSE 0.03 kPa), ce qui est en parfait accord avec la simulation. Pour le milieu hétérogène (b), l'inclusion est clairement visible et plus rigide que le milieu environnant, avec un module d'Young de 30.08±1.08 kPa vs 10.00±0.07 kPa pour le milieu environnant. Ces modules correspondent, là encore, parfaitement aux valeurs introduites dans la simulation. La RMSE est de 0.07 kPa pour le milieu environnant et légèrement plus grande (1.08 kPa) pour l'inclusion. La distribution spatiale du module d'Young est donc correctement reconstruite et les différentes régions apparaissent homogènes.

<figure><img src="image_42.png" /><figcaption>FIGURE 4.4 – Cartes de module d'Young reconstruites pour des simulations en déformations planes d'un milieu a) homogène de 20 kPa et b) hétérogène constitué d'une inclusion de 30 kPa localisée au centre d'un milieu environnant de 10 kPa.</figcaption></figure>

Ces exemples montrent que la méthode proposée permet de reconstruire avec précision la distribution du module d'Young pour des milieux satisfaisant les hypothèses de contraintes et de déformations planes. Des tests supplémentaires ont été réalisés avec d'autres types de milieux comportant une ou plusieurs inclusions avec des tailles, positions mais aussi des modules d'Young différents. Les résultats, non présentés dans ce manuscrit de thèse, ont là encore montré une estimation précise du module d'Young [2, 3]. Il apparaît donc clairement que lorsque les conditions introduites dans la méthode sont remplies, la distribution spatiale du module d'Young est parfaitement reconstruite. Cependant, les hypothèses de contraintes planes ou de déformations planes ne peuvent être satisfaites expérimentalement, et c'est pourquoi, dans la suite de ce chapitre,
l'analyse des résultats portera spécifiquement sur les simulations 3D, plus proches des données réelles.

### III.3 Résultats de simulations 3D

Des tests ont ensuite été réalisés pour des simulations 3D reproduisant des milieux dans différentes configurations. Les résultats sont présentés ci-dessous.

#### III.3.1 Étude de l'influence des propriétés mécaniques

Dans un premier temps, nous nous intéressons à l'influence des propriétés mécaniques du milieu étudié sur les cartes reconstruites. Les valeurs utilisées sont récapitulées dans la table suivante :

<table><caption>TABLE 4.2 – Valeurs des propriétés mécaniques des milieux étudiés. Les deux lignes sont indépendantes.</caption><tr><td>Module d'Young [kPa] (E<sub>b</sub> – E<sub>i</sub>)</td><td>25-10</td><td>10-5</td><td>10-20</td><td>10-30</td><td>10-50</td><td></td><td></td></tr><tr><td>Coefficient de Poisson</td><td>0.499</td><td>0.495</td><td>0.49</td><td>0.47</td><td>0.45</td><td>0.4</td><td>0.1</td></tr></table>

La figure 4.5 présente les résultats obtenus avec la méthode utilisant les hypothèses de contraintes planes (première colonne) ou de déformations planes (deuxième colonne), pour des milieux 3D possédant différents modules d'Young et un coefficient de Poisson de 0.49. Dans l'ensemble des cas, l'inclusion est clairement visible sur les cartes reconstruites, et les deux régions apparaissent relativement homogènes. En contraintes planes, des variations légèrement plus importantes du module d'Young peuvent être observées sous les inclusions, restant malgré tout cohérent avec les valeurs théoriques. En revanche, en déformations planes, les cartes reconstruites présentent des valeurs beaucoup plus faibles, environ dix fois inférieures aux valeurs théoriques. Enfin, dans les deux cas, une diminution du contraste élastique peut être observée par rapport aux valeurs réelles. Il peut également être noté que la diminution du contraste élastique est plus importante pour les résultats obtenus à partir des hypothèses de déformations planes (fig. 4.5-f). Néanmoins, il reste possible de distinguer deux inclusions aux modules d'Young différents dans les deux cas.
<figure><img src="image_43.png" /><figcaption>FIGURE 4.5 – Résultats pour des milieux de différents modules d'Young. Les modules du milieu environnant et de l'inclusion sont respectivement ($E_b - E_i$) : 25 – 10 kPa (a1,a2), 10 – 5 kPa (b1,b2), 10 – 20 kPa (c1,c2), 10 – 30 kPa (d1,d2) et 10 – 50 kPa (e1,e2). Le coefficient de Poisson est égal à 0.49. A gauche : cartes de module d'Young reconstruites en utilisant les hypothèses de contraintes planes (a1-e1) et de déformations planes (a2-e2). A droite : courbe présentant le rapport entre le module d'Young moyen de l'inclusion et celui du milieu environnant, obtenu à partir des images a1-e1 (bleu), a2-e2 (rouge), et des valeurs théoriques (noir).</figcaption></figure>

<table><caption>TABLE 4.3 – Modules d'Young moyens [kPa] et écarts types estimés pour des régions sélectionnées dans l'inclusion et le milieu environnant sur les cartographies présentées figure 4.5. Pour le milieu environnant, la ligne grise correspond à une estimation du module d'Young moyen en tenant compte de l'ensemble du milieu (configuration 2) et l'autre en sélectionnant une ROI circulaire (configuration 1). Ces résultats ont été obtenus avec la méthode utilisant l'hypothèse des contraintes planes.</caption><thead><tr><th></th><th>a1</th><th>b1</th><th>c1</th><th>d1</th><th>e1</th></tr></thead><tbody><tr><td>E<sub>b</sub> - E<sub>i</sub> simulé</td><td>25-10</td><td>10-5</td><td>10-20</td><td>10-30</td><td>10-50</td></tr><tr><td rowspan="2">E<sub>b</sub> estimé</td><td>24.96 ± 0.10</td><td>9.99 ± 0.03</td><td>10.02 ± 0.04</td><td>10.03 ± 0.07</td><td>10.05 ± 0.10</td></tr><tr><td>25.35 ± 0.56</td><td>10.11 ± 0.17</td><td>9.90 ± 0.17</td><td>9.84 ± 0.26</td><td>9.79 ± 0.34</td></tr><tr><td>E<sub>i</sub> estimé</td><td>15.45 ± 0.34</td><td>6.84 ± 0.11</td><td>16.09 ± 0.23</td><td>22.05 ± 0.59</td><td>33.87 ± 1.70</td></tr></tbody></table>
<table><thead><tr><th></th><th>a2</th><th>b2</th><th>c2</th><th>d2</th><th>e2</th></tr></thead><tbody><tr><td>E<sub>b</sub> - E<sub>i</sub> simulé</td><td>25-10</td><td>10-5</td><td>10-20</td><td>10-30</td><td>10-50</td></tr><tr><td>E<sub>b</sub> estimé</td><td>2.65 ± 0.04</td><td>1.07 ± 0.01</td><td>1.15 ± 0.02</td><td>1.17 ± 0.03</td><td>1.20 ± 0.05</td></tr><tr><td></td><td>2.78 ± 0.15</td><td>1.11 ± 0.05</td><td>1.10 ± 0.05</td><td>1.10 ± 0.08</td><td>1.11 ± 0.11</td></tr><tr><td>E<sub>i</sub> estimé</td><td>1.99 ± 0.02</td><td>0.85 ± 0.01</td><td>1.60 ± 0.01</td><td>2.10 ± 0.02</td><td>3.08 ± 0.06</td></tr></tbody></table>

TABLE 4.4 – Modules d'Young moyens [kPa] et écarts types estimés pour des régions sélectionnées dans l'inclusion et le milieu environnant sur les cartographies présentées figure 4.5. Pour le milieu environnant, la ligne grise correspond à une estimation du module d'Young moyen en tenant compte de l'ensemble du milieu (configuration 2) et l'autre en sélectionnant une ROI circulaire (configuration 1). Ces résultats ont été obtenus avec la méthode utilisant l'hypothèse des déformations planes.

Les modules d'Young moyens estimés dans l'inclusion et le milieu environnant sont présentés en tables 4.3 et 4.4. Concernant les résultats en contraintes planes, si dans le milieu environnant le module d'Young est correctement estimé, avec des valeurs moyennes comprises dans la gamme [9.79-10.11] vs 10 kPa pour la valeur théorique, et pour le cas du milieu environnant de plus forte valeur de module 24.96 ou 25.35 kPa vs 25 kPa, il est sur- ou sous-estimé (en fonction des cas) dans l'inclusion. Par exemple, dans le cas a1, où le module d'Young de l'inclusion dans le modèle est de 10 kPa, sa valeur reconstruite a une valeur moyenne de 15.45 kPa. Pour le cas c1, où le module d'Young théorique est 20 kPa, sa valeur reconstruite moyenne est de 16.09 kPa. Une estimation moins précise du module d'Young dans l'inclusion est donc à l'origine de la diminution du contraste élastique. Cette observation est confirmée par des valeurs plus élevées de la RMSE concernant les inclusions, avec erreurs comprises entre [1.84-16.22] kPa vs [0.04-0.66] kPa pour le milieu environnant. Enfin, il peut être noté que les écarts types estimés sont faibles, avec des valeurs dans la gamme [0.03-0.56] kPa pour le milieu environnant et [0.11-1.70] kPa pour l'inclusion. Bien que ces valeurs soient légèrement plus élevées pour l'inclusion, elles restent suffisamment faibles pour que les régions observables sur les cartes reconstruites apparaissent homogènes.

Concernant les résultats en déformations planes, les modules d'Young estimés sont environ dix fois plus faibles. Cette observation est confirmée par une plus grande RMSE dans les deux régions, avec des valeurs comprises dans la gamme [4.15-46.92] kPa. Cependant, les valeurs d'écarts types estimées restent faibles à la fois dans le milieu environnant et dans l'inclusion. Celle-ci apparaît plus ou moins rigide par rapport au milieu environnant en accord avec les simulations.

Il peut être noté que les deux configurations utilisées pour sélectionner les différentes régions produisent des résultats similaires. Cela s'explique par le fait que le milieu environnant apparaît très homogène sur les cartes reconstruites. Par la suite, dans de tels cas, seule la première configuration sera utilisée. Enfin, que ce soit dans le cas où la méthode utilise les hypothèses de contraintes planes ou de déformations planes, la diminution du contraste élastique pouvait être attendue, car les simulations 3D ne satisfont aucune de ces deux hypothèses. Cependant, les résultats obtenus sous l'hypothèse de contraintes planes sont plus proches des valeurs réelles, à la fois en termes de module d'Young moyen et de rapport élastique. Pour ces raisons, cette hypothèse sera privilégiée dans la suite des tests menés.

Dans la méthode développée, le coefficient de Poisson est supposé connu et a été fixé à 0.49, considérant que les tissus mous sont quasi-incompressibles. Il ne reste alors que la distribution du module d'Young à déterminer. Cependant, on peut s'interroger sur l'impact sur le résultat d'un coefficient de Poisson qui serait différent de la valeur introduite dans l'algorithme. Afin d'étudier l'influence de ce coefficient, des simulations 3D ont été réalisées en faisant varier sa valeur (tab.
4.2) alors que les modules d'Young sont restés constants, respectivement égaux à 10 kPa et 30 kPa pour le milieu environnant et l'inclusion. Bien entendu, dans la méthode, le coefficient de Poisson est resté inchangé (fixé à 0.49).

<figure><img src="image_44.png" /><figcaption>FIGURE 4.6 – Influence du coefficient de Poisson sur les résultats. a) exemple de carte de module d'Young reconstruite pour un milieu ayant un coefficient de Poisson de 0.49, b) profils horizontaux et verticaux traversant l'inclusion en son centre, pour l'ensemble des milieux étudiés.</figcaption></figure>

Dans tous les cas examinés, les cartes de module d'Young reconstruites sont similaires et seul un exemple est proposé en figure 4.6-a. Les profils horizontaux et verticaux (fig. 4.6-b) obtenus pour les différents milieux sont tous quasi-identiques, sauf pour la valeur de $\nu = 0.1$ pour laquelle une faible diminution du module d'Young par rapport aux autres cas peut être remarquée au niveau de l'inclusion. Les valeurs moyennes du module Young et les écarts types estimés sont présentés en figure 4.7-a. Il peut être observé que le coefficient de Poisson du milieu étudié n'a que peu d'impact sur les résultats. De même, l'erreur quadratique moyenne reste relativement stable dans tous les cas observés (fig. 4.7-b). Comme discuté précédemment, une valeur d'erreur plus importante peut être observée dans l'inclusion du fait de l'utilisation de l'hypothèse des contraintes planes. Un module d'Young légèrement plus faible peut malgré tout être observé pour $\nu = 0.1$, avec une augmentation de l'erreur quadratique moyenne. Il est à noter qu'une telle valeur de coefficient de Poisson n'est pas attendue pour les tissus étudiés en élastographie, mais elle a été utilisée ici afin d'évaluer la méthode développée.
<figure><img src="image_45.png" /><figcaption>FIGURE 4.7 – Module d'Young moyen et écart type (a), et erreur RMSE (b) estimés dans l'inclusion (rouge) et le milieu environnant (bleu) en fonction du coefficient de Poisson des milieux étudiés.</figcaption></figure>

### III.3.2 Étude de l'influence de la géométrie interne

Dans le cadre de l'évaluation de la méthode, nous avons souhaité étudier si la taille ou la position des inclusions au sein des milieux pouvait avoir un impact sur les cartes reconstruites. Afin d'étudier ces paramètres, deux configurations ont été simulées reproduisant un milieu hétérogène composé d'une inclusion de module d'Young 30 kPa dans un milieu environnant de 10 kPa, et de coefficient de Poisson 0.49. Pour la première configuration, l'inclusion possède un diamètre variable de 12, 16 ou 20 mm. Pour la deuxième configuration, la position de l'inclusion au sein du milieu varie.

La figure 4.8 présente les résultats pour des milieux possédant une inclusion dont le diamètre varie. La présence de l'inclusion est clairement visible sur les cartes de module d'Young reconstruites et une augmentation de sa taille peut être observée (fig. 4.8-a). De plus, les valeurs moyennes de module d'Young sont similaires dans chacun des cas, avec des valeurs comprises entre [22.05 – 22.77] kPa pour l'inclusion et [10.02 – 10.03] kPa pour le milieu environnant (tab. 4.5). De même, les valeurs d'écart type restent faibles, traduisant l'homogénéité des régions dans l'ensemble des cas. Concernant la taille de ces inclusions, les profils tracés en figure 4.8-b montrent l'augmentation du diamètre aussi bien horizontalement que verticalement. Estimer le diamètre des inclusions à partir de ces cartes est une tâche peu précise. En effet, comme cela a été indiqué pour la sélection de ROI, déterminer la frontière entre l'inclusion et le milieu environnant est complexe et source d'erreurs. Cependant, afin d'obtenir une valeur indicative, le parti pris ici a été de segmenter les deux régions en fonction de la valeur moyenne entre le module d'Young estimé dans le milieu environnant et celui estimé dans l'inclusion. Les diamètres ainsi estimés sont proches des valeurs réelles aussi bien horizontalement que verticalement (fig. 4.8-c).
<figure><img src="image_46.png" /><figcaption>FIGURE 4.8 – Résultats obtenus à partir de milieux simulés contenant une inclusion de 12, 16 ou 20 mm de diamètre : a) cartes de module d'Young reconstruites, b) profils horizontaux (bleu) et verticaux (rouge), c) diamètre de l'inclusion estimé à partir des profils horizontaux (bleu) et verticaux (rouge), et valeurs théoriques (noir).</figcaption></figure>

<table><caption>TABLE 4.5 – Module d'Young [kPa] estimé dans l'inclusion et le milieu environnant en fonction du diamètre de l'inclusion (fig. 4.8).</caption><thead><tr><th>diamètre [mm]</th><th>12</th><th>16</th><th>20</th></tr></thead><tbody><tr><th><i>E</i><sub>b</sub> estimé</th><td>10.03 ± 0.07</td><td>10.03 ± 0.14</td><td>10.02 ± 0.37</td></tr><tr><th><i>E</i><sub>i</sub> estimé</th><td>22.05 ± 0.59</td><td>22.35 ± 0.98</td><td>22.77 ± 0.86</td></tr></tbody></table>

Concernant la position des inclusions, les cartes reconstruites en figure 4.9 révèlent que ces dernières sont clairement visibles et homogènes quelle que soit leur position dans le milieu. Encore une fois, les valeurs d'écart type sont relativement faibles, bien que légèrement plus élevées dans l'inclusion (tab. 4.6). Enfin, que ce soit pour l'inclusion ou le milieu environnant, les modules d'Young estimés sont similaires dans toutes les configurations, même si une légère diminution de la valeur du module d'Young estimé peut être observée avec la profondeur.

<figure><img src="image_47.png" /><figcaption>FIGURE 4.9 – Cartes de module d'Young reconstruites pour des milieux contenant une inclusion à différentes positions.</figcaption></figure>

<table><caption>TABLE 4.6 – Module d'Young estimé [kPa] en fonction de la position de l'inclusion (fig. 4.9).</caption><thead><tr><th></th><th>a</th><th>b</th><th>c</th></tr></thead><tbody><tr><th><i>E</i><sub>b</sub> estimé</th><td>10.03 ± 0.07</td><td>10.11 ± 0.06</td><td>9.96 ± 0.07</td></tr><tr><th><i>E</i><sub>i</sub> estimé</th><td>22.05 ± 0.59</td><td>22.69 ± 0.96</td><td>21.45 ± 0.66</td></tr></tbody></table>

Pour conclure, la reconstruction des cartes de module d'Young n'est pas affectée par des variations de tailles ou de positions des inclusions.
### III.3.3 Multi-hétérogénéités

Enfin, la méthode développée est évaluée grâce à des simulations contenant plusieurs hé- térogénéités/inclusions différentes. Ce type de simulations a l'avantage d'être plus complexe et de combiner les études précédentes. Le milieu environnant est homogène et possède un module d'Young de 10 kPa, ainsi qu'un coefficient de Poisson de 0.49. Une première inclusion de 12 mm, située en haut à gauche du milieu, possède un module d'Young de 30 kPa et un coefficient de Poisson de 0.49. Une deuxième inclusion, dont les paramètres varient et sont récapitulés en table 4.7, est située en bas à droite du milieu.

<table><caption>TABLE 4.7 – Récapitulatif des propriétés de l'inclusion située en bas à droite du fantôme.</caption><thead><tr><td></td><td>Module d'Young [kPa]</td><td>Coefficient de Poisson</td><td>Diamètre [mm]</td></tr></thead><tbody><tr><td>a</td><td>30</td><td>0.49</td><td>12</td></tr><tr><td>b</td><td>30</td><td>0.45</td><td>12</td></tr><tr><td>c</td><td>20</td><td>0.49</td><td>12</td></tr><tr><td>d</td><td>30</td><td>0.49</td><td>8</td></tr></tbody></table>

La figure 4.10 présente les cartes reconstruites pour les différents cas indiqués dans la table 4.7. Dans tous les cas, les inclusions sont clairement visibles et les différentes régions apparaissent homogènes. Pour les cas (c) et (d), les inclusions du bas sont visuellement différentes de celles du haut. Elles sont respectivement moins rigide (c) et plus petite (d), ce qui est cohérent avec les simulations réalisées. En revanche, pour les deux premiers cas, aucune différence n'est observable. Le module d'Young moyen estimé dans les inclusions et le milieu environnant, ainsi que les écarts types, sont présentés en figure 4.11. Dans tous les cas, les valeurs estimées dans le milieu environnant sont similaires et environ égales à 10 kPa. Cependant, des petites différences peuvent être notées pour les inclusions. En effet, pour les cas (a), (b), et (d) la valeur moyenne du module d'Young de l'inclusion basse est légèrement plus faible que celle du haut, ce qui est en accord avec les observations précédentes. Cette différence est légèrement plus marquée pour le cas (b), ce qui peut être lié au coefficient de Poisson plus faible dans cette inclusion (0.45) que dans le reste du milieu (0.49). Pour finir, concernant le cas (c), les deux inclusions de module d'Young différents sont clairement distinguables. De plus, leurs valeurs moyennes estimées (environ 22 et 16 kPa) sont similaires aux valeurs observées précédemment.

<figure><img src="image_48.png" /><figcaption>FIGURE 4.10 – Cartes de module d'Young reconstruites pour des milieux contenant plusieurs inclusions aux propriétés mécaniques différentes (tab. 4.7).</figcaption></figure>
<figure><img src="image_49.png" /></figure>

FIGURE 4.11 – Modules d'Young moyens et écarts types estimés dans les inclusions (rouge, l'inclusion supérieure est représentée par un cercle et celle inférieure par un diamant) et dans le milieu environnant (bleu).

### III.3.4 Étude de l'influence des conditions aux frontières

Pour finir, l'influence des conditions aux frontières des milieux sur les cartes reconstruites a été étudiée. Pour cela, un premier milieu 3D soumis à différentes valeurs de force (5, 2, 0.5, -0.5, -2, -5 N) a été simulé. Les valeurs positives correspondent à une décompression et celles négatives à une compression. La force est appliquée verticalement sur la surface supérieure et comme nous l'avions décrit lors de la présentation des simulations 3D (fig. 4.13-a), la surface inférieure est fixe dans la direction verticale et libre de s'étendre horizontalement. Quant aux surfaces verticales, elles sont libres de tout mouvement. Le milieu étudié est toujours composé d'une inclusion sphérique de 30 kPa et d'un milieu environnant de 10 kPa, les deux régions ayant un coefficient de Poisson de 0.49.

La figure 4.12 présente les valeurs de module d'Young moyen et d'écart type reconstruites dans l'inclusion et le milieu environnant pour les différentes simulations. Il est possible d'observer que ces valeurs sont stables/similaires pour tous les cas étudiés. Ni la valeur de la force ni son sens (compression ou décompression) ne semble donc avoir d'impact sur les cartes de module d'Young reconstruites. Il doit être noté ici que ces valeurs de force mènent à des champs de déformations dans des ordres de grandeur observables en élastographie.

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

FIGURE 4.12 – Modules d'Young moyens et écarts types estimés dans l'inclusion (rouge) et dans le milieu environnant (bleu) en fonction de la force appliquée sur le milieu.

Le même milieu (mêmes propriétés mécaniques) est ensuite repris avec des conditions aux limites différentes. Ainsi, deux configurations additionnelles ont été simulées. Pour l'une d'elle, la surface inférieure qui était initialement uniquement contrainte verticalement, est maintenant
fixe dans les 3 directions (fig. 4.13-a). Enfin, dans la dernière configuration considérée, la force est appliquée sur une zone restreinte afin d'imiter une compression réalisée à l'aide de la sonde échographique (fig. 4.13-b). La surface inférieure est uniquement contrainte verticalement.

<figure><img src="image_51.png" /><figcaption>FIGURE 4.13 – Schématisation des milieux simulés en 3D, soumis à un chargement vertical et uniforme appliqué en $x_2 = L_y$ (bleu) sur a) l'ensemble de la surface, b) un rectangle positionné au centre de la surface afin de simuler une compression/décompression réalisée à l'aide d'une sonde échographique.</figcaption></figure>

La méthode développée permet d'ignorer les conditions aux limites des milieux étudiés grâce à l'utilisation de champs virtuels imposés comme nuls aux frontières. Cependant, il est possible de voir, pour les résultats présentés en figure 4.14, que ces conditions peuvent malgré tout avoir une influence sur les cartes de module d'Young reconstruites. En effet, pour la deuxième configuration (a2), une "zone rouge" est visible en bas de l'image ce qui n'était pas le cas dans la première configuration (a1). Néanmoins, l'inclusion est clairement visible et les valeurs estimées (module d'Young, écart type et RMSE) sont similaires dans les deux configurations (tab. 4.8). En revanche, pour la dernière configuration (a3), si l'inclusion est aussi visible sur la carte reconstruite, des variations du module d'Young sont observables dans le milieu environnant. Ces variations se traduisent par une augmentation de la valeur d'écart type estimée (tab. 4.8). De plus, les valeurs de module estimées sont beaucoup plus fortes que celles attendues, entrainant une augmentation de la RMSE. Néanmoins, le rapport de module entre l'inclusion et le milieu environnant reste cohérent avec celui obtenu dans les autres configurations (table 4.8). Pour ce cas plus particulier, la surface sur laquelle s'applique la force est plus petite que la surface du milieu, et il a pu être noté que des valeurs de module proches des valeurs réelles pouvaient être obtenues si une surface légèrement plus grande était considérée. Pour finir, il est possible de remarquer que les conditions aux limites imposées lors de la simulation ont une forte influence sur les images de déformation (fig. 4.14,b-d), entrainant dans certaines zones des contrastes importants. Des variations du module d'Young localisées dans ces mêmes zones sont visibles sur les cartes reconstruites.

<table><caption>TABLE 4.8 – Résultats sur les valeurs de module d'Young estimées : moyenne et écart type [kPa], rapport R (4.1) et erreur quadratique moyenne (RMSE, [kPa]) en fonction des conditions aux limites du milieu simulé (fig. 4.14). La ligne en gris correspond aux valeurs estimées en prenant en compte l'ensemble de milieu environnant (fig. 4.1-b).</caption><thead><tr><th></th><th>E<sub>b</sub> estimé</th><th>E<sub>i</sub> estimé</th><th>Rapport R</th><th>RMSE milieu environnant</th><th>RMSE inclusion</th></tr></thead><tbody><tr><td>a1</td><td>10.03 ± 0.07</td><td>22.05 ± 0.59</td><td>2.20</td><td>0.07</td><td>7.97</td></tr><tr><td>a2</td><td>10.28 ± 0.35</td><td>22.28 ± 0.84</td><td>2.17</td><td>0.45</td><td>7.77</td></tr><tr><td rowspan="2">a3</td><td>33.27 ± 2.72</td><td>66.51 ± 4.87</td><td>2.00</td><td>23.43</td><td>36.83</td></tr><tr><td>30.30 ± 12.88</td><td>66.51 ± 4.87</td><td>2.20</td><td>24.04</td><td>36.83</td></tr></tbody></table>
<figure><img src="image_52.png" /><figcaption>FIGURE 4.14 – Cartes de module d'Young reconstruites (a) et composantes axiales (b), en cisaille-ment (c) et latérales (d) des tenseurs de déformation 2D. Les cartes de module sont reconstruites pour un même milieu 3D soumis à des conditions aux limites différentes. Ces conditions sont : la surface supérieure est soumise à un chargement uniforme et vertical et celle inférieure est contrainte verticalement (gauche), la surface supérieure est soumise à un chargement uniforme et vertical et celle inférieure est fixe (centre), la force est appliquée sur une zone restreinte et la surface inférieure est contrainte verticalement (droite).</figcaption></figure>

Pour conclure, dans l'ensemble des cas étudiés, les cartes de module d'Young reconstruites permettent clairement de distinguer des régions aux propriétés mécaniques, tailles ou positions différentes. Cependant, les conditions de contraintes planes n'étant pas satisfaites pour les simula-tions 3D, une sur/sous-estimation du module d'Young dans les inclusions entraine une diminution du contraste élastique.
# IV Données expérimentales

La méthode a ensuite été évaluée à l'aide de données expérimentales acquises sur des fantômes et sur des tissus mammoires *in vivo*. Les données expérimentales sont utilisées pour déterminer si la méthode est toujours capable de produire une image des propriétés élastiques des milieux dans des conditions d'acquisition réelles. Les tests sur fantômes constituent ainsi une étape intermédiaire pour l'évaluation de la méthode, car ces objets sont moins complexes que les tissus biologiques. De plus, pour ceux examinés dans ces travaux, des informations sur leurs propriétés mécaniques sont données par le fabricant, permettant une analyse plus approfondie des cartes reconstruites. Pour les tissus biologiques *in vivo*, aucune information mécanique n'est évidemment disponible. Cependant, ils révèlent la complexité de l'élastographie mammaire, application pour laquelle notre méthode a été développée.

L'acquisition des données et les milieux utilisés sont présentés ci-dessous, ainsi que les résultats obtenus.

## IV.1 Milieux utilisés

Nous retrouvons dans ce chapitre des milieux décrits précédemment, lors de l'imagerie du teneur de déformation. Plus précisément, des reconstructions de la distribution du module d'Young vont être effectuées pour deux fantômes d'élastographie mammaire (CIRS modèles 059). Pour rappel, ces milieux imitent un sein de femme en position allongée, dans lequel se trouvent plusieurs inclusions sphériques de module d'Young environ 3 fois égal à celui du milieu environnant. Une information complémentaire va cependant être apportée dans ce chapitre. Si pour l'un des fantômes (fantôme 1), ce rapport de 3 entre les modules est l'unique information fournie par le constructeur, pour le second fantôme (fantôme 2), les mesures effectuées par le fabricant ont indiqué un module d'Young de 43.3 kPa pour les inclusions et de 13 kPa pour le milieu environnant. Trois régions ont été imagées, l'une contenant une seule inclusion (cas #1, fantôme 1), et les deux autres, deux inclusions superposées (cas #2, fantôme 2) ou positionnées côte à côte (cas #3, fantôme 2).

Trois cas de lésions mammoires vont également être testés (cas #4, cas #5, cas #6). Pour rappel, ces cas ont été extraits d'une base de données archivée, résultant d'une étude d'élastographie mammaire réalisée aux Hospices Civils de Lyon (Collaboration, Dr A. Coulon). Des informations détaillées sur cette étude sont disponibles dans l'article suivant [4]. Deux cas de carcinome canalaire infiltrant (CCI) (cas #4 et #6) et un cas de fibroadénome (lésion bénigne, cas #5) seront présentés.

Les données ultrasonores radiofréquences à partir desquelles les champs de déplacement sont estimés (puis régularisés, Chap. 2) ont été acquises lors d'expériences typiques d'élastographie quasi-statique. Il doit être noté que lors de ces expériences, la valeur de la force appliquée aux milieux n'était pas disponible. Cette valeur a donc été fixée arbitrairement et seules des cartes relatives du module d'Young pourront être obtenues.

## IV.2 Résultats obtenus sur fantômes

Les résultats obtenus sous les hypothèses de contraintes planes, ainsi que les images du teneur de déformation utilisées, sont d'abord présentés pour les données fantômes cas #1 (fig. 4.15), cas #2 (fig. 4.16) et cas #3 (fig. 4.17). Dans tous les cas examinés, les inclusions sont clairement visibles sur les cartes relatives de module d'Young reconstruites, et elles apparaissent plus rigides que le milieu environnant, comme cela est attendu. Il peut être noté que les cartes reconstruites ont toutes un aspect similaire bien que les milieux et la force appliquée soient différents. Pour le cas #1, l'inclusion et le milieu environnant apparaissent relativement homogènes (fig. 4.15-d).
Néanmoins, une "zone rouge" (donc supposément plus rigide) peut être observée dans la partie supérieure de l'image. Cette zone est de petite taille et n'affecte pas l'interprétation du résultat. Pour le cas #2 (fig. 4.16-d1) et le cas #3 (fig. 4.17-d1), cette zone est légèrement plus étendue (en particulier pour le cas #3) avec des valeurs plus élevées, et des variations verticales plus prononcées du module d'Young relatif peuvent être observées, conduisant à des régions d'aspect moins homogènes. Des zones de rigidité plus élevée ont pu être observées précédemment sur les images obtenues à partir de simulations 3D lorsque les conditions aux limites variaient, et semblent correspondre à des zones où le contraste de déformation est important. Pour le cas #2 et le cas #3, des cartes relatives du module d'Young ont été reconstruites en supprimant au préalable les premières lignes des champs de déformation (fig. 4.16-4.17-d2.). Les résultats obtenus sont alors de qualité similaire à celui obtenu pour le cas #1. Pour le cas #2, l'inclusion supérieure apparaît plus petite sur les images, bien que qu'elle semble être de même taille que l'inclusion basse, lorsqu'on examine le milieu à l'échographie. Cependant, lors de l'acquisition, ces deux inclusions n'étaient pas exactement situées dans le même plan, et explique la plus petite taille de l'inclusion haute. Enfin, cette inclusion apparaît également moins rigide que celle inférieure, bien que le rapport des modules est supposé être le même. La encore, le fait que les deux inclusions n'étaient pas tout à fait dans le même plan lors de l'acquisition des données peut avoir contribué aux différences observées.

Ces observations sont confirmées par le rapport de module d'Young R estimé entre la/les inclusions et le milieu environnant (tab. 4.9). Pour rappel, le module d'Young relatif du milieu environnant est soit estimé dans une ROI circulaire de même taille que celle de l'inclusion (configuration 1), soit sur toute l'image en excluant l'inclusion et les bords (configuration 2). Il faut noter pour la première configuration, que les ROI du milieu environnant ont été sélectionnées dans des zones visuellement homogènes. En raison d'une variabilité des valeurs dans le milieu environnant plus fortes que pour les résultats de simulations, les ROI doivent être soigneusement sélectionnées. Pour toutes les inclusions, le rapport R est compris entre 2 et 3, et plus exactement dans la gamme [2.04 – 2.92], sauf pour l'inclusion haute du cas #2, où ce rapport est inférieur, compris entre 1.61 et 2.21. Cette diminution du rapport des modules par rapport aux valeurs fournies par le fabricant rejoint les observations faites précédemment. Pour finir, si supprimer les premières lignes des images permet d'améliorer l'aspect visuel des cartes reconstruites, le rapport élastique estimé peut être légèrement modifié (plus faible dans le cas #2 et plus élevé dans le cas #3).
<figure><img src="image_53.png" /></figure>

FIGURE 4.15 – Résultats expérimentaux cas #1 - obtenus avec un fantôme CIRS, modèle 059 (fantôme 1) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young reconstruite.

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

FIGURE 4.16 – Résultats expérimentaux cas #2 - obtenus avec un fantôme CIRS, modèle 059 (fantôme 2) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young reconstruite. La carte (d1) est reconstruite en utilisant l'intégralité des images du tenseur des déformations et (d2) en supprimant les premières lignes de ces images.
<figure><img src="image_55.png" /><figcaption>FIGURE 4.17 – Résultats expérimentaux cas #3 - obtenus avec un fantôme CIRS, modèle 059 (fantôme 2) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young reconstruite. La carte (d1) est reconstruite en utilisant l'intégralité des images du tenseur des déformations et (d2) en supprimant les premières lignes de ces images.</figcaption></figure>

<table><caption>TABLE 4.9 – Rapports élastiques *R* estimés pour les cas de fantômes à l'aide des deux configurations de sélection des régions. La configuration 1 utilise deux ROI circulaires de même taille et la configuration 2 sélectionne le milieu environnant comme étant l'ensemble de l'image moins la/les inclusions et les bords. Les lignes grises correspondent aux résultats obtenus en supprimant les premières lignes des images de déformation, et les autres en utilisant l'intégralité des images.</caption><thead><tr><th rowspan="2"></th><th colspan="2">Rapport élastique R</th></tr><tr><th>configuration 1</th><th>configuration 2</th></tr></thead><tbody><tr><td>Cas #1</td><td>2.42</td><td>2.46</td></tr><tr><td rowspan="2">Cas #2 - inc 1</td><td>2.21</td><td>1.75</td></tr><tr><td>1.95</td><td>1.61</td></tr><tr><td rowspan="2">Cas #2 - inc 2</td><td>2.92</td><td>2.86</td></tr><tr><td>2.53</td><td>2.55</td></tr><tr><td rowspan="2">Cas #3 - inc 1</td><td>2.31</td><td>2.04</td></tr><tr><td>2.50</td><td>2.15</td></tr><tr><td rowspan="2">Cas #3 - inc 2</td><td>2.30</td><td>2.33</td></tr><tr><td>2.46</td><td>2.37</td></tr></tbody></table>
IV.3 Résultats sur tissus mammaires *in vivo*

Pour finir, les résultats obtenus à partir des données *in vivo*, sous les hypothèses de contraintes planes, sont présentés. Pour tous les cas, la présence dans ces images d'une région clairement plus rigide que les tissus environnants concorde avec la position de la lésion. Ces premiers résultats sont compatibles avec les lésions analysées puisqu'il s'agit ici de deux carcinomes canalaires infiltrants, cancers connus pour être généralement beaucoup plus rigides que les tissus sains [5], ainsi que d'un fibroadénome pour lequel il était précisé qu'il était palpable. En comparaison avec les résultats présentés précédemment, les cartes relatives de module d'Young reconstruites semblent moins homogènes. Cependant, cela peut être attendu au vu de la plus grande hétérogénéité et complexité des tissus biologiques. Le rapport des modules *R* entre la lésion et le milieu environnant est estimé autour de 6.6 pour le cas #4, 4.5 pour le cas #5 et 3.6 pour le cas #6 (tab. 4.10). Si ce rapport est plus grand pour le cas #4, il est proche pour les cas #5 et #6, et ne permet pas de différencier les deux lésions, alors que l'une est un cancer et l'autre un fibroadénome. Des valeurs proches de ce rapport ne sont pas impossibles, d'autant plus que le fibroadénome imagé était palpable, et les diverses études cliniques publiées montrent qu'il existe toujours des lésions qui ne peuvent être séparées au vu des critères analysés [5]. Nous savons aussi que les hypothèses introduites au niveau de la méthode, notamment l'hypothèse de contraintes planes, conduisent à des erreurs lors de la reconstruction des images de modules, pouvant entraver la séparation entre lésions bénignes et lésions malignes. Cependant, les valeurs obtenues dans notre étude sont du même ordre que celles relevées dans des études cliniques utilisant l'élastographie par ondes de cisaillement. En effet, dans une étude sur 939 masses où ce rapport des modules avait notamment été calculé [5], sa valeur médiane était de 7.8 pour les cancers et de 2.0 pour les lésions bénignes. Une analyse plus poussée des résultats est néanmoins difficile, en raison de la grande complexité des milieux biologiques et du manque crucial d'informations sur les propriétés mécaniques des tissus examinés.

<figure><img src="image_56.png" /><figcaption>FIGURE 4.18 – Résultats obtenus sur des tissus mammaires *in vivo*, dans le cas d'un carcinome canalaire infiltrant (cas #4) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young.</figcaption></figure>
<figure><img src="image_57.png" /></figure>

FIGURE 4.19 – Résultats obtenus sur des tissus mammaires *in vivo*, dans le cas d'un fibroadénome (cas #5) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young.

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

FIGURE 4.20 – Résultats obtenus sur des tissus mammaires *in vivo*, dans le cas d'un carcinome canalaire infiltrant (cas #6) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale, (d) carte relative de module d'Young.

# V Discussion

Afin de reconstruire des cartes des propriétés élastiques (module d'Young) des tissus biologiques en élastographie ultrasonore quasi-statique, nous avons proposé une approche basée sur le principe des travaux virtuels, utilisant l'ensemble des composantes du tenseur de déformation 2D. Le tenseur de déformation a été obtenu grâce à l'introduction d'une méthode de régularisation (chap. 2), indispensable à l'obtention de données exploitables. En effet, sans régularisation, les champs de déformation latérale et en cisaillement sont dominés par le bruit et ils ne peuvent permettre une reconstruction correcte du module d'Young. Une illustration en est donnée avec
<table><thead><tr><th rowspan="2"></th><th colspan="2">Rapport élastique R</th></tr><tr><th>configuration 1</th><th>configuration 2</th></tr></thead><tbody><tr><td>Cas #4</td><td>7.19</td><td>6.04</td></tr><tr><td>Cas #5</td><td>5.12</td><td>4.02</td></tr><tr><td>Cas #6</td><td>3.53</td><td>3.84</td></tr></tbody></table>

TABLE 4.10 – Rapports élastiques *R* estimés pour les tissus mammaires *in vivo* à l'aide des deux configurations de sélection des régions. La configuration 1 utilise deux ROI circulaires de même taille et la configuration 2 sélectionne le milieu environnant comme étant l'ensemble de l'image moins la/les inclusions et les bords.

le fantôme *CIRS* modèle 059 (cas #1), en figure 4.21. Ici, la carte relative du module d'Young reconstruite à partir de champs de déformation non régularisés est bruitée et la présence de l'inclusion plus rigide que le milieu environnant n'est pas révélée.

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

FIGURE 4.21 – Résultat expérimental obtenu à partir de champs de déformation non régularisés, pour un fantôme *CIRS* modèle 059 (cas #1) : (a) déformation axiale, (b) déformation en cisaillement, (c) déformation latérale et (d) carte relative de module d'Young reconstruite sous les hypothèses de contraintes planes.

À l'inverse, les différents résultats présentés dans ce chapitre ont montré des cartes de module d'Young interprétables où des régions aux propriétés mécaniques différentes au sein d'un milieu étaient observables, que ce soit pour les simulations, les fantômes ou les tissus mammaires. Sur les cartes reconstruites, la présence d'inclusions/lésions aux propriétés élastiques distinctes sont clairement visibles. Pour les données expérimentales sur fantômes, les rapports des modules d'Young estimés sont inférieurs mais restent cohérents avec les valeurs fournies par le fabricant, et la diminution du contraste élastique est en accord avec les résultats de simulations. Cette diminution est liée à l'approche 2D utilisée dans notre méthode (avec l'hypothèse des contraintes planes) alors que le problème est tridimensionnel. Pour les résultats sur les tissus biologiques *in vivo*, ils doivent être interprétés avec précaution. Ils laissent cependant apparaître une zone plus rigide qui correspond à la région de la lésion. De plus, les valeurs de rapports obtenus sont dans la gamme de valeurs que l'on peut trouver dans la littérature. Pour finir, les divers résultats présentés, en particulier ceux issus de données *in vivo*, ont montré la robustesse de la méthode face à des milieux divers et complexes. En ce qui concerne la précision de la méthode, les résultats obtenus à partir de simulations, pour lesquels des conditions de contraintes planes ou de déformations planes sont respectées, ont montré une estimation précise de la distribution spatiale du module d'Young, avec une erreur quadratique moyenne (RMSE) très faible. En revanche, dès lors que ces conditions ne sont plus satisfaites, une diminution du contraste élastique peut être observée, avec une diminution plus prononcée sous les hypothèses de déformations planes (fig. 4.5). En élastographie ultrasonore quasi-statique, aucune de ces hypothèses ne permet de modéliser parfaitement le comportement des tissus biologiques, et les deux sont couramment
utilisées [6-9]. Au vu de nos résultats, le choix a été fait d'utiliser, dans nos travaux, principalement l'hypothèse des contraintes planes. Il n'est cependant pas exclu que dans certaines situations, l'hypothèse des déformations planes se montre plus adaptée. Afin de s'affranchir de ces hypothèses, et grâce aux développements actuels des systèmes d'acquisition ultrasonore 3D, de futurs travaux porteront sur l'extension de la méthode proposée à la troisième dimension. Cette extension apparaît primordiale afin d'obtenir à la fois des valeurs de module (en kPa) et des contrastes élastiques plus proches des valeurs réelles, avec les données expérimentales. Enfin, une limite peut être soulevée concernant certains cas, où des zones apparaissaient plus rigides sur les cartes reconstruites sans que les propriétés mécaniques du milieu ne soient différentes. Ces zones semblent être dues aux conditions aux limites qui créent un fort contraste de déformation dans certaines régions, impactant ensuite les cartes reconstruites. Toutefois, il peut être noté que les régions concernées se concentrent essentiellement sur les bords des images et peuvent donc être préalablement supprimées.

Actuellement, à partir des données expérimentales, seules des cartes du module d'Young relatif peuvent être reconstruites. En effet, dans la pratique de l'élastographie ultrasonore quasi-statique, le mouvement de compression/décompression est réalisé manuellement par le/la radiologue à l'aide de la sonde échographique. Par conséquent, la force appliquée sur le milieu est inconnue, empêchant ainsi une estimation quantitative des propriétés mécaniques des tissus. Dans nos travaux, la force a été fixée à une valeur arbitraire, mais des développements en cours, portant sur un système permettant d'acquérir la valeur de la force appliquée aux tissus en synchronisation avec les données ultrasonores, pourraient prochainement permettre d'accéder à des valeurs de module d'Young.

Il peut être souligné ici qu'aucun filtrage n'a été appliqué aux cartes (relatives) de module d'Young présentées. En revanche, nous rappelons que lors de l'estimation d'une valeur de module, une prise en compte pondérée du voisinage est opérée, ce qui peut produire un lissage plus ou moins important selon la taille du voisinage et les coefficients de pondération choisis. Dans nos travaux, ces coefficients sont déterminés par l'utilisation d'un masque gaussien 3 × 3. Cependant, si dans les travaux futurs, des tests supplémentaires révèlent la nécessité d'un lissage plus important, un masque plus large pourra être utilisé. De même, un autre type de masque pourra être considéré en fonction des besoins.

Concernant le temps de calcul, la reconstruction d'une carte telle que nous en avons montrée dans ce travail (composée environ de 75 × 75 pixels) prend en moyenne 4 minutes. Il est à noter que pour cette première version de la méthode, les aspects temps réel n'ont pas été pris en considération. Certains choix ont malgré tout été faits afin de limiter le temps de calcul. La taille du maillage est l'un paramètre qui intervient dans le temps de calcul. En effet, un maillage large comporte peu d'éléments ce qui limite le nombre d'inconnues à déterminer et donc le temps de calcul. Cependant, comme nous l'avons évoqué dans le chapitre précédent, ce choix est aussi un compromis à faire avec la résolution des images reconstruites. Ces questions seront intéressantes à prendre en compte dans de futurs travaux.

Pour finir, dans ces travaux, le coefficient de Poisson a été défini comme étant constant et connu pour tous les milieux étudiés (égal à 0.49), réduisant ainsi le nombre d'inconnues. Par conséquent, seule une cartographie du module d'Young reste à reconstruite. Une étude réalisée avec des simulations de milieux 3D possédant différentes valeurs de coefficients de Poisson, a montré une influence limitée de ce paramètre. En effet, seule une légère diminution des modules d'Young estimés a pu être observée lorsque le coefficient de Poisson du milieu s'éloigne de manière importante de la valeur fixée (par exemple 0.1 vs 0.49, fig. 4.6). Toutefois, certaines applications pourraient bénéficier d'une estimation de la distribution spatiale du coefficient de Poisson. Pour
cela, la méthode développée peut être adaptée afin de réintégrer ce terme comme étant une incon-
nue à estimer. Cela suppose alors de résoudre un problème contenant deux fois plus d'inconnues.
Au delà, d'une estimation du module d'Young et du coefficient de Poisson, des modèles plus
complexes, permettant de modéliser le comportement mécanique des milieux étudiés, pourraient
aussi être envisagés afin de se rapprocher des propriétés des tissus biologiques. En effet, le choix
du modèle modifie simplement l'écriture de la relation contrainte-déformation utilisée dans notre
méthode. Néanmoins, cela soulève encore une fois la question du nombre d'inconnues, ce qui
pourra être le sujet de travaux futurs.

VI Conclusion

Dans ce chapitre, l'imagerie des propriétés élastiques des tissus biologiques a été étudiée grâce à l'introduction d'une méthode de reconstruction spécifique, basée sur l'utilisation des champs virtuels et l'estimation de l'ensemble des composantes du tenseur de déformation 2D. Nos résultats montrent la capacité de notre méthode à reconstruire des cartes du module d'Young des milieux, notamment pour les données expérimentales. Plus particulièrement, la méthode proposée montre des résultats prometteurs pour la reconstruction de cartes relatives du module d'Young dans les tissus mammaires in vivo. Si ces résultats ne sont pas parfaits, les lésions peuvent être détectées visuellement et les rapports des modules estimés entre ces lésions et le milieu environnant sont cohérents avec les données de la littérature.
Bibliographie

[1] *Comsol Multiphysics®*. Comsol AB, Stockholm, Sweden, 2019. URL http://www.comsol.com. III

[2] A.-L. Duroy, O. Basset, and E. Brusseau. Virtual fields based-method for reconstructing the elastic modulus in quasi-static ultrasound elastography. In *IEEE International Ultrasonics Symposium*, 2022. III.2.2

[3] A.-L. Duroy, O. Basset, and E. Brusseau. Elastic modulus imaging for breast application using a virtual fields based-method in quasi-static ultrasound elastography. IEEE TUFFC, 2023. III.2.2

[4] E. Brusseau, V. Detti, A. Coulon, E. Maissiat, N. Boublay, Y. Berthezène, J. Fromageau, N. Bush, and J. Bamber. In Vivo Response to Compression of 35 Breast Lesions Observed with a Two-Dimensional Locally Regularized Strain Estimation Method. Ultrasound in Medicine & Biology, 40:300-312, 2014. IV.1

[5] W. A. Berg, D. O. Cosgrove, C. J. Doré, F. K. W. Schäfer, W. E. Svensson, R. J. Hooley, R. Ohlinger, E. B. Mendelson, C. Balu-Maestro, M. Locatelli, C. Tourasse, B. C. Cavanaugh, V. Juhan, A. T. Stavros, A. Tardivon, J. Gay, J.-P. Henry, C. Cohen-Bacrie, and For the BE1 Investigators. Shear-wave Elastography Improves the Specificity of Breast US : The BE1 Multinational Study of 939 Masses. Radiology, 262(2):435-449, 2012. IV.3

[6] P. E. Barbone and J. C. Bamber. Quantitative elasticity imaging : what can and cannot be inferred from strain images. Physics in Medicine and Biology, 47(12):2147-2164, 2002. V

[7] M. M. Doyley. Model-based elastography : a survey of approaches to the inverse elasticity problem. Physics in Medicine and Biology, 57:R35-R73, 2012.

[8] O. Babaniyi, A. Oberai, and P. Barbone. Recovering vector displacement estimates in quasi-static elastography using sparse relaxation of the momentum equation. Inverse Problems in Science and Engineering, 25(3):326-362, 2017.

[9] D. T. Seidl, A. A. Oberai, and P. E. Barbone. The Coupled Adjoint-State Equation in forward and inverse linear elasticity : Incompressible plane stress. Computer Methods in Applied Mechanics and Engineering, 357:1125-1138, 2019. V
<figure><img src="image_60.png" /></figure>
Conclusion générale et perspectives

Les travaux présentés dans ce manuscrit de thèse ont porté sur l’élastographie ultrasonore quasi-statique. Aujourd’hui, cette technique d’imagerie est implantée sur des échographes cliniques et les critères utilisés pour le diagnostic sont issus principalement de la déformation axiale. L’une des raisons à l’utilisation prédominante de cette composante de la déformation est due à l’estimation plus précise du déplacement axial comparé au déplacement latéral. Afin d’analyser les propriétés élastiques des tissus biologiques, deux axes de recherche ont été développés dans ces travaux :

1. l'imagerie du tenseur de déformation 2D

2. la reconstruction de cartographies du module d'Young

Afin d'accéder à une imagerie du tenseur de déformation 2D, une méthode de régularisation basée sur la propriété de quasi-incompressibilité des tissus biologiques a été développée, permettant un lissage des champs de déplacement. Plus exactement, un problème d'optimisation a été formulé comportant, outre les termes d'attache aux données, des termes de pénalisation des fortes variations sur les composantes du déplacement ainsi que le terme de divergence du déplacement. Des tests sur des fantômes et sur des tissus mammaires in vivo ont montré la forte contribution de cette méthode face à des données complexes et bruitées, notamment pour la composante latérale du déplacement. Il en résulte des images de la déformation latérale et de cisaillement qui deviennent interprétables, et les inclusions/lésions deviennent ainsi aussi clairement visibles sur les images de déformation latérale qu'elles l'étaient sur les images de déformation axiale, ce qui n'était pas le cas avant la régularisation. Ces résultats ont été confirmés par une augmentation significative des valeurs de $CNR_{lateral}$. En plus d'une imagerie complète du tenseur de déformation 2D, des élastogrammes de rotation pour lesquels des motifs de rotation étaient clairement détectables autour des inclusions/lésions ont aussi été obtenus. Enfin, il est important de soulagner que, bien que notre méthode soit une approche 2D, les champs de déformation régularisés sont cohérents avec des milieux se déformant en 3D. Le tenseur de déformation ainsi obtenu a ensuite été utilisé pour la reconstruction des propriétés mécaniques des tissus, ce qui a fait l'objet du deuxième axe de ces travaux de recherche.

Afin de reconstruire des cartographies des propriétés mécaniques des tissus biologiques, une méthode basée sur les champs virtuels a été développée. Plus exactement, pour ces premiers travaux, nous nous sommes intéressés à la reconstruction du module d'Young pour des milieux élastiques, linéaires, isotropes et quasi-incompressibles (c.-à-d. avec un coefficient de Poisson égal à 0.49). De plus, pour pallier le manque de données 3D, les hypothèses de déformations et de contraintes planes ont été étudiées. Ces hypothèses sont souvent utilisées en élastographie, bien qu'en pratique, aucune d'entre elles ne permette de modéliser parfaitement le comportement des tissus biologiques. De premiers tests réalisés à l'aide de simulations numériques ont montré qu'une reconstruction précise de la distribution spatiale du module d'Young était obtenue lorsque les hypothèses introduites dans la méthode étaient vérifiées. Cependant, comme cela pouvait être attendu, une diminution du contraste élastique (estimé entre la/les inclusion(s) et le milieu envi-
ronnant) a été observée quand les hypothèses de contraintes ou de déformations planes n'étaient plus satisfaites, et ce, plus fortement dans le deuxième cas. Pour cette raison, dans ces travaux, ce sont les hypothèses de contraintes planes qui ont été principalement considérées. Des tests sur des fantômes ont par la suite permis de reconstruire des cartes relatives de module d'Young sur lesquelles les inclusions sont clairement visibles et où les différentes régions apparaissent relativement homogènes. Enfin, des résultats prometteurs ont été montrés pour des tissus mammaires *in vivo*, pour lesquels une lésion plus rigide que les tissus environnants pouvait être détectée. Pour ces premiers travaux, seules des cartes relatives du module d'Young ont pu être reconstruites dans le cas des fantômes et des tissus biologiques, car la valeur de la force appliquée sur le milieu est inconnue. Cependant, des travaux en cours portant sur le développement d'un système permettant de mesurer la force en parallèle de l'imagerie, pourraient permettre d'obtenir prochainement des cartes quantitatives (en kPa). Ce système est destiné à enregistrer les données issues d'un capteur de force en synchronisation avec l'acquisition des images ultrasonores radiofréquences. De nombreux tests devront être effectués, avec notamment pour les fantômes, une comparaison des modules estimés avec des modules de référence fournis par le constructeur ou issus d'essais mécaniques. Pour les données d'examens sur patients, les propriétés mécaniques des tissus biologiques étant inconnues, seule une comparaison avec une autre technique d'élastographie pourra être effectuée, afin d'obtenir une première appréciation des performances du système développé. Enfin, la méthode proposée ici a l'avantage de permettre une grande liberté d'utilisation, et de futurs travaux pourront porter sur la reconstruction d'autres paramètres mécaniques.

Pour finir, les développements réalisés se sont appuyés sur l'utilisation d'images et non pas de volumes ultrasonores acquis au cours du temps. Avec les recherches actuelles portant sur le développement de systèmes d'acquisition de données volumiques, il pourra être envisagé pour de futurs travaux, d'étendre les méthodes développées à la troisième dimension. En effet, pour la méthode de régularisation, l'introduction de la composante en élévation du déplacement est immédiate et aucune difficulté ne devrait être rencontrée. Quant à la méthode de reconstruction, cette dernière bénéficiera de l'estimation du tenseur de déformation 3D, mais surtout, l'utilisation de données volumiques permettra de s'affranchir des hypothèses de déformations et de contraintes planes, qui comme nous l'avons vu avec les résultats de simulations numériques, sont une source majeure de perte de contraste dans les cartes de module reconstruites.
Publications et communications

Revues avec comité de lecture

* **A-L. Duroy**, V. Detti, A. Coulon, O. Basset, E. Brusseau. Regularization-based 2D tissue strain tensor imaging in quasi-static ultrasound elastography. *Ultrasonic Imaging*, 2023 - accepté

* **A-L. Duroy**, O. Basset, E. Brusseau. Elastic modulus imaging for breast application using a Virtual fields based-method in quasi-static ultrasound elastography. *IEEE TUFFC*, 2023 - soumis

Proceedings

* **A-L. Duroy**, V. Detti, A. Coulon, O. Basset and E. Brusseau. 2D tissue strain tensor imaging in quasi-static ultrasound elastography. *43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)*, pp. 2847-2851, 2021
https://ieeexplore.ieee.org/abstract/document/9630570

* **A-L. Duroy**, O. Basset and E. Brusseau. Virtual fields based-method for reconstructing the elastic modulus in quasi-static ultrasound elastography. *IEEE International Ultrasonics Symposium (IUS)*, 2022
https://hal.archives-ouvertes.fr/hal-03793187

* E. Brusseau, **A.L. Duroy**, O. Basset, A. Coulon, and V. Detti. Breast lesion examination with ultrasound elastography, *International Conference in Gynaecology and Obstetrics*, 2023

Communications

* Présentation (virtuel) pour *43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)*, 2021

* Poster (présentiel) pour *IEEE International Ultrasonics Symposium (IUS)*, 2022

* Intervention auprès de lycéennes dans le cadre de la journée de la femme 2021 et 2022 :
  *Sciences, un métier de femmes!*
https://www.femmesetsciences.fr/news/sciences%2C-un-m%C3%A9tier-de-femmes-!-2022