Problèmes d’Optimisation Non Linéaire avec Contraintes en Tomographie de Réflexion 3D

Seismic reflection tomography allows the determination of a subsurface model from the traveltimes of seismic waves reflecting on geologic interfaces. The solution model of this inverse problem is a model that minimizes a least-squares functional measuring the mismatch between observed traveltimes and traveltimes calculated by raytracing in this model. My thesis subject consists in studying the resolution of this non linear minimization problem. The main difficulties come first from the size of the data space (around 50000 traveltimes) and of the model space (around 50000 unknowns), secondly from the conditioning problems due to the structure of the model vector (containing different sorts of unknowns : velocities in the geologic layers, transmission interfaces or reflecting interfaces). Thirdly, the non linearity of the modelling operator of traveltimes and the complexity of wave propagation in subsurface often imply an ill-posed minimization problem (Delprat-Jannaud, F., and Lailly, P., 1993). To ensure its well-posedness, curvature regularization is used. In addition, introduction of geological a priori information thanks to constrained optimization reduces uncertainties on the model solution.A classical line search Gauss-Newton method was used to solve unconstrained optimization problems in seismic reflection tomography (Chauvier, L., Masson, R., and Sinoquet, D., 2000). At each Gauss-Newton step, the (approximate) solution of a quadratic model of the objective function is computed using a conjugate gradient algorithm. This method generally allows to solve the minimization problem, but in some cases we noticed that the line search method is neither able to find a solution nor to decrease significantly the objective function. This difficulty is likely to be due to a too ill-conditioned Hessian matrix. To overcome this, we have proposed and implemented an original trust-region Gauss-Newton method coupled with a truncated conjugate gradient algorithm (Delbos, F., Sinoquet, D., Gilbert, J. C., and Masson, R., 2001) where the trust-region parameters (for example the trust-region radius) are automatically tuned by the solver. This method provides more stable and accurate results. In fact, it better solves intricate examples where the Hessian matrix is strongly ill-conditioned.The second part of my thesis consists in studying constrained optimization problems in seismic reflection tomography. The constraints we want to introduce in the optimization problem are of multiple types. They could be non linear (for example we could constrain the impact points of the rays on one interface to be located in a particular area) but in a first approach we prefer to limit ourselves with linear constraints. Even if the linearity brings simplifications, our optimization problems are still very difficult to solve because of the huge number of constraints (around 10000) and their different types.A widely used method to solve constrained non linear optimization problem is an interior point approach. However, we suspect this method to require significantly more resolutions of the forward problem, a raytracing solver, which is a drawback when this one is expensive in CPU time. Hence, we preferred to choose a Sequential Quadratic Programming (SQP) method because we had supposed that, as (Gauss-) newtonnian method, it does not require more iterations than the Gauss-Newton algorithm in the unconstrained case. This assumption proved to be exact. In the context of seismic tomography problems, a version of the augmented Lagrangian (AL) algorithm has been proposed by Glowinski and Tran (1993) to solve the tangent quadratic problem (TQP) of the SQP method. The solver QPAL we have developped in this thesis takes inspiration from that work and goes further by improving the efficiency of its augmented Lagrangian QP solver : it uses the classical multiplier method of Hesteness (1969) and Powell (1969) to minimize the dual fonction and it uses a GP-AC-GC algorithm to solve the Lagrange problem of the AL method. We have shown on a various concrete inversion that our nonlinear optimization method gives nice results :• the number of Gauss-Newton iterations has the same order of magnitude than the number of Gauss-Newton iterations necessary in the unconstrained case (around 10 iterations),• our nonlinear optimization method is efficient even for a large number of constraints,• the introduction of geological a priori information thanks to constrained optimization reduces uncertainties on the model solution.
La tomographie de réflexion de données sismiques permet la d´détermination d’un modèle de sous-sol à partir des temps de trajet des ondes sismiques se réfléchissant sur les interfaces géologiques. Le modèle solution de ce problème inverse est le modèle qui minimise une fonctionnelle moindres-carrés mesurant les écarts entre les temps de trajet mesurés et les temps de trajet calculés par tracé de rayons dans ce modèle. Le sujet de ma thèse est l’étude de la résolution numérique de ce problème non linéaire de minimisation.Les principales difficultés proviennent tout d’abord de la taille de l’espace des données (de l’ordre de 500000) et de l’espace des modèles (de l’ordre de 50000 inconnues), puis des problèmes de conditionnement liés à la structure du vecteur modèle (différents types d’inconnues interviennent : les vitesses dans les couches géologiques, les interfaces en transmission ou en réflexion). Enfin, la non linéarité de l’opérateur de modélisation des temps de trajet et la complexité de la propagation des ondes dans le sous-sol conduisent souvent à un problème de minimisation mal posé. Une régularisation sur la courbure du modèle permet d’avoir un problème de minimisation mathématique bien posé. De plus, l’introduction d’informations géologiques a priori par résolution sous contraintes du problème inverse permet de réduire les incertitudes sur le modèle solution.Classiquement, une méthode de Gauss-Newton couplée avec une méthode de recherche linéaire était utilisée pour résoudre les problèmes d’optimisation sans contrainte en tomographie de réflexion sismique (Chauvier, L., Masson, R., and Sinoquet, D., 2000). A chaque itération de Gauss-Newton, une solution (approchée) d’un modèle quadratique de la fonction coût non linéaire est calculée en utilisant l’algorithme du gradient conjugué. Cette méthode permet en général de résoudre les problèmes de minimisation traités en tomographie. Cependant, nous avons remarqué pour certains cas que la méthode de recherche linéaire était non seulement incapable de trouver une solution, mais encore de faire décroitre de manière significative la fonction coût. Ces difficultés sont principalement dues à une matrice Hessienne mal conditionnée. Afin de pallier à ces difficultés, lapremière partie de ma thèse a consisté à remplacer la recherche linéaire par une méthode de région de confiance, et à développer un algorithme de gradient conjugué tronqué pour minimiser le modèle quadratique de la fonction coût non linéaire. Les différents paramètres de la méthode des régions de confiance (par exemple le rayon de la région de confiance) sont automatiquement réglés par le solveur. Cette méthode donne des résultats plus stables et plus précis. Elle permet de mieux résoudre les exemples complexes où le Hessien est très mal conditionné.La deuxième partie de ma thèse consiste à étudier la résolution des problèmes d’optimisation avec contraintes en tomographie de réflexion sismique. Les contraintes à ajouter dans la résolution de nos problèmes d’optimisation peuvent être de types très différents. Elles peuvent être non linéaires (nous pourrions par exemple contraindre les points d’impact des rayons sur une interface à être localisés dans une région particulière), mais dans une première approche nous préférons nous limiter à des contraintes linéaires. Même si la linéarité apporte des simplifications, nos problèmes d’optimisation restent encore très difficiles à résoudre notamment à cause du nombre élevé de contraintes (de l’ordre de 10000) et de leurs types différents.Actuellement, une méthode souvent utilisée pour résoudre les problèmes d’optimisation non linéaire avec contraintes est l’approche par points intérieurs. Cependant, nous soupçonnons que cette méthode demande significativement plus de résolutions du problème direct, le tracé de rayon, ce qui est un inconvénient lorsque celui-ci coûte cher en temps CPU. Ainsi, nous avons préféré une méthode de Programmation Quadratique Successive (PQS) car nous supposions que, comme méthode (Gauss-) Newtonnienne, elle ne demanderait pas plus d’itérations (donc d’évaluation de fonctions), que l’algorithme de Gauss-Newton dans le cas sans contrainte. Cette précision s’est avérée exacte. Glowinski et Tran, en 1993, ont développé et testé une méthode PQS dans laquelle le problème quadratique tangent (PQT) est résolu via une méthode de Lagrangien Augmenté. Nous nous sommes inspirés de cette première approche en améliorant la résolution des PQT. Le code QPAL que nous avons développé pendant cette thèse a montré son efficacité pour résoudre des problèmes quadratiques convexes de grandes tailles : il utilise d’une part la méthode classique des multiplicateurs de Hesteness et Powell pour minimiser la fonction duale et d’autre part l’algorithme GP-AC-GC pour résoudre les problèmes de Lagrange. Nous avons montré sur de nombreux cas tests ou réels que cette méthode donne des résultats convaincants :• le nombre d’évaluations de la fonction coût reste du même ordre de grandeur que celui obtenu par le solveur Gauss-Newtonnien sans contrainte,• cette méthode est capable de résoudre des problèmes de grande taille,• l’introduction d’informations a priori par l’inversion avec contraintes permet de réduire les incertitudes sur le modèle solution.