# Une application de la théorie du contrôle a un problème inverse de sismique

Guy Chavent, Patrick Lailly, Alain Bamberger

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

## To cite this version:

Guy Chavent, Patrick Lailly, Alain Bamberger. Une application de la théorie du contrôle a un problème inverse de sismique. [Rapport de recherche] IRIA-RR-172, IRIA. 1976, pp.55. hal-04716340

HAL Id: hal-04716340

https://inria.hal.science/hal-04716340v1

Submitted on 1 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
IRIA

laboria

Institut de Recherche
d'Informatique
et d'Automatique

Domaine de Voluceau
Rocquencourt
B.P. 5 78150 - Le Chesnay
France
Tél.: 954 90 20

laboratoire de recherche
en informatique
et automatique

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

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

# UNE APPLICATION DE LA THÉORIE DU CONTRÔLE A UN PROBLÈME INVERSE DE SISMIQUE

Guy CHAVENT
Patrick LAILLY
Alain BAMBERGER

Rapport de Recherche N° 172

Juin 1976
UNE APPLICATION DE LA THÉORIE DU CONTROLE
A UN PROBLEME INVERSE DE SISMIQUE

G. Chavent*, P. Lailly*, A. Bamberger**

Résumé :

Une méthode nouvelle de détermination des impédances acoustiques à partir des enregistrements réalisés en sismique réflexion est exposée dans ce rapport. Elle consiste essentiellement dans une application de la théorie du contrôle à la détermination des paramètres inconnus d'une équation aux dérivées partielles utilisant une méthode de gradient. On étudie d'abord comment se pose le problème. On met en évidence quelques propriétés de l'équation des ondes et leurs conséquences pour le problème étudié. On propose une résolution numérique et on montre des résultats obtenus par traitement numérique.

Abstract :

This paper deals with a new method for the determination of acoustical impedance from recordings obtained by reflection seismic. The aim of the method is to apply optimal control theory for the determination of unknown coefficients in PDE. After having defined precisely the inverse problem, we bring in light some properties of the wave equation and point out their consequences for the seismic inverse problem. A computational algorithm is given, and illustrated by numerical results.

<footnote>* IRIA-LABORIA - Rocquencourt.</footnote>

<footnote>** Centre de Mathématiques Appliquées. Ecole Polytechnique - Massy.</footnote>
# PLAN

## INTRODUCTION

### 1 - POSITION DU PROBLEME

1.1 - Notations

1.2 - Equation retenue pour représenter le phénomène

1.3 - Transformation de l'équation

1.4 - Position du problème : ses difficultés

1.4.1 - Le problème inverse

1.4.2 - Formulation du problème d'inversion comme un problème de minimisation

1.4.3 - Les difficultés du problème

### 2 - QUELQUES PROPRIETES DE L'EQUATION DES ONDES ET LEURS CONSEQUENCES POUR LE PROBLEME INVERSE

2.1 - Un résultat de convergence

2.1.1 - Enoncé du résultat

2.1.2 - Exemples d'application

2.1.3 - Simulation numérique d'un exemple

2.1.4 - Conséquences sur le problème d'identification

2.2 - Un résultat de continuité

2.3 - Un résultat de comportement lorsque t tend vers l'infini

2.3.1 - Enoncé du résultat

2.3.2 - Illustration numérique : conséquences sur le problème d'identification

### 3 - RESOLUTION NUMERIQUE DU PROBLEME

3.1 - Résolution de l'équation d'état

3.1.1 - Discrétisation en espace

3.1.2 - Discrétisation en temps

3.1.3 - Conditions de stabilité

3.1.4 - Commentaires sur le choix de la discrétisation
## 3.2 - Résolution du problème d'identification

3.2.1 - Méthode d'optimisation

3.2.2 - Calcul du gradient de J

3.2.3 - Mise en oeuvre

3.2.4 - Organigramme

# 4 - RESULTATS NUMERIQUES

4.1 - Choix des différents paramètres et description des exemples traités

4.1.1. - Modèle et algorithme retenus

4.1.2 - Description des exemples étudiés

4.1.3 - Essais effectués

4.2 - Essais à épaisseur de couches connue

4.2.1 - Exposé des résultats

4.2.2 - Commentaires sur le cheminement de l'algorithme

4.2.3 - Temps de calculs - vitesse de convergence

4.3 - Essais à épaisseur de couches inconnue : stabilité de l'algorithme par rapport à la finesse de résolution

4.3.1 - Exposé des résultats

4.3.2 - Conclusions sur la stabilité par rapport à la discrétisation des impédances

4.4 - Conclusions
INTRODUCTION

En sismique réflexion, on enregistre des mesures d'amplitude ou de pression et l'information que l'on peut déduire de ces enregistrements concerne principalement la structure des interfaces géologiques et l'impédance acoustique caractérisant les différentes couches. Les structures sont mises en évidence par des traitements appropriés mais leur détermination fait inter-venir des distances non négligeables. Les impédances acoustiques peuvent au contraire, moyennant certaines précautions, être évaluées à la verticale d'un point donné de la surface du sol. Il s'agit là, en fait, des deux problèmes fondamentaux de la géophysique pétrolière. En particulier la solution du deuxième problème consisterait à caractériser la nature des terrains que traversent les ondes acoustiques.

Les solutions les plus utilisées dans l'état actuel du traitement sismique consistent à déterminer les vitesses moyennes de propagation d'où l'on déduit, connaissant les temps de parcours des ondes le long des rayons, les vitesses de propagation à l'intérieur de chaque couche. Ces méthodes fournissent des résultats dont l'interprétation est délicate. Elles supposent tout d'abord, que les interfaces n'ont pas de pendage. Elles sont ensuite basées sur des critères de meilleure mise en phase des arrivées d'énergie plutôt que sur des critères de meilleure restitution des interfaces à leur vraie position dans l'espace. L'erreur commise sur les résultats est enfin d'autant plus grande que l'épaisseur des couches est faible.

Une autre tentative de détermination des vitesses a été proposée par G. Kunetz [6]. Le procédé reproduit le calcul inverse des sismogrammes synthétiques à une dimension. On suppose que certains traitements effectués sur les enregistrements réels permettent de définir l'enregistrement que l'on aurait eu en un point donné de la surface si on avait pu générer des ondes planes en ce point et si un récepteur avait été placé en ce même point. On peut alors théoriquement calculer les impédances acoustiques qui ont produit l'enregistrement, à condition que ce dernier obéisse à certaines contraintes bien définies. La méthode, très séduisante, s'applique très bien dans le cas de sismogrammes obtenus par calcul. L'application à des cas réels se heurte à des difficultés qui semblent bien provenir de défauts inhérents à la méthode
elle-même. En particulier, la détermination des impédances acoustiques se fait de proche en proche en commençant par celles qui caractérisent les terrains les plus voisins de la surface. On conçoit alors facilement que les erreurs commises sur les premières valeurs se propagent et se répercutent sur toutes les autres valeurs. Or le début de l'enregistrement réel est souvent le moins bien connu, on a donc une probabilité d'erreur importante sur la détermination des impédances acoustiques des premières couches et, par voie de conséquence, sur la détermination de l'ensemble des valeurs.

Plutôt que d'utiliser une méthode de calcul des impédances acoustiques de proche en proche, il semble préférable d'employer un procédé itératif de détermination globale des impédances qui répartit les erreurs sur l'ensemble des valeurs à calculer. C'est une méthode de ce genre qui est développée dans cet article. On étudie d'abord comment se pose le problème. On donne ensuite des théorèmes de convergence et de continuité qui permettent d'expliquer certaines difficultés rencontrées dans la recherche de la solution. On propose une résolution numérique du problème et l'on montre des résultats obtenus par traitement numérique.
# 1 - POSITION DU PROBLEME

## 1.1 - Notations

On adopte les notations habituelles suivantes.

$z$ : désigne la profondeur ; $z$ est positif ou nul. (1.1)

$t$ : désigne le temps, $t$ est positif ou nul. (1.2)

$y(z,t)$ : désigne l'ébranlement à la profondeur $z$, à l'instant $t$. (1.3)

$\rho(z)$ : désigne la densité du milieu à la profondeur $z$. (1.4)

$\mu(z)$ : désigne le coefficient d'élasticité du milieu à la profondeur $z$. (1.5)

$\sigma(z)$ : $\sigma(z) = (\rho(z) \mu(z))^{1/2}$ désigne l'impédance acoustique du milieu à la profondeur $z$. (1.6)

$v(z)$ : $v(z) = (\mu(z)/\rho(z))^{1/2}$ désigne la célérité de propagation de l'ébranlement dans le milieu à la profondeur $z$. (1.7)

$x(z)$ : $x(z) = \int_0^z dz' \text{ représente le temps de parcours :}$ (1.8) c'est le temps mis par un ébranlement pour se propager de la profondeur $0$ à la profondeur $z$.

$\sigma(x)$ : désigne l'impédance acoustique du milieu en fonction du temps de parcours. (1.9)

## 1.2 - Equation retenue pour représenter le phénomène

Le modèle retenu ici est le plus simple possible. On suppose en effet que l'ébranlement est la solution de :

$$ \rho(z) \frac{\partial^2 y}{\partial t^2} - \frac{\partial}{\partial z}(\mu(z) \frac{\partial y}{\partial z}) = 0 \quad \text{pour } z>0 \text{ et } t>0 \qquad (1.10) $$

