La formulation classique

La formulation classique

La mise en œuvre pratique

 Nous venons de présenter la formulation mathématique conventionnelle des problèmes direct et inverse associés à la tomographie de temps de trajet de première arrivée. Nous nous intéressons maintenant à la mise en oeuvre pratique d’un algorithme de tomographie des temps de première arrivée défini à partir de cette formulation. Nous présentons dans cette partie les méthodes classiquement utilisées pour le calcul des temps de première arrivée t, pour la construction de la matrice des dérivées de Fréchet Ln et pour la résolution du système linéaire qui permet la détermination de la perturbation de lenteur δsn. 

La discrétisation du modèle de vitesse

 Le choix de la discrétisation retenue pour décrire le modèle du vitesse dépend de l’application considérée et des algorithmes mis en œuvre. Dans le cadre de la tomographie des temps de première arrivée, le modèle de vitesse est classiquement discrétisé sur une grille cartésienne. Cette discrétisation permet la prise en compte de variations latérales de vitesse et elle ne nécessite aucune connaissance a priori du milieu, contrairement à la discrétisation par couches où un nombre donné de couches est considéré. De plus, elle est totalement adaptée aux algorithmes mis en œuvre pour le calcul des temps de première arrivée et pour la résolution du problème inverse.

Le calcul des temps de première arrivée

Le calcul des temps de première arrivée repose sur une résolution numérique de l’équation eikonale (2.9). L’approche la plus couramment utilisée est la résolution par différences finies introduite par Vidale [1988]. Les algorithmes ainsi définis calculent les temps de première arrivée sur des grilles cartésiennes, à partir d’un modèle de vitesse lui aussi défini sur une grille cartésienne avec le même espacement de grille. L’implémentation d’un solveur eikonal peut être décomposée en deux problèmes : • un problème local, qui a pour objectif de définir le schéma aux différences finies utilisé pour calculer le temps de première arrivée en un point donné de la grille en fonction des valeurs de temps des points voisins ; • un problème global, qui permet de déterminer dans quel ordre les points de la grille doivent être considérés. L’article de Vidale [1988] a rapidement été suivi par des dizaines d’articles introduisant des variantes de schémas aux différences finies développés pour améliorer la rapidité, la robustesse, la stabilité, ou la précision, par exemple [van Trier and Symes, 1991], [Podvin and Lecomte, 1991], [Qin et al., 1992], [Pica, 1997]. De nombreux articles continuent d’être publiés chaque année sur le sujet donnant ainsi naissance à des nouveaux noms de méthodes telles que la Fast Marching Method [Sethian, 1996], [Popovici and Sethian, 1998], ou la Fast Sweeping Method [Boué and Dupuis, 1999], [Zhao, 2005]. Des variantes de ces méthodes permettent aussi de prendre en compte par exemple l’anisotropie du milieu [Lecomte, 1993], [Qian et al., 2007] ou bien encore de réaliser une implémentation parallèle de l’algorithme [Zhao, 2006]. 

Le tracé de rais a posteriori 

La résolution par différences finies de l’équation eikonale fournit une grille des temps de première arrivée en tout point du modèle. Il est alors ensuite possible de tracer les rais qui relient la source aux récepteurs de manière a posteriori [Podvin and Lecomte, 1991], [Baina, 1998]. Une méthode de calcul de la trajectoire des rais a posteriori repose sur le principe de stationnarité proposé par Vidale [1988]. On part du récepteur, pour une meilleure stabilité numérique, et on avance itérativement vers la source en suivant la direction opposée au gradient des temps. Ce gradient est estimé localement par un schéma aux différences finies centré sur un nœud du réseau ou bien sur le centre de la maille selon leur proximité de la position du point courant. Le tracé de rais a posteriori réalisé pour toutes les positions des couples {source – récepteur}permet alors de calculer la matrice des dérivées de Fréchet Ln. On calcule pour cela en chaque maille du modèle (Fig. 2.4) la longueur du segment du rai qui la traverse. 

