Brevet FR2691811A1 — publié le 19931203

PROCEDE DE DETERMINATION D'UN MODELE DE DISTRIBUTION OPTIMALE DES VITESSES DE PROPAGATION D'ONDES SISMIQUES DANS LE SOUS-SOL. La présente invention concerne un procédé pour la détermination d'un modèle de distribution optimale des vitesses de propagation d'ondes sismiques dans le soussol. Des ondes sismiques se propageant à travers un terrain donnent naissance à des signaux de sortie qui caractérisent le terrain et représentent en quelque sorte sa signature. Ces signaux de sortie peuvent donc théoriquement servir à déterminer la distribution des vitesses à travers le sous-sol et à connaître à partir de là, la composition des différentes couches géologiques. Le procédé qui consiste à trouver la distribution des vitesses à partir des signaux sismiques s'appelle l'inversion sismique. Comme l'inversion sismique est un problème fortement non-linéaire, il n'existe pas en général de solution analytique. Par conséquent, la solution de ce problème doit être déterminée par des méthodes numériques. La méthode la plus simple pour effectuer une inversion sismique est celle qui consiste à tester tous les signaux d'entrée possibles, et à observer les sorties qu'ils engendrent. La réponse modélisée qui s' approche le plus de la réponse observée est considérée comme étant la solution de l'inversion. Mais cette méthode est trop lente et trop coûteuse, comme on va l'expliquer sur un exemple. La distribution des vitesses peut être représentée par une matrice dont chaque élément est égal à la vitesse mesurée en une pluralité de points équidistants formant un réseau s 'étendant en profondeur et en surface. En supposant que les vitesses sont mesurées sur une profondeur de 5 km tous les 25 mètres, soit en 200 points, et sur une longueur en surface de 20 km tous les 200 mètres, soit en 100 points, la matrice représentant un modèle de vitesses aura 20 000 éléments, ce qui représente un échantillonnage grossier pour une surface tellement étendue. La modélisation par inversion consisterait à essayer pour chacun des éléments de la matrice, plusieurs valeurs, par exemple 10 valeurs, ce qui donnerait un nombre total de modèles égal à 10200000.Aucun ordinateur actuel n'est capable de traiter tant de modèles en un temps raisonnable. Une autre méthode couramment utilisée pour résoudre le problème de l'inversion sismique est celle dite du gradient. On en expliquera le principe en regard des figures 1 et 2 annexées. Sur la figure 1, la référence 10 désigne le vrai système dont on veut déterminer la distribution des vitesses, x désigne l'entrée réelle et y la sortie observée. Le but est de retrouver la valeur de x. La figure 2 montre un système 12 qui est un modèle de 10 et dans lequel on introduit une entrée de test xi. A la sortie du système 12 on recueille une sortie modélisée De la même façon, pour toutes les entrées modélisées possibles xl, x2...xk, on peut mesurer les sorties modélisées correspondantes y1, Y2,...yk. On peut ainsi calculer pour chaque valeur yak, la valeur de l'écart IIY-Ykll entre l'observation réelle y et l'observation due à la modélisation yak. Cet écart sera désigné par la suite par distance ou coût C. On pourra donc tracer la courbe C (xk) I tY-Ykl I Naturellement, le modèle correct xk est celui qui conduit à un coût nul ou minimum.Le but de la méthode du gradient est de trouver une suite de xk qui diminue successivement le coût jusqu a ce qu'on arrive à un minimum de la fonction C. De cette manière, on évite de tester toutes les entrées possibles. Pour cela, la méthode du gradient calcule la tangente à la courbe C par rapport à l'entrée xk. Comme cette tangente indique la direction de la pente la plus forte de la fonction du coût, on peut toujours trouver un xk+l pour lequel C(Xk+l)4 C(xk) Le procédé du gradient peut être appliqué, de façon connue en soi, au problème de l'inversion sismique. On l'expliquera en regard de la figure 3. Celle-ci représente une installation classique d'acquisition de données sismiques. L'installation comprend un émetteur d'ondes sismiques 14 monté sur la surface de la terre et une pluralité de récepteurs 16 disséminés sur la surface de la terre. L'émetteur envoie un train d'ondes incidentes 18 qui est réfléchi vers la surface de la terre par des réflecteurs 20 constitués par les interfaces entre les différentes couches géologiques.Les ondes sismiques réfléchies 22,24 sont enregistrées sous forme de traces sismiques correspondantes 26, 28. Si les réflecteurs sont horizontaux et la vitesse constante, les traces sont des portions d'hyperbole, mais dans le cas général, on obtient des enregistrements différents de l'hyperbole, que l'on appelera indicatrices. Ici également, on doit trouver par le procédé de l'inversion, le modèle de distribution des vitesses qui explique les traces sismiques enregistrées par les récepteurs 16. Comme précédemment, la méthode du gradient appliquée au problème de l'inversion sismique conduit à définir une fonction de coût C, c'est-à-dire une mesure de la différence entre les données sismiques calculées ou modélisées Pcal et les données sismiques observées PObs en fonction de la vitesse v. La fonction de coût s'exprime par Toutefois, la méthode du gradient appliquée à l'inversion sismique pose un certain nombre de problèmes - comme on le sait, l'inversion sismique est non linéaire.Il en résulte que la courbe du coût présente de nombreux minima locaux où la tangente est horizontale, de sorte que le minimum que la méthode du gradient permet de trouver dépend du choix de l'entrée initiale. Or, ce que l'on recherche c'est le minimum global, c'est-à-dire le plus petit de tous les minima locaux. Un exemple de génération d'un minimum local pour l'inversion sismique est illustré à la figure 4 ci-jointe où l'on voit une indicatrice modélisée 30 et une indicatrice observée 32 qui ne coïncident pas. On comprend qu'il n'existe pas dans ce cas de perturbation du modèle de vitesses qui puisse modifier la fonction de coût qui, on le rappelle, est égale à la différence entre les données observées et les données modélisées. Le gradient est donc égal à zéro et on se retrouve dans un minimum local.Pourtant, il est clair que ce minimum n1 est pas global, ce qui arrive quand les deux indicatrices se superposent - les composantes à haute fréquence du modèle de vitesses contribuent fortement au calcul du gradient et les composantes à faible fréquence n'y contribuent que faiblement, ce qui signifie que la méthode du gradient appliquée à l'inversion converge rapidement pour les composantes à haute fréquence et lentement pour les composantes à faible fréquence. Or, ce sont ces dernières qui contribuent le plus à la solution du problème de l'inversion sismique. La présente invention a pour objet de remédier aux inconvénients des méthodes exposées ci-dessus, et concerne un procédé de détermination d'un modèle de distribution des vitesses de propagation d'ondes sismiques dans un terrain, basé sur la migration et la remodélisation des données sismiques et utilisant l'équation d'onde paraxiale. La migration sismique donne des images de réflecteurs en fonction du modèle estimé de distribution des vitesses. La fonction de coût est formulée de telle manière qu'elle ne dépende pas de ces réflecteurs imagés. I1 en résulte que la méthode de l'inversion sera plus sensible aux composantes à basse fréquence du modèle de vitesses. De plus, même si la fonction de coût dépend des amplitudes des signaux sismiques, le problème de 1 ' inversion est formulé de telle manière que l'indicatrice observée et l'indicatrice modélisée coincident presque toujours, quelle que soit l'estimation qui est faite pour le modèle de vitesses. On peut donc supposer que la fonction de coût selon cette méthode, aura moins de minima locaux et se prêtera mieux à l'inversion. Avant de décrire la nouvelle formulation de l'inversion sismique basée sur l'équation d'onde paraxiale, on rappelera ci-après ce qu'est cette équation. Une équation d'onde paraxiale convient idéalement pour la migration de l'hyperbole modélisée du fait qu'elle ne concerne pas les réflexions se dirigeant vers le haut. Pour cette raison, on l'appelle aussi équation d'onde "à un sens", par opposition à l'équation d'onde acoustique qui est une onde à deux sens puisqu'elle fait intervenir l'onde réfléchie et l'onde transmise par un réflecteur. L'onde réfléchie peut donner naissance à des relations temporelles avec la propagation vers le bas des données sismiques, et produire ainsi des artéfactes dans les images migrées. On comprend donc qu'une équation à un sens est exactement l'outil qu'il faut pour effectuer la migration sismique. On expliquera à présent le développement mathématique de l'approximation paraxiale et la raison pour laquelle elle conduit à une équation à un sens. L'équation acoustique à deux sens x,z est de la forme où p est le champ d'ondes et v est la vitesse de propagation. En lui appliquant la transformation de Fourier directe, cette équation devient d'où l'on tire l'équation de dispersion La représentation graphique de kzv/w en fonction de kxv/ est un cercle 38 de rayon 1, comme représenté à la figure 5 ci-jointe. D'après les explications données précédemment, ce que lton désire c'est une relation de dispersion à un sens qui ne comprenne que le demi-cercle inférieur 40 représenté à la figure 6, et qui concerne donc uniquement l'énergie qui se propage vers le bas, de manière à éviter les réflexions. La relation de dispersion de ce demi-cercle est Malheureusement, on ne sait pas calculer la transformée de Fourier inverse de ce demi-cercle.Mais, par un artifice de calcul, on peut faire une bonne approximation de l'équation (1) en la remplaçant par une équation polynomiale rationnelle de la forme La figure 7 ci-jointe montre sur un même graphique la représentation graphique des équations (1) et (2) pour un choix des coefficients an, bm donnant Comme on le voit sur cette figure, l'approximation du demi-cercle par l'équation (2) est bonne pour La transformée de Fourier inverse de l'équation (2) est illustrée sur la figure 8 par 42. Elle correspond à une source de Dirac située à l'origine du plan distance horizontale-distance verticale x-z et elle a la forme d'une cardiolde se trouvant dans la moitié inférieure du cercle 44 qui correspond à la transformée de l'équation d'onde acoustique.La figure 8 montre bien que cette approximation de l'équation d'onde est vraiment à un sens puisque l'énergie de la source de Dirac ne se propage que dans la direction des z négatifs. Les branches voûtées de la cardioïde sont des artéfactes mais ne sont pas gênantes parce qu'elles peuvent être filtrées. Ayant ainsi expliqué l'équation d'onde paraxiale, on exposera à présent le procédé selon l'invention. La présente invention concerne un procédé de détermination d'un modèle de distribution optimale des vitesses de propagation d'ondes sismiques dans un sous-sol qui comprend au moins deux couches séparées par un interface, caractérisé en ce qu'il consiste - à émettre des ondes sismiques à partir d'une source placée successivement en au moins deux points de tir de la surface du sol ; - à enregistrer pour chacune des positions de la source les traces des ondes sismiques réfléchies par l'interface entre les formations géologiques sous-jacentes à l'aide de récepteurs disposés en une pluralité d'emplacements à la surface du sol - à imaginer à partir des traces obtenues un modèle initial de distribution de vitesses à travers le sous-sol;; - à migrer les points de tir sismiques par l'équation paraxiale en utilisant ce modèle de distribution de vitesses afin de créer une image sismique de l'intérieur du sous-sol - à modéliser chaque point de tir avec l'équation paraxiale en propageant en profondeur une onde due à une source placée en surface du sol - à remodéliser les données sismiques de chaque point de tir avec équation paraxiale en extrapolant une onde vers la surface, initialisée (ou déclenchée) par l'interaction de l'onde provenant de la source et l'image sismique - à comparer les données sismiques enregistrées de chaque point de tir avec les données remodélisées de chaque point de tir et, dans le cas où ils ne coincident pas, à calculer le gradient de la fonction de coût par rapport au modèle de vitesses initiales - à modifier le modèle de vitesses avec le gradient afin de diminuer la différence entre les données sismiques observées et remodélisées - et à itérer de cette manière jusqu a ce que le gradient soit zéro ou suffisamment faible pour satisfaire les besoins de l'inversion, ce qui définit le modèle optimal de vitesses. Ainsi, selon l'invention, on introduit une information supplémentaire qui facilite grandement la réussite de l'inversion, à savoir que les données sismiques obéissent au principe de la cohérence géométrique. En d'autres termes, lorsqu'un modèle de vitesses correct est utilisé, et que l'on déplace la source d'ondes sismiques, les zones de l'interface qui sont imagées à la fois par deux points de tir doivent coïncider. Ce principe est basé sur l'hypothèse logique de la nature géométrique invariable de la terre : la terre est supposée ne pas changer d'une expérience sismique à une autre. Le procédé selon l'invention sera expliqué dans le cas où les données sismiques proviennent d'une source que l'on déplace en deux points différents, mais il s'applique également au cas où l'on fait un grand nombre d'expériences. Lorsque le modèle de vitesses imaginé est incorrect, le réflecteur est imagé différemment par les deux expériences et l'on obtient une image incohérente formée de deux courbes croisées. On peut utiliser cette image incohérente pour remodéliser les données de la seconde position de la source. Pour cela, on fait propager vers le bas les ondes sismiques émises par la source placée au second point de tir, puis on génère des ondes se dirigeant vers le haut, chaque fois que l'onde descendante touche une partie de l'image. Cette étape de modélisation génère des données observées à la surface de la terre qui contiennent à la fois les données originales de la source placée au second point de tir (puisque la migration et la modélisation sont des procédés inverses l'un de l'autre) et les parties de l'image générée par la source placée au premier point de tir qui ne coincident pas avec 1 image formée par la source placée au second point de tir. On applique ensuite à ce résultat le procédé de l'inversion afin de minimiser la différence entre les données du second point de tir et les données remodélisées de ce point de tir obtenues après avoir imagé les points de tir avec le modèle de vitesses. Par rapport aux procédés de l'inversion exposés précédemment, à savoir l'inversion de la forme d'onde, l'inversion selon l'invention présente les avantages suivants - elle tire partie de la cohérence géométrique du terrain, - même si cette inversion est basée sur la comparaison des formes d'onde, une confidence entre les données observées et les données modélisées est toujours assurée, étant donné que les données remodélisées sont générées à partir du même modèle de vitesses que celui que l'on utilisait pour imager les réflecteurs - étant donné que les réflecteurs sont imagés et ne constituent pas une partie du modèle de vitesses, l'inversion ne nécessite pas d'inverser ladite partie, et permet plus facilement de récupérer les composantes à fréquence basse et intermédiaire. Pour pouvoir effectuer l'inversion, il est nécessaire de calculer le gradient du coût par rapport au modèle de vitesses. Pour cela, le problème est formulé au moyen des équations de Lagrange suivantes où M(p,v) est l'opérateur paraxial, L(p,v) est le Lagrangien défini en fonction de p et v, k représente le champ d'ondes, pO les conditions initiales du champ d'ondes, v le modèle de vitesses et C le coût. Le développement de Lagrange du gradient fait apparaître un opérateur adjoint X et une variable associée q. Ce développement s' exprime comme suit où Zf est la profondeur d'investigation finale. On décrira à présent l'invention en regard des dessins annexés dans lesquels Les figures 9 et 10 sont des vues schématiques en élévation expliquant le principe de la cohérence géométrique La figure 11 est une vue schématique en élévation montrant l'influence de la cohérence géométrique sur les traces d'ondes sismiques modélisées ou observées La figure 12 est un exemple de modèle de vitesses à trois points de tir La figure 13 représente les données fournies par un dispositif de modélisation utilisant le champ de vitesses de la figure 12 Les figures 14 et 15 montrent respectivement les courbes d'évolution du coût et du gradient en fonction du nombre d 'itérations de l'algorithme d inversion du gradient La figure 16 montre l'image migrée des trois points de tir en fonction du nombre d'itérations de l'algorithme d'inversion ; et La figure 17 montre les trois points de tir remodélisés, en haut, après la première itération, au milieu, tels qu'observés à l'origine, et en bas, après la centième itération. Sur les figures 9 et 10 on distingue un interface 44 entre deux couches de terrain. La source placée en un premier point de tir 14 (figure 9) émet un faisceau incident d'ondes sismiques 18 qui se réfléchit sur l'interface dont elle illumine une zone d'impact 46. La source est ensuite déplacée en un second point de tir 14' (figure 10) et son faisceau incident 18' illumine une zone impact 48 décalée par rapport à la zone d'impact 46. Les deux zones d'impact ont une partie d'illumination commune 50. Sur la figure 11, on distingue les deux images 52,54 de l'interface migrées à partir des traces enregistrées aux points de tir 14 et 14' respectivement. Les données remodélisées pour le point de tir 14' sont enregistrées par les récepteurs 16 et montrent les réflections dues aux images 52, 54. La réfection 58 correspond aux données observées du point de tir 14', étant la suite d'une migration et d'une modélisation des mêmes données. Par contre, la réflection 56 correspond aux données du point de tir 14 ramenées au point de tir 14' On se référera à présent à l'exemple de réalisation illustré par la figure 12. Les données synthétiques sont générées par un programme de modélisation acoustique utilisant le modèle de vitesses tel que montré à la figure 12. Ce modèle comprend deux couches 60,62 séparées par un interface 64 pourvu d'une dépression 66.La couche supérieure du modèle présente un gradient de vitesse allant de 1600 m.s -1 à 2200 m.s -1, tandis que la couche inférieure a une vitesse constante de 3500 m.s Avec ce modèle, on a utilisé trois points de tir 14,14',14" disposés comme sur la figure 12 et soixante quatre récepteurs, non représentés, associés à chaque point de tir. Les récepteurs ont été distribués à des intervalles égaux sur toute la surface du modèle. Les données générées par le dispositif de modélisation acoustique sont représentées à la figure 13. Les indicatrices relatives aux trois points de tir ont des amplitudes qui varient en raison de la présence de la dépression et du gradient de vitesse latéral. Sur les figures 14 et 15, l'itération initiale de l'algorithme d'inversion utilise un champ de vitesse homogène de 1500 m.s-l. On constate bien sur la figure 14 que le coût diminue en fonction des itérations et que le gradient tend vers zéro, ce qui indique que le coût converge vers un minimum (qui peut être local). Sur l'image migrée de la figure 16, les itérations retenues sont les itérations 1,4,7,9,13,16,19,22 et 100, l'itération 1 étant montrée dans le coin inférieur gauche de la figure 18 et l'itération 100 étant visible dans le coin supérieur droit de la figure. Comme on le constate sur cette figure, la cohérence de l'image migrée s'améliore à mesure que l'algorithme d'inversion converge. On note également que la dépression du réflecteur est de mieux en mieux imagée. Enfin sur la figure 17, on peut noter qu'après la première itération, l'image migrée contient des indicatrices fantômes qui ne correspondent pas aux indicatrices observées, mais qu'après la centième itération, les indicatrices fantômes ont disparu. L'invention a permis de proposer un nouveau procédé d'inversion basé sur l'équation d'onde paraxiale qui établit une correspondance entre un point de tir et un autre point de tir et qui utilise des informations géométriques sur la nature fixe du terrain, dans le processus d'inversion. REVENDICATION Procédé de détermination d'un modèle de distribution optimale des vitesses de propagation d'ondes sismiques dans un sous-sol, qui comprend au moins deux couches séparées par un interface, caractérisé en ce qu'il consiste - à émettre des ondes sismiques à partir d'une source placée successivement en au moins deux points de tir (14,14') de la surface du sol - à enregistrer pour chacune des positions de la source les traces des ondes sismiques réfléchies par l'interface entre les formations géologiques sous-jacentes à l'aide de récepteurs (16) disposés en une pluralité d'emplacements à la surface du sol - à imaginer à partir des traces obtenues un modèle initial de distribution de vitesses à travers le soussol - à migrer les points de tir sismiques (14,14') par une équation paraxiale en utilisant ce modèle de distribution de vitesses afin de créer une image sismique de l'intérieur du sous-sol - à modéliser chaque point de tir avec l'équation paraxiale en propageant en profondeur une onde due à une source placée en surface du sol - à remodéliser les données sismiques de chaque point de tir avec l'équation paraxiale en extrapolant une onde vers la surface, initialisée (ou déclenchée) par l'interaction de l'onde provenant de la source et l'image sismique - à comparer les données sismiques enregistrées (56) de chaque point de tir avec les données remodélisées (58) de chaque point de tir et, dans le cas où ils ne coincident pas, à calculer le gradient de la fonction de coût par rapport au modèle de vitesses initiales - à modifier le modèle de vitesses avec le gradient afin de diminuer la différence entre les données sismiques observées et remodélisées - et à itérer de cette manière jusqu a ce que le gradient soit zéro ou suffisamment faible pour satisfaire les besoins de l'inversion, ce qui définit le modèle optimal de vitesses.