avec les conditions initiales

$$ \left\{ \begin{array}{l} y(z,0) = 0 \\ \frac{\partial y}{\partial t}(z,0) = 0. \end{array} \right. \quad \text{pour } z \ge 0 \qquad (1.11) $$

On suppose donc le milieu au repos à l'instant initial $t = 0$ et avec la condition aux limites

$$ - \mu(0) \frac{\partial y}{\partial z}(0,t) = g(t) \quad \text{pour } t > 0 \qquad (1.12) $$

où $g(t)$ représente l'impulsion en pression que l'on communique au milieu en surface à chaque instant $t$. Ce modèle est le plus simple possible, car c'est un modèle à une dimension d'espace et les termes d'amortissement sont négligés.
## 1.3 - Transformation de l'équation

Il est essentiel, pour le problème d'inversion de transformer l'équation (1.10) en changeant la variable z en la variable temps de parcours x, ar dans le plan (x,t) un ébranlement se propage avec une célérité égale à 1, indépendante de σ. Un ébranlement provoqué en x = 0 à l'instant t = 0, est réfléchi par un saut d'impédance au temps de parcours x' au temps t' = x' et revient donc en x = 0 au temps t' = 2 x'.

Le système (1.10), (1.11), (1.12) s'écrit alors :

$$ \sigma(x) \frac{\partial^2 y}{\partial t^2} - \frac{\partial}{\partial x} \left( \sigma(x) \frac{\partial y}{\partial x} \right) = 0 \quad x > 0 ; T > 0 \quad (1.13) $$

avec les conditions initiales

$$ \begin{cases} y(x,0) = y_0(x) = 0 \\ \frac{\partial y}{\partial t}(x,0) = y_1(x) = 0 \end{cases} \quad \text{pour } x \ge 0 \quad (1.14) $$

et avec la condition aux limites

$$ - \sigma(0) \frac{\partial y}{\partial x}(0,t) = g(t) \quad \text{pour } t > 0 \quad (1.15) $$

Il est utile d'introduire la formulation variationnelle du problème (1.13), (1.14), (1.15) : cf. Lions [8]. La fonction y vérifie pour toute fonction z régulière et pour t appartenant à [0,T]

$$ \frac{d^2}{dt^2} \int_0^x \sigma(x) y(x,t) z(x) \, dx + \int_0^x \sigma(x) \frac{\partial y}{\partial x}(x,t) \frac{\partial z}{\partial x}(x) \, dx = g(t) z(0) \quad (1.16) $$

avec $x_0 = \frac{T}{2}$, où T désigne la durée d'observation.

### Définition de l'application $H_g$

L'impulsion g(t) est supposée connue sur [0,T] (de même σ(0)). A toute fonction σ(x) définie sur [0,x], on associe la "réponse" du modèle correspondant à savoir la fonction $Y_σ(t)$ :

$$ Y_{\sigma}(t) = y(0,t) \quad \text{où} $$

$y(0,t)$ est la solution de (1.13), (1.14), (1.15) pour $σ(x)$.
Dans le cas général, (les conditions initiales $y_0(x)$ et $y_1(x)$ définies en (1.14) étant différentes de 0) la réponse $Y_\sigma$ sur l'intervalle $[0,t]$ ne dépend que de $\sigma(x)$ sur le même intervalle $[0,t]$. Dans le cas envisagé ici ($y_0(x) = y_1(x) = 0$), la réponse $Y_\sigma$ sur l'intervalle $[0,t]$ ne dépend que de $\sigma(x)$ sur l'intervalle $[0, \frac{t}{2}]$.

On note $H_g$ l'application qui à $\sigma$ fait correspondre la réponse $Y_\sigma$ :

$$ \sigma \xrightarrow{H_g} Y_\sigma \quad (\sigma \text{ et } Y_\sigma \text{ sont des fonctions définies sur } [0,t]) \quad (1.17) $$

L'indice g de $H_g$ vise à rappeler que cette application dépend de g : elle en dépend en fait linéairement.

$H_g(\sigma) + H_g(\sigma) = H_g(\sigma)$ pour toute fonction $\sigma$, cette propriété est bien sûr liée à la linéarité de l'équation (1.13), (1.14), (1.15).

## 1.4 - Position du problème : ses difficultés

### 1.4.1 - Le problème inverse

Pour un g donné, il s'agit d'inverser l'application $H_g$. On suppose donc connu $g(t)$ l'impulsion ; on observe la réponse $Y_\sigma(t)$ sur une période $[0,T]$ et le problème est de déterminer :

$$ \sigma = H_g^{-1}(Y_\sigma) \quad (1.18) $$

En écrivant (1.18), on fait implicitement l'hypothèse que le problème inverse est bien posé, à savoir :

$$ Y_{\sigma_1} = Y_{\sigma_2} \Longleftrightarrow \sigma_1 = \sigma_2 \quad (1.19) $$

En réalité, l'observation n'est pas la réponse $Y_\sigma(t)$ mais un sismogramme bruité $Y_d(t)$ et le problème est alors :

Quelles informations peut-on tirer sur $\sigma$ de la connaissance de $g$ et $Y_d(t)$ ?

Pour résoudre numériquement le problème d'inversion, on l'a posé comme un problème de minimisation.
### 1.4.2 - Formulation du problème d'inversion comme un problème de minimisation

Trouver une répartition $\sigma_{opt}$ dans la classe $\Sigma$ (voir 1.4.3) qui minimise

$$J(\sigma) = \int_{0}^{T} |Y_{\sigma}(t) - Y_{d}(t)|^{2} dt \qquad (1.20)$$

c'est-à-dire l'écart quadratique entre le sismogramme synthétique $Y_\sigma(t)$ défini en (1.17) et le sismogramme observé $Y_d(t)$.

$T$ représente la durée d'enregistrement du sismogramme

### 1.4.3 - Les difficultés du problème

#### i) Les difficultés de nature physique

Comme on l'a déjà indiqué, l'observation est bruitée : quelle est la nature du bruit ?

On ne connaît pas en fait avec exactitude l'impulsion $g(t)$ : peut-on espérer identifier $g$ et $\sigma$ à la fois ?

Enfin, le modèle n'est qu'une approximation de la réalité. On a par exemple négligé dans l'équation tous les termes d'amortissement.

#### ii) Difficulté mathématique

La difficulté essentielle est de déterminer une classe $\Sigma$ de fonctions $\sigma$ pour laquelle le problème (1.4.2) admet une solution unique stable par rapport aux données. Cette bonne classe de fonctions $\Sigma$ doit bien sûr correspondre aussi à des critères physiques.

#### iii) Difficulté informatique

La difficulté est celle de la taille du problème. Pour traiter un exemple réel, on est amené à introduire une discrétisation (300 x 1000) et donc à calculer l'ébranlement en $3 \times 10^5$ points. Les méthodes numériques d'optimisation utilisées exigent au moins 200 itérations, c'est-à-dire en fait 1000 calculs environ de l'état ou de l'état adjoint soit globalement environ $3 \times 10^8$ valeurs à calculer et donc $3 \times 10^9$ opérations élémentaires.
## 2 - QUELQUES PROPRIETES DE L'ÉQUATION DES ONDES ET LEURS CONSEQUENCES POUR LE PROBLEME INVERSE

On présente dans ce paragraphe quelques résultats théoriques de nature à faciliter la compréhension du problème et de ses difficultés. Ces résultats sont démontrés dans [3].

### 2.1 - Un résultat de convergence

#### 2.1.1 - Enoncé du résultat

**Théorème 1**

Soit une suite de fonctions $\sigma_N(x)$ vérifiant

$$0 < C \le \sigma_N(x) \quad \forall x \in [0, x_0] \quad \forall N \quad (2.1)$$

$$\lim_{N \to +\infty} \int_{x_1}^{x_2} \sigma_N(x) dx = \int_{x_1}^{x_2} \sigma^*(x) dx \quad \forall x_1, x_2; 0 \le x_1 \le x_2 \le x_0 \quad (2.2)$$

$$\lim_{N \to +\infty} \int_{x_1}^{x_2} \frac{dx}{\sigma_N(x)} dx = \int_{x_1}^{x_2} \frac{dx}{\bar{\sigma}(x)} \quad \forall x_1, x_2; 0 \le x_1 \le x_2 \le x_0 \quad (2.3)$$

Alors la suite $y_{\sigma_N}(x,T)$ converge uniformément sur $[0,x_0] \times [0,T]$ vers $y(x,T)$ solution de :

$$\left\{
\begin{array}{l}
\begin{aligned}
& \bullet \sigma^*(x) \frac{\partial^2 y}{\partial t^2} - \frac{\partial}{\partial x}(\bar{\sigma}(x) \frac{\partial y}{\partial x}) = 0 \\
\\
& \bullet - \bar{\sigma}(0) \frac{\partial y}{\partial x}(0,t) = g(t)
\end{aligned} \\
\\
\bullet \text{CI} \left[ \begin{array}{l} y(x,0) = 0 \\ \frac{\partial y}{\partial t}(x,0) = 0 \end{array} \right.
\end{array}
\right.
\quad (2.4)$$

en particulier $Y_{\sigma_N}(t)$ converge uniformément sur $[0,T]$ vers $Y(t) = y(0,t)$.
On note que l'équation (2.4) est du même type que (1.10) avec :

$$
\left\{
\begin{array}{ll}
. \text{ une densité : } & \rho(x) = \sigma^*(x) \\
. \text{ un coefficient d'élasticité : } & \mu(x) = \bar{\sigma}(x) \\
. \text{ une vitesse de propagation : } & C_a(x) = \left( \frac{\bar{\sigma}(x)}{\sigma^*(x)} \right)^{\frac{1}{2}}
\end{array}
\right.
\qquad (2.5)
$$

L'équation (2.4) reformulée en variable temps de parcours x', avec

$$
\frac{dx'}{dx} = \frac{1}{C_a(x)} \qquad (2.6)
$$

s'écrit :

$$
\left\{
\begin{array}{l}
\displaystyle \sqrt{\sigma^* \bar{\sigma}} \frac{\partial^2 y}{\partial t^2} - \frac{\partial}{\partial x'} (\sqrt{\sigma^* \bar{\sigma}} \frac{\partial y}{\partial x'}) = 0 \\[1em]
\text{CI} \left\{ \begin{array}{l} y(x',0) = 0 \\ \frac{\partial y}{\partial t}(x',0) = 0 \end{array} \right.
\end{array}
\right.
\qquad (2.7)
$$

Il est essentiel de remarquer que l'on a :

$$
C_a(x) \le 1 \quad \forall x \in [0, x_0] \qquad (2.8)
$$

ceci résulte de l'inégalité

$$
\bar{\sigma}(x) \le \sigma^*(x) \quad \forall x \in [0, x_0] \qquad (2.9)
$$

Pour que l'on ait $C_a = 1$ (c'est-à-dire égalité dans (2.9)), il faut supposer
que $\sigma_N$ converge dans un sens plus fort que (2.1), (2.2), (2.3), par exemple

$$
\sigma_N(x) \to \sigma^*(x) \quad \forall x \in [0, x_0] \text{ et aussi } 0 < C \le \sigma_N(x) \le M \quad \left( \begin{array}{l} \forall x \in [0, x_0] \\ \forall N \end{array} \right)
$$

Cette vitesse $C_a(x)$ constitue pour une répartition $\sigma_N$, du moins
lorsque N est assez grand, la vitesse apparente de propagation.
2.1.2 - Exemples de suites $\sigma_N$ vérifiant (2.1), (2.2), (2.3)

i) Soit une répartition $\sigma(x)$ vérifiant (2.1) sur $[0,1]$. On définit
$$ \sigma_N(x) = \sigma(N \cdot (x - \frac{k}{N})) \text{ pour } \frac{k}{N} \le x < \frac{k+1}{N} $$

$\sigma_N$ est une fonction périodique, de période $\frac{1}{N}$.

<figure><img src="image_5.png" /><figcaption>Figure 1</figcaption></figure>

Il est facile de montrer que l'on a :

$$ \lim_{N \to +\infty} \int_{x_1}^{x_2} \sigma_N(x) \, dx = (x_2 - x_1) \left( \int_{0}^{1} \sigma(x) \, dx \right) \neq x_1, x_2 \quad 0 \le x_1 \le x_2 \le 1 $$

$$ \lim_{N \to +\infty} \int_{x_1}^{x_2} \frac{dx}{\sigma_N(x)} \, dx = (x_2 - x_1) \left( \int_{0}^{1} \frac{dx}{\sigma(x)} \right) \neq x_1, x_2 \quad 0 \le x_1 \le x_2 \le 1 $$

et donc :

$$ \sigma^*(x) = \sigma^* = \int_0^1 \sigma(x) \, dx $$

$$ \bar{\sigma}(x) = \bar{\sigma} = \left( \int_{0}^{1} \frac{dx}{\sigma(x)} \right)^{-1} $$

En général, $\left(\int_{0}^{1} \frac{dx}{\sigma(x)}\right)^{-1} \le \int_{0}^{1} \sigma(x)dx$, l'égalité n'ayant lieu que si $\sigma$ est une fonction constante.
ii) On illustre en (2.1.3) les résultats du théorème 1, par une simulation numérique sur un exemple du type précédent, avec $\sigma(x)$ comme indiquée sur la figure 2.

<figure><img src="image_6.png" /><figcaption>Figure 2</figcaption></figure>

Dans ce cas, la vitesse apparente s'exprime en fonction du rapport :

$$ 
\left\{ 
\begin{array}{l} 
q = \frac{\sigma_2}{\sigma_1} \quad \text{, par} \\ 
C_a(q) = \frac{2\sqrt{q}}{1+q} 
\end{array} 
\right. 
$$

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

Vitesse apparente
en fonction de q

<table><caption>Figure 3</caption><thead><tr><td>q</td><td>C<sub>a</sub></td></tr></thead><tbody><tr><td>0.1</td><td>0.58</td></tr><tr><td>0.2</td><td>0.76</td></tr><tr><td>0.3</td><td>0.84</td></tr><tr><td>0.4</td><td>0.90</td></tr><tr><td>0.5</td><td>0.94</td></tr><tr><td>0.6</td><td>0.96</td></tr><tr><td>0.7</td><td>0.98</td></tr><tr><td>0.8</td><td>0.98</td></tr><tr><td>0.9</td><td>1.00</td></tr><tr><td>1.0</td><td>1.00</td></tr></tbody></table>

$$ \text{On a } C_a(\frac{1}{q}) = C_a(q) \quad \forall q > 0 $$

On constate que la vitesse $C_a$ reste voisine de 1 pour une très large bande de q : ainsi on a pour :

$$ 0,4 < q \le 2,5 $$

$$ 0,9 \le C_a(q) \le 1,0 $$
### 2.1.3 - Simulation numérique d'un exemple

La simulation porte sur l'exemple décrit au ii) du 2.1.2. Elle permet d'illustrer les résultats du théorème 1, à savoir :

#### 1 - La convergence de la suite $Y_{\sigma_N}(t)$

Sur la figure 4, on constate bien la convergence des réponses $Y_{\sigma_N}(t)$ vers une réponse limite $Y(t)$ lorsque la fréquence augmente.

#### 2 - La vitesse apparente $C_a$

On la met en évidence d'une première façon (figure 5) en comparant les réponses en des points d'espace différents, ce qui permet de constater un retard de l'arrivée de l'excitation, retard proportionnel à la profondeur du point considéré.

On la met en évidence d'une seconde façon (figures 6 et 7) en observant la réponse au bord, pour la répartition précédente modifiée par l'introduction d'un coefficient de réflexion égal à + 1 au point d'abscisse $x_0$. On constate alors un retard à l'arrivée de la réflexion totale sur le mur en $x_0$.

### 2.1.4 - Conséquences sur le problème d'identification

L'exemple de la figure 6 montre que la réponse sur $[0,80]$ ne dépend (pratiquement) pas des valeurs de $\sigma$ supérieures à 50. En observant $Y(t)$ sur $[0,80]$ il n'est donc pas possible d'identifier les valeurs de $\sigma$ supérieures à 50.

On peut donc dire que de fortes oscillations jouent le rôle de "cache".

On conçoit bien maintenant l'importance de la **classe des fonctions** dans laquelle on cherche une solution au problème d'identification. Il est naturel d'introduire la classe $\Sigma$ ($c_1, c_2$, var).

$$ \Sigma(c_1, c_2, \text{var}) = \{\sigma|0 < c_1 \le \sigma(x) \le c_2 \quad \forall x \in [0, x_0] ; \text{var}(\sigma) \le \text{var}\} \quad (2.10) $$

où $\text{var}(\sigma)$ désigne la variation de $\sigma$ sur $[0, x_0]$ soit :

$$ \text{var}(\sigma) = \sup_{(x_i)} \left\{ \sum_{i=0}^{N-1} |\sigma(x_{i+1}) - \sigma(x_i)| \right\} \quad \text{pour } (x_i) ; 0 \le x_1 \le \dots \le x_i \le x_{N+1} \quad (2.11) $$
<figure><img src="image_8.png" /><figcaption>ILLUSTRATION DE LA CONVERGENCE</figcaption></figure>

Figure 4
RÉPARTITION $\sigma = (\sigma_1, \sigma_2, \sigma_1, \sigma_2, \dots)$ $r = \frac{\sigma_1 - \sigma_2}{\sigma_1 + \sigma_2} = +0.6$

IMPULSION $G(t)$

VITESSE APPARENTE $C_a = 0.8$

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

ORIGINE DES TEMPS : différente pour chacune des 3 réponses. Elle est prise égale à $T_o = X_o$
<figure><img src="image_10.png" /><figcaption>Figure 6</figcaption></figure>
$$ \text{soit } n = \frac{\sigma_1 - \sigma_2}{\sigma_1 + \sigma_2} = 0.6 $$

$$ \left\{ \begin{array}{l} \text{•de 0 à 100; } \sigma_1, \sigma_2, \sigma_1, \sigma_2, \dots \\ \sigma_1 = 1.0 \\ \sigma_2 = 0.25 \end{array} \right. $$

•pour $X \ge 100$ soit un coefficient de réflexion égal à +1 en $X = 100$

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

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

RÉPONSE EN $X = 0$

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

VITESSE APPARENTE sur $[0,100]$
$C_d = 0.8$

L'échelle en temps est double de l'échelle en espace.

EXCITATION G(t) Voir 4.1.2
figure 9.

Figure 7

50

100

150

200
Le problème d'identification ainsi formulé est un problème d'optimisation avec contrainte.

## 2.2 - Un résultat de continuité

Si on se restreint à la classe $\Sigma$ définie par (2.10), on a le résultat de continuité uniforme suivant. Ce résultat est déjà énoncé dans [5] mais avec d'autres hypothèses sur $g$ et est démontré dans [3] par des techniques différentes.

**Théorème 2**

$$ \forall \delta > 0 ; \exists \varepsilon (\delta, \Sigma) \text{ tel que } \nabla_{\sigma_1, \sigma_2} \in \Sigma \text{ vérifiant } \left| \int_0^\Sigma \sigma_2(x')dx' - \int_0^\Sigma \sigma_1(x')dx' \right| \le \varepsilon \quad \nabla x \in [0, x_0] \quad (2.12) $$

$$ \text{on a : } \quad \int_0^T |Y_{\sigma_1}(t) - Y_{\sigma_2}(t)|^2 dt \le \delta $$

Ce théorème démontre la continuité de $H_g$ sur $\Sigma$ pour la topologie définie par (2.12). Sous réserve d'un résultat d'unicité, il permet de montrer aussi la **stabilité** du problème inverse, lorsqu'on se restreint à cette classe $\Sigma$.

## 2.3 - Un résultat de comportement lorsque $t$ tend vers l'infini

### 2.3.1 - Enoncé du résultat

On suppose que $\sigma(x)$ a un comportement pseudo périodique à l'infini et que de façon précise on a :

$$ \left\{ \begin{array}{ll} \lim_{X_1 \to +\infty} & \displaystyle \int_{X_2-X_1}^{X_2} \sigma(x) dx = \sigma^* \\ \lim_{X_2 \to +\infty} & \displaystyle \int_{X_1}^{X_2} \frac{dx}{\sigma(x)} dx = \frac{1}{\sigma} \\ \lim_{X_1 \to +\infty} & \displaystyle \int_{X_2-X_1}^{X_2} \sigma(x) dx = \sigma^* \\ \lim_{X_2 \to +\infty} & \displaystyle \int_{X_1}^{X_2} \frac{dx}{\sigma(x)} dx = \frac{1}{\sigma} \end{array} \right. \quad (2.13) $$
Alors, la solution y(x,t) de

$$
\left\{
\begin{array}{l}
\sigma(x) \frac{\partial^2 y}{\partial t^2} - \frac{\partial}{\partial x} \left(\sigma(x) \frac{\partial y}{\partial x}\right) = 0 \\
-\sigma(0) \frac{\partial y}{\partial x}(0,t) = g(t) \\
y(x,0) = y_0(x) \\
\frac{\partial y}{\partial t}(x,0) = y_1(x)
\end{array}
\right.
\qquad (2.14)
$$

converge lorsque t tend vers l'infini et l'on a

$$
\lim_{t \to +\infty} y(x,t) = y^* = \left( \int_0^{+\infty} g(t)dt + \int_0^{+\infty} y_1(x)dx \right) \frac{1}{\sqrt{\bar{\sigma}}} \sigma^* ; \forall x \ge 0 \quad (2.15)
$$

Dans le problème envisagé, on a :

$$
y^* = 0 ; \quad \text{en effet :} \tag{2.16}
$$

$$
i) \quad \int_{0}^{+\infty} g(t)dt = 0
$$

l'impulsion (en pression) a une valeur moyenne nulle, c'est-à-dire que l'ébran-
lement revient à zéro.

$$
ii) \int_{0}^{+\infty} y_{1}(x)dx = 0
$$

le milieu est initialement au repos, et donc :

$$
\lim_{t \to +\infty} y(x,t) = 0 \quad \forall x \ge 0, \quad \forall \sigma^* \text{ et } \bar{\sigma} \text{ définis par (2.13).}
$$

2.3.2 - Illustration numérique : conséquence sur le problème d'identification

Sur les exemples traités, dont les résultats sont représentés sur
la figure 8, on constate une convergence rapide* des réponses Y(t) vers leur
valeur asymptotique donnée par (2.15).

<footnote>* *Remarque* : Ici se pose un problème d'échelle. On a le résultat important suivant. Pour une excitation g(t) égale à un Dirac si Y(t) est la réponse à σ(x), Y(λt) est la réponse à σ(λx). Une contraction en espace se traduit donc par une contraction en temps de même rapport.</footnote>
Par conséquent, pour une période d'observation "assez longue", la valeur asymptotique devenant prépondérante, la fonctionnelle définie en (1.20) dépend essentiellement du produit des valeurs moyennes aux grandes profondeurs de $\sigma$ et de $\frac{1}{\sigma}$ (cf. 2.15), et de la valeur moyenne de l'excitation $g(t)$. Dans le cas où cette dernière est nulle, ces valeurs moyennes de $\sigma$ n'interviennent pas et seront donc difficiles à identifier.
<figure><img src="image_14.png" /></figure>

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

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

1- RÉPARTITION 1
$$r = \frac{\sigma_1 - \sigma_2}{\sigma_1 + \sigma_2} = +0.4$$

2- RÉPARTITION 2
$$r = \frac{\sigma_1 - \sigma_2}{\sigma_1 + \sigma_2} = +0.6$$

<caption>Figure 8</caption>
# 3 - RESOLUTION NUMÉRIQUE DU PROBLÈME

On donnera successivement des indications sur la façon dont a été résolue l'équation des ondes, puis le problème de contrôle associé.

## 3.1 - Résolution de l'équation d'état

On trouvera ci-dessous l'approximation du modèle en temps de parcours (1.13), (1.14), (1.15). L'approximation du modèle (1.10), (1.11), (1.12) en variable de profondeur est similaire.

### 3.1.1 - Discrétisation en espace

On utilise une méthode d'éléments finis de Lagrange de degré 3. L'intervalle $[0, x_0]$ est découpé en N intervalles égaux de longueur $h = \frac{x_0}{N}$ : on pose $x_{3i} = ih$ pour $i = 0, 1, ... N$. On va chercher, à chaque instant $t$, à approcher la fonction $y$ représentant l'état vibratoire du milieu à cet instant par une fonction $y_h$ appartenant à l'espace $V_h$ des fonctions continues sur $[0, x_0]$ et dont la restriction à chaque intervalle $[x_{3i}, x_{3i+3}]$ est un polynôme de degré 3.

$$ \text{Posons : } x_i = i \frac{h}{3} \quad i = 0, 1, \dots, 3N \qquad (3.1) $$
$$ \begin{array}{c} \begin{array}{ccccccc} x_0 & x_1 & x_2 & x_3 & & & \\ & \ddots & \ddots & \ddots & X_{3i} & X_{3i+3} & X_{3N} \\ \hline 0 & & & & h & & \\ \cline{2-7} & & & & & x_0 \end{array} \\ \rightarrow 0 \end{array} $$

Il est alors facile de construire des fonctions $w_k$ appartenant à cet espace $V_h$ et telles que :

$$ w_k(x_i) = \begin{cases} 0 & \text{si } k \neq i \\ 1 & \text{si } k = i \end{cases} \quad \forall k, i = 0, 1, \dots, 3N \quad (3.2) $$

Les fonctions $w_k$ forment une base de l'espace $V_h$ et par conséquent la détermination de l'état vibratoire approché à l'instant $t$, $y(x,t)$, se réduit
à celle de ses 3 N + 1 composantes $y_{ih}(t)$ $i = 0, 1, ..., 3N$ sur la base précédente :

$$ y_h(x,t) = \sum_{i=0}^{3N} y_{ih}(t) w_i(x) \quad (3.3) $$

[et par construction de la base $w_k$ - cf. (3.2) - on a $y_{ih}(t) = y_h(x_i, t)$].

Ces composantes $y_{ih}(t)$ sont déterminées en écrivant que la fonction $y_h(x,t)$ doit satisfaire la formulation variationnelle (1.16) pour toute fonction test appartenant à l'espace $V_h$. D'où en prenant successivement toutes les fonctions de base $w_k$ comme fonctions test

$$ \cdot \sum_{i=0}^{3N-1} \frac{d^2 y_{ih}(t)}{dt^2} \int_0^{x_0} \sigma(x) w_i(x) w_k(x) dx + \sum_{i=0}^{3N-1} y_{ih}(t) \int_0^{x_0} \sigma(x) \frac{dw_i}{dx} \frac{dw_k}{dx} dx = \\ g(t) w_k(0) \quad \forall k = 0, 1, ..., 3N-1 \quad (3.4) $$

$$ \cdot y_{ih}(0) = \frac{dy_{ih}}{dt}(0) = 0 \quad i = 0, 1, ..., 3N - 1 \quad (3.5) $$

$$ \cdot y_{3N,k}(t) = 0 \quad (\text{condition aux limites en } x = x_0) \quad (3.6) $$

L'impédance $\sigma(x)$ est supposée continue et linéaire par morceaux sur des intervalles de longueur $N_0h$, où $N_0$ est un entier à notre disposition(*) - plus $N_0$ est grand, plus la représentation de $\sigma(x)$ est grossière -

$$ \left\{ \begin{array}{ll} \forall n = 0, 1, ..., N/N_0; & \forall x \in [x_{3N_0n}, x_{3N_0(n+1)}] \\ \text{on a :} & \sigma(x) = \sigma_n + \frac{x - x_{3N_0n}}{N_0h} (\sigma_{n+1} - \sigma_n) \quad (3.7) \\ \text{où } \sigma_n = \sigma(x_{3N_0n}) & \end{array} \right. $$

On peut alors calculer les intégrales apparaissant dans (3.4) (par une formule d'intégration numérique pour la lère, exactement pour la 2ème). On obtient ainsi pour le vecteur $y_h(t) \in \mathbb{R}^{3N}$ de composantes $y_{ih}(t)$ $i = 0, 1, ..., 3N-1$, le système différentiel suivant :

<footnote>(*) On suppose que N est un multiple de $N_0$.</footnote>
$$
\begin{cases}
h A \frac{d^2 y_h(t)}{dt^2} + \frac{1}{h} B y_h(t) = g(t) e_1 \quad \forall t \in [0, T] \\
y_h(0) = \frac{dy_h(0)}{dt} = 0
\end{cases}
\quad (3.8)
$$

où $e_1$ est le 1er vecteur de la base canonique de $\mathbb{R}^{3N}$ et où A et B sont deux
matrices symétriques définies positives 3 N x 3 N respectivement diagonales
et heptadiagonales, indépendantes de h, mais dont les éléments non nuls sont
des fonctions linéaires de l'impédance $\sigma$.

3.1.2 - Discrétisation en temps

On divise l'intervalle $[0,T]$ en $N_t$ intervalles de longueur $k=\frac{T}{N_t}$,
et on approche la solution $y_h(t)$ du système différentiel (3.8) par une
suite de vecteurs $y_h^j \in \mathbb{R}^{3N}$, $j = 0, 1, \dots, N_t$ définie par

$$
\begin{cases}
A \frac{y_h^{j+1} - 2 y_h^j + y_h^{j-1}}{k^2} + \frac{1}{h^2} B y_h^j = \frac{1}{h} g^j e_1, & j = 1, 2, \dots, N_t - 1 \\
y_h^0 = y_h^1 = 0
\end{cases}
\quad (3.9)
$$

où $g^j$ est une approximation de $g(jk)$.

Ce schéma explicite en temps permet un calcul facile de $y_h^{j+1}$ à partir
de $y_h^{j-1}$ et $y_h^j$ puisque la matrice A est diagonale.

3.1.3 - Conditions de stabilité

La matrice A étant diagonale et B heptadiagonale, le schéma numérique
(3.9) fait intervenir au plus, dans le calcul de $y_{ih}^{j+1}$, les valeurs $y_{i-3,h}^j$ à
La célérité de l'ébranlement étant égale à l'unité, il en résulte une condition
<u>nécessaire de stabilité</u>

$$
k / h \le 1 \tag{3.10}
$$

En utilisant le schéma numérique (3.9), on peut démontrer une
<u>condition suffisante</u> asymptotique de <u>stabilité</u>.

$$
k / h < 1/6,5 \tag{3.11}
$$

Le rapport k/h utilisé lors des expériences numériques (cf. § 4) était
intermédiaire entre (3.10) et (3.11) et plus proche de (3.11) !
### 3.1.4 - Commentaires sur le choix de la discrétisation

Des essais comparatifs (cf. [3]) portant sur des approximations du modèle (1.10), (1.11), (1.12).

* par différences finies avec schéma explicite en temps.

* par éléments finis de Lagrange, avec schéma explicite en temps et intégration numérique (décrite ci-dessus).

* par éléments finis de Hermite cubique, avec un schéma explicite ou implicite en temps.

sont résumés dans le tableau suivant :

<table><caption>Tableau 1</caption><thead><tr><th>Approximation</th><th>Lagrange expli-cite + integ.num.</th><th>Hermite explicite</th><th>Hermite implicite</th><th>Différences finies</th></tr></thead><tbody><tr><td>Pas d'espace</td><td>20 m</td><td>20 m</td><td>20 m</td><td>5 m</td></tr><tr><td>Pas de temps</td><td>2 ms</td><td>2 ms</td><td>4 ms</td><td>1 ms</td></tr><tr><td>Temps de calcul<br>IBM 360.91</td><td>8 s</td><td>10 s</td><td>5,6 s</td><td>4 s</td></tr><tr><td></td><td>Bonne</td><td>Bonne</td><td>Mauvaise</td><td>Moyenne</td></tr><tr><td></td><td colspan="4">Impulsions obtenues après une réflexion partielle de l'impulsion de gauche</td></tr></tbody></table>

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

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

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

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

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

La nécessité d'une simulation précise pour pouvoir identifier correctement le paramètre $σ(x)$ a conduit à retenir le schéma de Lagrange bien que n'étant pas le plus rapide, car c'est le plus précis ; il présente en outre l'avantage de ne pas nécessiter de résolution de systèmes linéaires à chaque
pas de temps ce qui est avantageux lorsqu'on en vient à traiter des exemples réels de taille importante...

Par contre le choix d'approcher $\sigma(x)$ par une fonction continue et linéaire sur un certain nombre d'éléments finis n'est pas lié à l'approximation du problème par éléments finis de Lagrange ; tout autre choix serait loisible (par exemple, $\sigma(x)$ constant sur chaque élément fini).

### Remarque

A une profondeur $X$ donnée, l'ébranlement $y$ est nul pour $t \le x$ (phénomène de propagation). On a tenu compte de ce fait dans le programme numérique pour ne pas effectuer des calculs inutiles.

## 3.2 - Résolution du problème d'identification

Soit $y_d^j, j = 1, 2, ... N_t$ le sismogramme. Le critère (1.20) s'écrit :

$$J(\sigma) = k \sum_{j=1}^{N_t} (y_{oh}^j - y_d^j)^2 \quad (3.12)$$

où $y_h^j$ est calculé à partir de $\sigma$ à l'aide des équations d'état (3.9).

$\sigma$ est le vecteur $\sigma_0, \sigma_1, ..., \sigma_N^{(*)}$ (cf. (3.7))

On est donc ramené à minimiser une fonction $J$ de $N + 1$ variables.

### 3.2.1 - Méthode d'optimisation

Le critère $J$ pouvant dépendre d'une centaine de variables, toute méthode d'optimisation de type exploratoire est exclue. On a utilisé une méthode de gradient conjugué (Fletcher Reeves). On est ainsi assuré d'obtenir à chaque itération $m$ de gradient, une erreur $J^m$ sur le sismogramme inférieure à la valeur $J^{m-1}$ à l'itération précédente.

La répartition d'impédance $\sigma^m$ convergerait vers une valeur réalisant le minimum de $J$ en $N + 1$ itérations exactement si $J$ était une fonction quadra-tique de $\sigma$.

<footnote>(*) On suppose dans ce § $N_0 = 1$ pour simplifier l'exposé.</footnote>
Or le critère considéré ici n'est pas, par construction, une fonction quadratique de $\sigma_0, \sigma_1, \dots, \sigma_N$ : on peut alors seulement espérer la convergence, lorsque le nombre m d'itérations de gradient tend vers l'infini, de $\sigma^m$ vers une répartition d'impédance $\sigma^*$ correspondant seulement à un minimum local du critère J, ou même à un point stationnaire quelconque.

En pratique, une itération de gradient, c'est-à-dire le passage de $\sigma^m$ à $\sigma^{m+1}$, nécessite :

- le calcul du gradient $\nabla J$ de $J$ : cf. § 3.2.2

- une minimisation monodimensionnelle : elle est réalisée de façon approchée en balayant $\mathbb{R}^+$ par une progression géométrique de raison deux, ce qui nécessite quelques évaluations du critère J.

### 3.2.2 - Calcul du gradient de J

On utilise ici le calcul des variations (ou la théorie du contrôle) pour le calcul de ce gradient. C'est la méthode la plus précise et la plus économique : elle ne requiert que la résolution de l'équation d'état (3.9) et de l'équation adjointe.

$$ \left\{ \begin{array}{l} A \displaystyle\frac{p^{j-1} - 2p^j + p^{j+1}}{k^2} + \frac{1}{h^2} Bp^j = -2k(y_{oh}^j - y_{d}^j) e_1, \quad j = N_t, \dots, 2, 1 \\ p^{N_t+1} = p^{N_t} = 0 \end{array} \right. \qquad (3.13) $$

où $p^0, p^1, \dots, p^{N_t+1}$ sont des vecteurs de $\mathbb{R}^{3N}$ formant l'état adjoint.

On a alors en notant $\delta J, \delta \sigma, \delta A, \dots$ les variations au 1er ordre des quantités correspondantes :

$$ \delta J = \sum_{j=1}^{N_t} \left\langle \delta A \frac{y_h^{j+1} - 2y_h^j + y_h^j}{k^2} + \frac{1}{h^2} \delta B y_h^j, p_j \right\rangle \quad (3.14) $$

Les matrices $\delta A$ et $\delta B$ étant des fonctions linéaires de $\delta \sigma_0, \dots, \delta \sigma_N$ la formule (3.14) donne après une simple substitution le gradient de $J$ par rapport à $\sigma_1, \dots, \sigma_N$ : $\frac{\partial J}{\partial \sigma_i}$ est le coefficient de $\delta \sigma_i$ au 2ème membre de (3.14).
Démontrons (3.14). En différentiant (3.9) et (3.12), il vient :

$$
A \frac{\delta y_h^{j+1} - 2 \delta y_h^j + \delta y_h^{j-1}}{k^2} + \frac{1}{h^2} B \delta y_h^j = - A \delta \frac{y_h^{j+1} - 2 y_h^j + y_h^{j-1}}{k^2} - \frac{1}{h^2} \delta B y_h^j \\
\text{pour } j = 1, 2, \dots, N_t - 1
$$

(3. 15)

$$
\delta y_h^0 = \delta y_h^1 = 0
$$

(3. 16)

$$
\delta J = 2 k \sum_{j=1}^{N_t} (y_{oh}^j - y_{d}^j) \delta y_{oh}^j = < 2 k \sum_{j=1}^{N_t} (y_{oh}^j - y_{d}^j) e_j, \delta y_h^j > \quad (3.17)
$$

Faisons le produit scalaire de (3.15) par $p^j$ dans $\mathbb{R}^{3N}$ et sommons
de $j=1$ à $N_t - 1$. En remarquant que :

$$
\sum_{j=1}^{N_t-1} \langle A \delta y_h^{j+1}, p^j \rangle = \sum_{j=0}^{N_t-1} \langle A \delta y_h^{j+1}, p^j \rangle \text{ puisque } \delta y_h^1 = 0
$$

$$
\sum_{j=1}^{N_t-1} \langle A \delta y_h^j, p^j \rangle = \sum_{j=1}^{N_t} \langle A \delta y_h^j, p^j \rangle \quad \text{puisque } p^{N_t} = 0
$$

$$
\sum_{j=1}^{N_t-1} \langle A \delta y_h^{j-1}, p^j \rangle = \sum_{j=2}^{N_t+1} \langle A \delta y_h^{j-1}, p^j \rangle \quad \text{puisque } \delta y_h^0 = p^{N_t} = p^{N_t+1} = 0
$$

et en translatant les indices, il vient puisque A et B sont des matrices
symétriques

$$
\sum_{j=1}^{N_t} \langle A \frac{p^{j-1} - 2p^j + p^{j+1}}{k^2} + \frac{1}{h^2} B p^j, \delta y^j \rangle = \\
- \sum_{j=1}^{N_t} \delta A \frac{y_h^{j+1} - 2y_h^j + y_h^{j-1}}{k^2} + \frac{1}{h^2} \delta B y_h^j, p^j \rangle
$$

ce qui donne (3.14) compte tenu de (3.13).
### 3.2.3 - Mise en oeuvre

Le calcul du gradient J à partir de (3.14) nécessite la connaissance de tous les vecteurs $y^j$ et $p^j$ ; $j = 1, 2, ... N_t$. On évite le stockage de ces deux familles de vecteurs (soit deux tableaux de 300 x 1000 dans les exemples numériques) de la façon suivante :

- on résoud les équations d'état (3.9) et on ne stocke que les deux derniers vecteurs $y_h^{N_t}$ et $y_h^{N_t+1}$ (300 x 2).

- on résoud dans le sens rétrograde (l'équation des ondes étant réversible) les mêmes équations d'état (3.9) ainsi que les équations adjointes (3.13). On dispose ainsi, à chaque pas de temps j, des vecteurs $y_h^j$ et $p^j$ nécessaires au calcul du gradient.

Le temps supplémentaire requis pour la 2ème intégration du système direct est nettement inférieur à celui dû à la manipulation d'un tableau de très grande taille qui ne tient pas en mémoire centrale.

### 3.2.4 - Organigramme

On donne, pour simplifier l'organigramme correspondant à une méthode d'optimisation de type gradient pur. Le cas d'une méthode de gradient conjugué est similaire en un peu plus complexe.
ORGANI GRAMME

<figure><img src="image_22.png" /></figure>
# 4 - RESULTATS NUMERIQUES

On donne ici les résultats obtenus sur 3 exemples. Sur un premier exemple, de caractère académique, on vérifie la propriété fondamentale de l'algorithme qui était escomptée : la stabilité par rapport au bruit. Un second exemple, plus proche de la réalité, fait apparaître de nouvelles difficultés. Enfin, sur le dernier exemple, on aborde le problème de la stabilité par rapport à la discrétisation, en vue du traitement de sismogrammes non synthétiques.

## 4.1 - Choix des différents paramètres et description des exemples traités

### 4.1.1 - Modèle et algorithme retenus

#### 1 - Choix du modèle

Des essais ont montré que le problème est beaucoup mieux posé en variable temps de parcours qu'en variable profondeur ; la vitesse de convergence de l'algorithme est beaucoup plus grande. En plus, grâce à un choix optimal du rapport k/h, l'intégration numérique de l'équation des ondes, est quatre fois plus rapide : dans ce modèle, la vitesse de propagation étant égale à 1, la condition de stabilité est indépendante de la répartition d'impédances.

#### 2 - Algorithme d'optimisation retenu

La convergence est très lente avec l'algorithme du gradient simple. L'utilisation d'un algorithme du gradient conjugué (type Fletcher-Reeves) permet d'augmenter la vitesse de convergence dans un rapport de 1 à 20 (et de diminuer d'autant le temps calcul !).

Avec le modèle choisi, l'algorithme retenu est performant sur des exemples tests. On l'a appliqué à des exemples de difficulté croissante qu'on décrit ci-dessous.
### 4.1.2 - Description des exemples traités

#### 1 - L'impulsion

On utilisera la même pour tous les essais. Il s'agit de la dérivée du signal classique de Ricker d'une largeur de 56 ms. Ce signal est représenté en figure 9.

<figure><img src="image_23.png" /><figcaption>Figure 9</figcaption></figure>

De cette impulsion, on déduit une condition sur le pas de discrétisation en espace h : il doit être au plus égal à une valeur $h_0$ déterminée par le contenu en fréquence de l'impulsion (4.1). Pratiquement, on choisit h voisin de $h_0$. Le pas de discrétisation en temps k est alors choisi tel que le rapport k/h soit voisin de celui donné par la condition de stabilité (3.11).
On a pris :

$$h = 10 \text{ ms} \qquad (4.2)$$

$$k = 2 \text{ ms} \qquad (4.3)$$

## 2 - Les répartitions à identifier

Tous les essais ont porté sur l'identification des 3 répartitions suivantes. On les note $\sigma_I, \sigma_{II}$, et $\sigma_{III}$.

On a choisi des répartitions $\sigma_i, i = I, II, III$ continues et linéaires (cf. § 3.1.1) sur des intervalles de longueur égale $N_{\sigma_i}$ h qu'on appellera couches :

$$ \begin{align*} N_{\sigma_I} &= 3 \\ N_{\sigma_{II}} &= 2 \end{align*} \qquad (4.4) $$

$$N_{\sigma_{III}} = 1$$

Le tableau ci-dessous résume les caractéristiques de ces trois répartitions et donne, à titre de comparaison, celles d'une répartition réelle.

<table><thead><tr><th>Tableau 2</th><th>σ<sub>I</sub></th><th>σ<sub>II</sub></th><th>σ<sub>III</sub></th><th>Exemple réel</th></tr></thead><tbody><tr><td>Longueur totale en temps de parcours</td><td>0,51 s</td><td>1,12 s</td><td>1,12 s</td><td>2 à 3 s</td></tr><tr><td>Epaisseur des couches</td><td>30 ms</td><td>20 ms</td><td>10 ms</td><td>2 à 4 ms</td></tr><tr><td>Nombre des couches</td><td>17</td><td>56</td><td>112</td><td>500 à 1500</td></tr><tr><td>Temps de parcours d'une couche/durée de l'impulsion</td><td>0,21</td><td>0,14</td><td>0,07</td><td>0,02</td></tr><tr><td>Variation (cf. 2.11)</td><td>2110</td><td>15200</td><td>22200</td><td>50000</td></tr></tbody></table>
<figure><img src="image_24.png" /><figcaption>Figure 10
Répartition $\sigma_1$</figcaption></figure>

Temps de parcours en s.

<figure><img src="image_25.png" /><figcaption>Figure 11
Répartition $\sigma_2$</figcaption></figure>

Temps de parcours en s.

<figure><img src="image_26.png" /><figcaption>Figure 12
Répartition $\sigma_3$</figcaption></figure>
# Commentaires sur les répartitions

La répartition $\sigma_I$ est un exemple académique ; elle est présentée ici à titre d'exemple test.

Les répartitions $\sigma_{II}$ et $\sigma_{III}$ sont relativement proches de répartitions réelles.

La répartition $\sigma_{III}$ est la répartition $\sigma_{II}$ à laquelle on a appliqué une légère perturbation.

## 3 - Les sismogrammes observés

### a) Sismogrammes non bruités :

Ce sont les réponses $Y_i$ (i = I, II, III) définies en (1.17). Elles ont été calculées ici par le schéma d'intégration numérique défini au § 3.1.

<figure><img src="image_27.png" /><figcaption>Figure 13</figcaption></figure>
<figure><img src="image_28.png" /><figcaption>Figure 14</figcaption></figure>

<figure><img src="image_29.png" /><figcaption>Figure 15</figcaption></figure>
b) Sismogrammes bruités

Ils sont générés en ajoutant aux sismogrammes non bruités un bruit gaussien de moyenne nulle que l'on a filtré par convolution avec l'impulsion. On obtient ainsi un bruit proche de celui observé expérimentalement (I.F.P.). Pour mesurer le niveau de bruit, on introduit le rapport bruit/signal, défini comme le rapport de l'écart-type du bruit à l'écart-type du signal (les écarts-types sont calculés à partir de la fin de l'impulsion). On représente ci-dessous des sismogrammes $Y_{dI}^{\alpha}$ correspondant chacun à une répartition $\sigma_i(i = I ou II ou III)$ et à un rapport bruit/signal égal à $\alpha$.

<figure><img src="image_30.png" /><figcaption>Sismogramme $Y_{dI}^{0,7}$</figcaption></figure>

<figure><img src="image_31.png" /><figcaption>Sismogramme $Y_{dI}^{1,15}$</figcaption></figure>

<caption>Figure 16</caption>

<figure><img src="image_32.png" /><figcaption>Sismogramme $Y_{dII}^{0,6}$</figcaption></figure>

<caption>Figure 18</caption>
### 4.1.3 - Essais effectués

Le problème est de trouver une répartition $\sigma$ approchant "au mieux" la répartition exacte $\sigma_i$, $i = I, II$ ou $III$.

Les résultats sont fonction du nombre $N_\sigma$ d'éléments finis consécutifs sur lesquels on suppose la répartition $\sigma$ linéaire.

Le cas $N_\sigma = N_{\sigma_i}$, $i = I, II$ ou $III$, cas peu réaliste est traité au § 4.2. (4.6)

Le cas $N_\sigma \neq N_{\sigma_i}$, $i = I, II$ ou $III$ est traité au § 4.3. (4.7)

Le tableau 3 présente l'ensemble des essais effectués.
<table><thead><tr><th>Tableau 3</th><th>Répartition à identifier</th><th>Sismogramme observé</th><th>N<sub>σ</sub></th><th>N<sub>σi</sub></th><th>Initialisation</th></tr></thead><tbody><tr><td>Essai 1</td><td>σ<sub>I</sub></td><td>γ<sub>σI</sub> = γ<sup>0</sup><sub>dI</sub></td><td>3</td><td>3</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 2</td><td>σ<sub>I</sub></td><td>γ<sup>0,7</sup><sub>dI</sub></td><td>3</td><td>3</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 3</td><td>σ<sub>I</sub></td><td>γ<sup>1,15</sup><sub>dI</sub></td><td>3</td><td>3</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 4</td><td>σ<sub>II</sub></td><td>γ<sub>σII</sub> = γ<sup>0</sup><sub>dII</sub></td><td>2</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 4 bis</td><td>σ<sub>II</sub></td><td>γ<sub>σII</sub> = γ<sup>0</sup><sub>dII</sub></td><td>2</td><td>2</td><td>σ<sup>0</sup><sub>2</sub></td></tr><tr><td>Essai 5</td><td>σ<sub>II</sub></td><td>γ<sup>0,6</sup><sub>dII</sub></td><td>2</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 5 bis</td><td>σ<sub>II</sub></td><td>γ<sup>0,6</sup><sub>dII</sub></td><td>2</td><td>2</td><td>σ<sup>0</sup><sub>2</sub></td></tr><tr><td>Essai 6</td><td>σ<sub>II</sub></td><td>γ<sub>σII</sub> = γ<sup>0</sup><sub>dII</sub></td><td>3</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 7</td><td>σ<sub>II</sub></td><td>γ<sup>0,6</sup><sub>dII</sub></td><td>3</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 8</td><td>σ<sub>II</sub></td><td>γ<sub>σII</sub> = γ<sup>0</sup><sub>dII</sub></td><td>4</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 9</td><td>σ<sub>II</sub></td><td>γ<sup>0,6</sup><sub>dII</sub></td><td>4</td><td>2</td><td>σ<sup>0</sup><sub>I</sub></td></tr><tr><td>Essai 10</td><td>σ<sub>III</sub></td><td>γ<sub>σIII</sub> = γ<sup>0</sup><sub>dII</sub></td><td>2</td><td>1</td><td>σ<sup>0</sup><sub>2</sub></td></tr></tbody></table>

Dans le tableau 3, $σ_1^0$ et $σ_2^0$ désignent les deux répartitions initiales utilisées pour l'algorithm d'optimisation.
<figure><img src="image_33.png" /></figure>

Initialisation $\sigma_1^0$

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

Initialisation $\sigma_2^0$
## 4.2 - Essais à épaisseur de couches connue

On donne ici les résultats obtenus au cours des essais 1 à 5 qui correspondent au cas $N_σ = N_{σ_i}$.

### 4.2.1 - Exposé des résultats

- Essai 1 : il s'agit d'identifier la répartition $σ_I$ à partir du sismogramme non bruité $Y_{σ_I}$ correspondant, et en partant de la valeur initiale $σ_I^0$.

Les résultats obtenus sont très bons (figure 19). On constate d'autre part que le sismogramme généré par la répartition obtenue et le sismogramme $Y_{σ_I}$ sont superposables.

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

#### Essai 1

<caption>Figure 19</caption>
## Essais 2 et 3 :

Comparables à l'essai 1, mais à partir de sismogrammes comportant des rapports de 70 et 115 % de bruit.

Le bruit a altéré les résultats : il y a eu formation de pics (Fig. 20 et 21).

<figure><img src="image_36.png" /><figcaption>Essai 2</figcaption></figure>

<figure><img src="image_37.png" /><figcaption>Essai 3</figcaption></figure>

<caption>Figure 21</caption>
Cependant l'allure générale reste satisfaisante, malgré le niveau de bruit important. Ceci montre la bonne stabilité de l'algorithme par rapport au bruit.

Les sismogrammes générés par les répartitions obtenues (figures 22 et 23)sont intermédiaires entre les sismogrammes bruité $Y_{dI}^{\alpha}$ et non bruité $Y_{\sigma I}$.

<figure><img src="image_38.png" /><figcaption>Sismogramme généré par
la répartition obtenue
à l'itération 60</figcaption></figure>

Figure 22

<figure><img src="image_39.png" /><figcaption>Sismogramme généré par
la répartition obtenue
à l'itération 80</figcaption></figure>

Figure 23
- Essai 4 : on identifie la répartition (plus complexe) $σ_{II}$ à partir d'un sismogramme non bruité. La convergence est très lente. La figure 24 montre que la valeur moyenne n'est pas ajustée mais que la position des sauts et leur amplitude relative sont assez bien rendues.

- Essai 4 bis : comme essai 4, mais en partant d'une valeur initiale différente, on obtient de très bons résultats ; la convergence est plus rapide (figure 25).

Dans les essais 4 et 4 bis, les sismogrammes générés par les répartitions obtenues sont superposables avec le sismogramme $Y_{σ_{II}}$.

<figure><img src="image_40.png" /><figcaption>Figure 24</figcaption></figure>
<figure><img src="image_41.png" /><figcaption>Essai 4 bis</figcaption></figure>

Figure 25
- Essai 5 : comme essai 4, mais avec 60 % de bruit. Les résultats obtenus après 500 itérations de gradient conjugué sont assez bons (fig. 26 et 27). Ils appellent les mêmes commentaires qu'aux essais 2 et 3.

<figure><img src="image_42.png" /><figcaption>Essai 5

Figure 26</figcaption></figure>

<figure><img src="image_43.png" /><figcaption>Essai 5

Figure 27</figcaption></figure>
- Essai 5 bis : comme essai 4 bis, mais avec 60 % de bruit. La répartition obtenue après 90 itérations (fig. 28) approche la solution cherchée, mais elle est perturbée par des pics qui se développent lorsqu'on poursuit les itérations. Ainsi la variation croît rapidement : de 24000 à l'itération 50, elle passe à 33000 à l'itération 90.

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

Essai 5 bis

Figure 28
<figure><img src="image_45.png" /><figcaption>Figure 29</figcaption></figure>

<figure><img src="image_46.png" /><figcaption>Figure 30</figcaption></figure>
### 4.2.2 - *Commentaires sur le cheminement de l'algorithme*

#### a) *Cas non bruité* (essais 1, 4 et 4 bis)

La fonction coût (3.12) décroît et tend vers 0 (figure 31).

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

### Essai 1

#### Fonction coût au cours des itérations

<caption>Figure 31</caption>

On obtient plus facilement l'allure de la répartition que le positionnement autour de la valeur moyenne. C'est une constatation générale.

On présente sur la figure 32 l'évolution des écarts $L^1$ sur l'impédance et sur les coefficients de réflexion$^{(1)}$ au cours des itérations de gradient.

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

### Essai 1

#### Ecart $L^1$ sur l'impédance

<caption>Figure 32</caption>

(1) On définit l'écart $L^1$ sur l'impédance à l'itération n par :

$$ \frac{\displaystyle\sum_{j=1}^{n} | \sigma^n(x_j) - \sigma(x_j) |}{\displaystyle\sum_{j=1}^{n} \sigma^n(x_j)} $$

où $\sigma$ est la répartition exacte cherchée.
L'initialisation $\sigma_2^0$ donne de meilleurs résultats pour un même nombre d'itérations.

On remarque que l'impédance ne bouge pratiquement pas de sa valeur d'initialisation sur la dernière couche. En effet, on a pris, par sécurité, un terrain à identifier un peu plus profond que le sismogramme ne le permet ($x_0 > \frac{T}{2}$).

b) **Cas bruité** (essais 2, 3, 5 et 5 bis)

La fonction coût décroît mais ne tend plus vers 0 (figure 33).

<figure><img src="image_49.png" /><figcaption>Figure 33</figcaption></figure>

La convergence est moins lente que dans le cas non bruité.

Il y a formation de pics : ceci est dû au fait que l'algorithme n'impose pas de borne supérieure à la variation de $\sigma$ (définie en 2.11). D'après le théorème 2, on n'est donc pas assuré de l'existence d'une répartition optimale.

L'initialisation $\sigma_2^0$ donne une convergence plus rapide, mais semble augmenter les risques d'"explosion" de l'algorithme, causée par une croissance rapide de la variation.

Dans l'essai 2, il faut noter que l'algorithme se bloque à l'itération 60.
### 4.2.3 - Temps de calcul, vitesse de convergence

L'ordinateur utilisé est le CDC 7600 de la FRANLAB(&). On a fait figurer dans le tableau les paramètres importants.

<table><thead><tr><th>Tableau 4</th><th>Nombre d'itérations par minute</th><th>Nombre d'itérations effectuées</th><th>Qualité des résultats</th><th>Temps total Unité Centrale</th><th>Longueur du sismogramme en secondes</th><th>Niveau de bruit</th></tr></thead><tbody><tr><td>Essai 1</td><td>10 itér/min</td><td>80</td><td>Très bons perfectibles</td><td>9 min.</td><td>0,51 s</td><td>0 %</td></tr><tr><td>Essai 2</td><td>10 itér/min</td><td>60</td><td>Très bon algorithme bloqué</td><td>6 min.</td><td>0,51 s</td><td>70 %</td></tr><tr><td>Essai 3</td><td>10 itér/min</td><td>80</td><td>Très bons</td><td>9 min.</td><td>0,51 s</td><td>115 %</td></tr><tr><td>Essai 4</td><td>1,7 itér/min</td><td>400</td><td>Assez bons : perfectibles</td><td>4 h. 30</td><td>1,12 s</td><td>0 %</td></tr><tr><td>Essai 4 bis</td><td>1,7 itér/min</td><td>170</td><td>Très bons perfectibles</td><td>2 h.</td><td>1,12 s</td><td>0 %</td></tr><tr><td>Essai 5</td><td>1,7 itér/min</td><td>500</td><td>Assez bons</td><td>5 h. 30</td><td>1,12 s</td><td>60 %</td></tr><tr><td>Essai 5 bis</td><td>1,7 itér/min</td><td>90</td><td>Assez bons</td><td>1 h.</td><td>1,12 s</td><td>60 %</td></tr></tbody></table>

Le temps calcul est très important ; en revanche l'encombrement mémoire est insignifiant (moins de 10 000 mots, soit le cinquième de la mémoire centrale du 7600).

Il est important de remarquer que la durée d'une itération croît comme le carré de la longueur du sismogramme à traiter.

<footnote>(&) Le compilateur FORTRAN est l'optimiseur FTN de niveau 1 (version 4.3 + 393 mai 1975).</footnote>
## 4.3 - Essais à épaisseur de couches inconnues

Dans ce paragraphe, on aborde la question suivante : une répartition quelconque peut-elle être identifiée correctement à l'aide d'impédances continues et linéaires par morceaux de longueur donnée ? Peut-on se contenter d'une discrétisation grossière ? En bref, on pose le problème de la stabilité de l'identification par rapport à la discrétisation des coefficients à identifier.

Soit donc une des répartitions $\sigma_i$ ($i = II$ ou $III$). On cherche à l'identifier à l'aide de répartitions $\sigma$ avec $N_\sigma > N_{\sigma_i}$.

### 4.3.1 - Exposé des résultats

**Essais 6 et 7** : il s'agit d'identifier la répartition $\sigma_{II}$, linéaire sur 2 éléments ($N_\sigma = 2$) finis, par une répartition linéaire sur 3 éléments finis ($N_\sigma = 3$).

Les répartitions obtenues après 200 itérations de gradient conjugué sont semblables dans les deux cas, bruité et non bruité (figure 34) ; elles n'approchent pas trop mal la solution cherchée dans la moitié la moins profonde mais dans la deuxième partie il n'y a plus aucune ressemblance. A l'essai 7, l'écart $L^1$ est minimum à l'itération 140.

**Essais 8 et 9** : comme essais 6 et 7, mais avec $N_\sigma = 4$. L'algorithme se bloque à l'itération 35 (essai 8) et à l'itération 33 (essai 9). Les répartitions obtenues appellent les mêmes commentaires qu'aux essais 6 et 7.

La mauvaise qualité des résultats obtenus au cours des 4 derniers essais montre que le problème est très instable par rapport à la discrétisation. Il est donc nécessaire, pour améliorer les résultats, de considérer des discrétisations plus fines.

**Essai 10** : il s'agit d'approcher la répartition $\sigma_{II}$, linéaire sur 1 élément fini, par une répartition linéaire sur 2 éléments finis. Les résultats obtenus sont meilleurs. On peut l'expliquer ainsi :

1. La discrétisation est plus fine.

2. La répartition $\sigma_{III}$ est une perturbation de $\sigma_{II}$.

Or pour $\sigma_{II}$, on a : $N_{\sigma_{II}} = 2 = N_\sigma$.

Le sismogramme généré par la répartition trouvée (figure 38) approche assez bien le sismogramme $Y_{\sigma_{III}}$ (figure 37).
<figure><img src="image_50.png" /><figcaption>Figure 34</figcaption></figure>

Essais
6 et 7

$N_6 = 3$

<figure><img src="image_51.png" /><figcaption>Figure 35</figcaption></figure>

Essais
8 et 9

$N_6 = 4$
<figure><img src="image_52.png" /></figure>

Essai 10

Figure 36
<figure><img src="image_53.png" /><figcaption>Figure 37</figcaption></figure>

<figure><img src="image_54.png" /><figcaption>Figure 38</figcaption></figure>
### 4.3.2 - Conclusions sur la stabilité par rapport à la discrétisation des impédances

Le problème est très instable par rapport à la discrétisation des impédances. Ceci est dû à la nature même du problème inverse et non à l'algorithme utilisé pour le résoudre. Il faudra donc, pour traiter des sismogrammes non synthétiques, utiliser des discrétisations très fines, ce qui augmentera la dimension du problème et donc le temps calcul.

## 4.4 - Conclusions

La méthode exposée, de par son caractère global possède une bonne stabilité par rapport au bruit. Cette propriété n'a été acquise qu'au prix d'un temps calcul très important.

Cette étude a fait apparaître deux difficultés :

* un niveau de bruit élevé, qui peut faire "exploser" l'algorithme d'optimisation sans contrainte utilisé. On essaie actuellement de surmonter cette difficulté en introduisant une contrainte sur la variation de l'impédance, contrainte suggérée par la théorie.

* un temps de calcul prohibitif. Pour réduire ce temps calcul, on se propose d'utiliser la méthode des caractéristiques pour l'intégration numérique de l'équation des ondes.

La résolution de ces deux problèmes est un préalable au traitement de sismogrammes non synthétiques.
BIBLIOGRAPHIE

[1] - P. BOIS, G. GRAU, Ch. HEMON et M. LA PORTE : "Calcul automatique des sismogrammes synthétiques en ondes planes à incidence normale". Revue de l'I.F.P., vol. XVII, n° 4, avril 1962.

[2] - G. CHAVENT : "Deux résultats sur le problème inverse dans les équations aux dérivées partielles du deuxième ordre en t et sur l'unicité de la solution du problème inverse de la diffusion". C.R Acad. Sc. Paris, T 270, 25-28 janvier 1970.

[3] - G. CHAVENT, P. LAILLY, A. BAMBERGER : Rapport Laboria (à paraître).

[4] - D. M. DETCHMENDY : "Attenuation and inverse problems in wave propagation" (SWIEECO 68), Esso Production Research Company, Houston, Texas.

[5] - M.L. GERVER : "The inverse problem for the vibrating string equation". Izv Acad. Sci. URSS Physic Solid Earth 1970-1971.

[6] - G. KUNETZ : "Quelques exemples d'analyse d'enregistrements sismiques". Geophysical Prospecting, 1963, 11, n° 4.

[7] - G. KUNETZ, BARANOV : "Film synthétique avec réflexions multiples. Théorie et calcul pratiques". Geophysical Prospecting, 1960, 8, n° 3.

[8] - J.L. LIONS : "Contrôle optimal de systèmes gouvernés par des équations aux dérivées partielles. Dunod.

[9] - N. G. ROMANOV : "Problème inverse pour l'équation des ondes". Dakl Acad. Sci. URSS, 1973, Tom 2 11, n° 5.

[10] - P.C. WUENSCHEL : "Seismogram synthesis including multiples and transmission coefficients". Geophysics, vol. XX, n° 1 (February 1960), pp. 106-129.