La résolution du système linéaire tomographique

 La détermination de la perturbation de lenteur δsn passe par la résolution au sens des moindres carrés du système linéaire tomographique (2.29). Plusieurs méthodes de résolution de système algébrique ou d’estimation de l’inverse généralisé correspondant existent. On peut en distinguer deux classes : les méthodes de résolution directe comme la factorisation de Cholesky, la décomposition LU ou la décomposition par valeurs singulières et les méthodes de résolution itérative comme les méthodes de projection et de reconstruction. Les caractéristiques de la matrice des dérivées de Fréchet Ln sont sa grande taille [nombre de données × nombre de mailles], le fait qu’elle soit creuse, c’est-à-dire qu’elle contient essentiellement des éléments nuls, et son mauvais conditionnement, du à la couverture inhomogène des rais. Ces caractéristiques font que les méthodes itératives sont les mieux adaptées à la résolution du système tomographique [van der Sluis and van der Vorst, 1987], [Baina, 1998]. Les méthodes de reconstruction Les méthodes de reconstruction ont pour principe général de rétropropager à chaque itération les résidus des temps sur les mailles traversées par les rais calculés pour un modèle donné. Cette rétropropagation est pondérée par des coefficients liés à la couverture des rais dans chaque maille. La méthode de reconstruction la plus répandue est la méthode Simultaneaous Iterative Reconstruction Technique (SIRT). van der Sluis et van der Vorst [1987] ont prouvé théoriquement la convergence de cette méthode vers une solution au sens des moindres carrés pondérés. A partir de cette formulation commune, de nombreuses variantes ont été développées pour améliorer la rapidité de convergence [Zelt and Barton, 1998] ou bien prendre en compte les zones de Fresnel [Watanabe et al., 1999], [Grandjean and Sage, 2004]. L’avantage principal des méthodes de reconstruction est leur simplicité de mise en œuvre. Par contre, ces méthodes souffrent d’une renormalisation intrinsèque qui modifie le système linéaire à résoudre et donc la solution obtenue [van der Sluis and van der Vorst, 1987]. Les méthodes de projection La méthode de projection la plus utilisée est la méthode Least Squares QR (LSQR) développée par Paige et Saunders [1982]. Dans son principe, cette méthode ressemble à la fois à la méthode itérative du gradient conjugué et à l’inversion par décomposition par valeurs singulières, de façon à tirer profit des avantages des deux méthodes simultanément [Baina, 1998]. La méthode LSQR, en calculant à chaque itération des estimateurs de la valeur du conditionnement, permet de prévenir les instabilités numériques. Nolet [1985] a montré, lors d’une étude comparative, que l’algorithme LSQR est supérieur aux algorithmes de reconstruction, tant du point de vue de la vitesse de convergence que de la stabilité numérique. De nombreux algorithmes de tomographie des temps de première arrivée utilisent l’algorithme LSQR pour déterminer la perturbation de lenteur δsn, par exemple [Spakman and Nolet, 1988], [Le Meur, 1994], [Baina, 1998], [Zelt and Barton, 1998]. 

Schéma de l’algorithme de tomographie

 Le diagramme (Fig. 2.5) reprend les principales étapes algorithmiques aboutissant à la détermination de la perturbation de lenteur δsn pour une itération de la méthode de Gauss-Newton. A partir d’un modèle de lenteur courant sn, la résolution du problème direct pour N points de tirs fournit les cartes des temps de première arrivée en tous points de la grille. Pour toutes les positions {source – récepteur}, un tracé de rais a posteriori permet la construction de la matrice des dérivées de Fréchet. Cette matrice et les résidus calculés à partir des temps observés à la surface sont ensuite utilisés pour la résolution itérative du système linéaire tomographique. Cette résolution permet de déterminer la perturbation de lenteur δsn à appliquer au modèle courant.