FR2843202A1

FR2843202A1 - Methode pour former un modele representatif de la distribution d'une grandeur physique dans une zone souterraine, affranchi de l'effet de bruits correles entachant des donnees d'exploration

Info

Publication number: FR2843202A1
Application number: FR0209919A
Authority: FR
Inventors: Patrick Lailly; Francois Renard; Laure Pelle; Jannaud Florence Delprat
Original assignee: IFP Energies Nouvelles IFPEN
Current assignee: IFP Energies Nouvelles IFPEN
Priority date: 2002-08-05
Filing date: 2002-08-05
Publication date: 2004-02-06
Anticipated expiration: 2022-08-05
Also published as: CN1503006A; GB2391665B; US7092823B2; GB2391665A; NO328835B1; CN1333267C; GB0318295D0; NO20033459D0; US20050075791A1; FR2843202B1

Abstract

- Méthode pour former, à partir de données obtenues par exploration d'une zone d'un milieu hétérogène (une zone souterraine par exemple), un modèle représentatif de la distribution dans la zone d'au moins une grandeur physique largement affranchi de la présence de bruits corrélés pouvant entacher les données.- La méthode comporte essentiellement les étapes suivantes :- acquisition de données renseignant sur certaines caractéristiques de la zone ;- spécification d'un opérateur de modélisation qui, à un modèle, associe la réponse du modèle (données synthétiques) ;- modélisation de chaque bruit corrélé par application, à une fonction génératrice du bruit (f.g.b.), d'un opérateur de modélisation spécifique ;- spécification d'une semi-norme dans l'espace des données ;- spécification, dans chaque espace de f.g.b., d'une norme pour laquelle chaque opérateur de modélisation du bruit constitue une isométrie ; et- recherche, par une méthode algorithmique tirant partie de ces propriétés d'isométrie, du modèle et des f.g.b. qui minimisent une fonction coût quantifiant, par l'intermédiaire de la semi-norme, l'écart entre les données mesurées et la superposition de la réponse du modèle et des bruits associés aux f.g.b.- Applications par exemple à la recherche de la distribution dans une zone souterraine, de l'impédance acoustique, de vitesses de propagation, de perméabilités, etc.

Description

La présente invention concerne une méthode pour former, à partir de

données obtenues par exploration d'une zone d'un milieu hétérogène, un modèle représentatif de la distribution dans la zone d'une grandeur physique affranchi, (au moins partiellement) de la

présence de bruits corrélés pouvant entacher les données.

La méthode s'applique par exemple à la quantification de l'impédance acoustique

dans une zone souterraine.

ETAT DE LA TECHNIQUE

La démarche consistant à rechercher un modèle qui s'ajuste sur des mesures expérimentales, s'est développée dans quasiment tous les domaines des sciences ou de la technologie. Une telle approche est connue sous différentes appellations: méthode des moindres carrés pour l' estimation de paramètres, pour la résolution du résolution du problème inverse. Pour une bonne présentation de cette approche dans le contexte des géosciences, on peut se référer par exemple à: - Tarantola, A.: " Inverse Problem Theory: Method for Data Fitting and Model Parameter Estimation ". Elsevier, Amsterdam, 1987; Notons que le terme " moindres carrés " fait référence au carré de la norme dans l'espace des données pour quantifier l'écart entre la réponse d'un modèle (laquelle est l'image du modèle par un opérateur de modélisation préalablement choisi) et les données, fonction cot qu'il s'agit de minimiser pour résoudre le problème. L'utilisation du carré de la norme pour définir la fonction cot ne constitue qu'une commodité pratique mais n'a fondamentalement rien d'essentiel. D'ailleurs de nombreux auteurs utilisent, pour des raisons diverses, une autre définition de la fonction cot mais cette définition reste basée sur l'exploitation de la norme, ou d'une semi-norme, dans l'espace des données. Enfin nous disposons d'une grande latitude pour choisir la norme (ou la semi-norme) dans l'espace des données (nous ne sommes en aucune façon contraints à utiliser la norme euclidienne). Dans le cas de données bruitées la solution peut dépendre sensiblement du choix fait à ce stade. Pour de plus amples développements au sujet de ce problème, on peut se référer par exemple aux publications suivantes Tarantola, A.: " Inverse Problem Theory: Method for Data Fitting and Model Parameter Estimation ". Elsevier, Amsterdam, 1987; Renard et Lailly, 2001; Scales et

Gersztenkorn, 1988; Al-Chalabu, 1992.

Les résultats de mesures sont souvent entachés d'erreurs. A ces erreurs de mesures se rajoutent, quand on veut confronter l'expérimentation avec des résultats de modélisation, des bruits de modélisation: les modélisations ne sont jamais parfaites et correspondent donc toujours à une vue simplifiée dé la réalité. Nous ferons donc apparaître dans la suite le bruit comme la superposition - de composantes non corrélées (bruit blanc par exemple) - de composantes corrélées, ce qui veut dire que l'existence de bruit sur un échantillon de mesure se traduira par l'existence d'un bruit de même nature sur certains des points de mesure voisins; typiquement les bruits de modélisation font partie de cette

catégorie.

Lorsque des bruits corrélés entachent les données, la qualité du modèle estimé par résolution du problème inverse peut s'en trouver gravement affectée. Comme nous l'avons dit aucun opérateur de modélisation n'est parfait. C'est donc toute la communauté des personnes impliquées dans l'identification de paramètres décrivant un modèle qui voit sa démarche entravée par l'existence de bruits corrélés. Parmi ces personnes, les praticiens de l'exploration sismique figurent parmi les plus concernés: leurs données ont en effet un rapport signal sur bruit médiocre voire très mauvais. Voilà pourquoi les techniques de filtrage des bruits corrélés constituent une partie importante des logiciels de traitement de données sismiques. Les techniques les plus classiques exploitent une transformation (transformation de Fourier par exemple) o signal et bruit se trouvent localisés dans des régions de l'espace différentes permettant ainsi une séparation du signal et du bruit. Pour un présentation générale des méthodes classiques d'élimination du bruit sur des données sismiques, on peut se référer au livre de Yilmaz: - Yilmaz, O. 1987: "Seismic data processing", Investigation in geophysics N02, Society

of Exploration Geophysicists, Tulsa, 1987.

Ces techniques ne sont cependant pas parfaites: elles supposent que l'on a trouvé une transformation permettant la séparation complète du signal et du bruit. Dans ce contexte, les bruits corrélés apparaissent comme particulièrement gênants car ils peuvent être difficiles à séparer du signal (qui, lui aussi, est corrélé) et peuvent être de fortes amplitudes. On est ainsi souvent conduit à la difficile gestion du compromis suivant: ou bien on préserve le signal mais il reste un fort résidu de bruit, on bien on élimine le bruit mais on a, du coup, distordu le signal. Ces techniques de filtrage peuvent être mises en oeuvre préalablement à la résolution du problème inverse: elles constituent alors un prétraitement appliqué au données. La qualité de la solution du problème inverse dépend

alors fortement de l'aptitude des filtres à éliminer le bruit sans distordre le signal.

L'approche introduite par Nemeth et al dans la publication suivante: Nemeth T. et al (1999), "Least-Square Migration of Incomplete reflection Data"; Geophysics, 64, 208-221, constitue une avancée importante pour l'inversion de données perturbées par des bruits corrélés de forte amplitude. Ces auteurs proposent d'éliminer le bruit par résolution d'un problème inverse: ayant défini l'espace des bruits corrélés comme l'espace image d'un espace vectoriel B (espace de fonctions génératrices de bruit) par un opérateur linéaire T et le signal comme l'espace image d'un espace vectoriel M (espace des modèles) par un opérateur de modélisation F (qu'ils supposent linéaire), ils recherchent le signal dans F(M) et le bruit dans T(B) dont la somme est la plus proche possible (au sens de la norme euclidienne ou, en continu, de la norme L2) de la donnée mesurée. Cette technique constitue une avancée importante dans la mesure o elle a permis l'élimination (ou du moins une réduction très substantielle) de bruits corrélés de forte amplitude (ondes de surface par exemple) difficilement séparables du signal par les techniques conventionnelles de filtrages. Cependant, de l'aveu même des auteurs, ces performances se payent par une augmentation d'un ordre de grandeur du temps de calcul requis pour la résolution du problème inverse classique c'est à dire sans rechercher la composante corrélée du bruit. Par ailleurs le résultat fourni par la méthode s'avère extrêmement sensible à toute inexactitude introduite lors de la définition de l'opérateur T: c'est là la contrepartie inéluctable de la

grande aptitude de la méthode à discriminer signal et bruit.

LA METHODE SELON L'INVENTION

La méthode selon l'invention permet d'estimer, à partir de données obtenues par exploration d'une zone souterraine, un modèle représentatif de la distribution, dans la zone, d'au moins une grandeur physique, ce modèle étant largement affranchi de la présence de bruits corrélés pouvant entacher les données. La méthode comporte essentiellement les étapes suivantes a) acquisition de mesures renseignant sur certaines caractéristiques physiques de la zone en suivant un protocole expérimental défini; b) spécification d'un opérateur de modélisation, qui, à un modèle de chaque grandeur physique associe des données synthétiques constituant la réponse du modèle, les mesures et les données synthétiques appartenant à un espace de données; c) pour chaque bruit corrélé repéré par un indice j allant de 1 à J, choix d'un opérateur de modélisation qui associe un bruit corrélé à une fonction génératrice du bruit appartenant à un espace préalablement défini de fonctions génératrices du bruit; d) spécification d'une norme ou d'une semi-norme dans l'espace des données; e) spécification, d'une semi-norme dans l'espace de fonctions génératrices de bruit pour laquelle chaque opérateur de modélisation du bruit établit sensiblement une relation isométrique entre l'espace des fonctions génératrices du bruit et l'espace des données; f) définition d'une fonction cot quantifiant l'écart entre les mesures d'une part et la superposition de la réponse du modèle et des bruits corrélés associés aux fonctions génératrices de bruit d'autre part; et g) ajustement du modèle et des fonctions génératrices du bruit en minimisant la fonction cot, par une méthode algorithmique tirant partie, des propriétés d'isométrie des

opérateurs de modélisation du bruit.

Suivant un premier mode de mise en oeuvre, on recherche la distribution en fonction de la profondeur de l'impédance acoustique dans le milieu, les bruits corrélés affectant les données sont des ondes de tube identifiées chacune par des paramètres caractérisant leur propagation, les données mesurées sont des données de type PSV en utilisant des capteurs adaptés à détecter le déplacement des particules dans le milieu en réponse à une excitation sismique localisée, l'emplacement des capteurs, la durée d'enregistrement et les points d'échantillonnage en temps étant définis, et l'opérateur de modélisation choisi associe les données synthétiques à une distribution d'impédance acoustique en fonction de la profondeur évaluée en temps de parcours et à la contrainte verticale mesurée en fonction

du temps à la profondeur du premier capteur.

Comme fonction cot quantifiant l'écart entre les mesures d'une part et la superposition de la réponse du modèle et des bruits corrélés associés aux fonctions génératrices de bruit d'autre part, on choisit par exemple le carré de la semi-norme de cet

écart dans l'espace des données.

Suivant un exemple de réalisation, l'ajustement du modèle et des fonctions génératrices du bruit est obtenu par exemple par une méthode de relaxation pour éliminer les inconnues correspondant à chaque fonction génératrice de bruits corrélés, cette méthode

de relaxation étant mise en oeuvre à l'intérieur des itérations d'un algorithme de quasiNewton pour le calcul du modèle.

Suivant un exemple de réalisation, le calcul numérique de l'image d'un modèle par cet opérateur est effectué par exemple par résolution numérique de l'équation des ondes à ID pour le modèle considéré, en sélectionnant des valeurs prises par le déplacement des particules aux emplacements de capteurs et aux points d'échantillonnage en temps préalablement spécifiés et en appliquant un opérateur propre à compenser les effets de

divergence sphérique et d'atténuation.

Suivant un exemple de réalisation, l'opérateur numérique de modélisation du bruit est un schéma numérique centré aux différences finies pour discrétiser l'équation de transport du bruit et la fonction génératrice de bruit intervenant comme condition initiale en bordure de la dite zone associe à chaque bruit corrélé, un espace de fonctions génératrices du bruit formé par les fonctions du temps à support dans un intervalle de

temps donné.

Différents exemples de normes ou semi-normes choisies pour l'espace des données

et l'espace des fonctions génératrices de bruit seront donnés dans la description détaillée

qui va suivre.

Suivant un deuxième mode de mise en oeuvre on recherche la distribution de perturbations, par rapport à un milieu de référence préalablement choisi, de l'impédance et de la vitesse dans la dite zone du milieu, les bruits corrélés affectant les données sont dus à des réflexions multiples dont la cinématique et les variations d'amplitude avec l'offset ont été préalablement estimés, les données mesurées sont captées par des capteurs sismiques de surface, l'emplacement des dits capteurs, le mode d'excitation sismique, la durée d'enregistrement et les points d'échantillonnage en temps étant définis, et l'opérateur de modélisation étant défini via une linéarisation de l'équation des ondes autour du milieu de référence. Comme fonction cot quantifiant l'écart entre les mesures d'une part et la superposition de la réponse du modèle et des bruits corrélés associés aux fonctions génératrices de bruit d'autre part, on choisit par exemple le carré de la semi-norme de cet

écart dans l'espace des données.

Suivant un exemple de réalisation, l'ajustement du modèle et des fonctions génératrices du bruit est obtenu par une méthode de relaxation par bloc pour éliminer les inconnues correspondant à chaque fonction génératrice de bruits corrélés, cette méthode de relaxation étant mise en oeuvre à l'intérieur des itérations d'un algorithme de gradient

conjugué pour le calcul du modèle.

PRESENTATION SUCCINCTE DES FIGURES

Les caractéristiques et avantages de la méthode selon l'invention, apparaîtront plus

clairement à la lecture de la description ci-après d'un exemple non limitatif de réalisation,

en se référant aux dessins annexés o: - la figure 1 montre un exemple de données obtenues par mise en oeuvre d'une méthode de prospection sismique verticale (PSV); - la figure 2 montre un modèle de distribution variation d'impédance acoustique avec la profondeur; - la figure 3 montre des données de PSV synthétiques obtenues sur la base du modèle d'impédance de la Fig.2; - la figure 4 montre un exemple de données de PSV contaminées par un seul bruit corrélé; - la figure 5 montre un exemple de données de PSV contaminées par deux bruits corrélés; - la figure 6 montre la distribution d'impédance en fonction de la profondeur, obtenue par inversion des données bruitées de la Fig.4; - la figure 7 montre les résidus d'inversion (écarts entre les données de la Fig. 4 et la réponse sismique du modèle de la Fig. 6); - la figure 8 montre la distribution d'impédance en fonction de la profondeur, obtenue par inversion des données bruitées de la Fig.5; - la figure 9 montre les résidus d'inversion correspondants (écarts entre les données de la Fig. 5 et la réponse sismique du modèle de la Fig.8) - la figure 10 montre la distribution d'impédance en fonction de la profondeur, obtenue par inversion des données bruitées de la Fig.4 en recherchant le bruit corrélé sous la forme de la superposition de deux bruits corrélés ayant des propriétés de propagation erronées c'est à dire différentes de celles du bruit apparaissant sur les données de la Fig. 4; - la figure 11 montre les résidus d'inversion correspondants d'inversion (écarts entre les données de la Fig. 4 et la réponse sismique du modèle de la Fig. 10) - les figures 12A, 12B montrent respectivement des données sismiques avec réflexions multiples, et la réponse sismique du modèle obtenu après inversion linéarisée classique; - les figures 13A, 13B montrent respectivement un exemple de distribution d'impédance et de vitesse de propagation auquel les données sismiques de la Fig.12A constituent la réponse sismique; et les figures 14A, 14B montrent respectivement la comparaison des réponses sismiques des modèles obtenus par inversion classique (14A) et par inversion suivant la technique proposée incluant un pointé du "moveout" du multiple et une estimation des variations

d'amplitude avec l'offset (14B).

DESCRIPTION DETAILLEE DE LA METHODE

La formulation du problème On dispose de données qui proviennent d'une mesure échantillonnée d'une

fonction. Cette fonction dépend de plusieurs variables (variables d'espace ou d'espacetemps par exemple).

Ces données, notées d, correspondent à des mesures effectuées pour nous renseigner sur un modèle m. Elles se trouvent entachées de différents bruits: - des bruits corrélés (ou une superposition de bruits corrélés);

- des bruits non corrélés.

Le problème est de déterminer quantitativement le modèle m (ou des fonctions de

ce modèle) à partir de la donnée d.

Pour ce faire nous choisissons un opérateur de modélisation F (linéaire ou non) qui, au modèle m, associe la réponse du modèle. Cet opérateur ne modélise en fait qu'imparfaitement le phénomène physique réel. C'est d'ailleurs essentiellement (mais pas exclusivement) pour cette raison qu'apparaissent, dans les données, des bruits corrélés: ces bruits correspondent à un signal lié au modèle mais la relation entre l'un et l'autre nous apparaît trop complexe pour être incluse dans l'opérateur de modélisation F que, pour des

raisons diverses, nous voulons garder relativement simple.

Dans ces conditions, nous modélisons le bruit en introduisant un ou plusieurs espaces générateurs de bruit Bj ( variant de 1 à J) et les opérateurs de modélisation de bruits associés que nous appelons Tj. Pour retrouver le modèle m à partir des données d,

nous proposons de rechercher ce modèle comme solution du problème d'optimisation.

min C(mfi11,i2.../ f)J F(m)+YTj,1/3j-d (1) ( P)eMxlB1 J=I D j=L

O D est une norme (ou éventuellement une semi-norme) dans l'espace des données.

Le formalisme ci-dessus se distingue de l'approche proposée par Nemeth et al (1999) de par le fait que: - l'opérateur de modélisation F n'est pas nécessairement linéaire; - nous nous autorisons à considérer une superposition de bruits corrélés de types différents (en ce sens qu'ils sont modélisés par des opérateurs différents référencés par l'indice j), ceci au prix d'une inflation des calculs tant en ce qui concerne le nombre d'inconnues (ce qui complique singulièrement la résolution du problème d'optimisation) qu'en ce qui concerne le nombre de modélisations de bruit à effectuer; - nous nous réservons la possibilité de choisir à notre convenance la norme il I|D ou

éventuellement la semi-norme, dans l'espace des données.

Le choix que nous allons faire pour i| ||D va être dicté par des considérations liées à l'efficacité de la méthode de résolution (décrite plus bas), cette efficacité permettant de traiter le cas d'une superposition de bruits corrélés. Il s'agit là d'une ouverture importante: outre la possibilité de traiter des données entachées de bruits corrélés aux caractéristiques complexes, on peut par ce biais surmonter les difficultés, rencontrées dans l'approche de Nemeth et al., liées à la grande sensibilité du résultat à une imprécision introduite lors de la définition de l'opérateur de modélisation du bruit: il suffira, pour surmonter cette difficulté, de rechercher le bruit par superposition de bruits associés à des opérateurs de

modélisation Tj approximatifs.

Notons que le choix du carré de la norme I ID intervenant dans le problème d'optimisation n'a rien d'essentiel: nous ne changeons pas la solution en choisissant à la place une fonction qui, comme la fonction carré, est une fonction croissante sur 91' (ou une

fonction décroissante à condition de changer min en max).

La méthode de résolution La méthode consiste à choisir convenablement: la norme (ou semi-norme) I Ji, dans l'espace des données;

- des normes || 1| dans chaque espace Bj.

de façon à permettre la mise en place d'une solution algorithmique très performante au problème d'optimisation (1). Par "choisir convenablement" nous entendons faire en sorte que chaque opérateur Tj constitue une isométrie (ou une approximation d'isométrie) pour

les normes Iij et i [||D respectivement dans l'espace de départ et dans l'espace d'arrivée.

L'idée directrice dans la détermination de normes permettant de rendre chaque opérateur Tj isométrique repose sur l'existence d'égalités de type à conservation d'énergie vérifiées par les solutions des équations aux dérivées partielles (ou d'égalités d'énergie discrètes pour les solutions de certains schémas numériques utilisés pour la discrétisation de ces équations). Si nous sommes capables de faire de tels choix, il est alors possible de rendre très simple la minimisation de la fonction cot en éliminant les inconnues associées aux fonctions génératrices de bruits corrélés (détermination du complément de Schur). Cette élimination est réalisée à l'aide d'une algorithmique adaptée telle que la méthode de relaxation par blocs, un bloc étant associé à une fonction génératrice de bruits, ou la méthode du gradient conjugué avec un préconditionnement de Gauss Seidel par bloc symétrique, etc., méthodes de mise en oeuvre algorithmique toutes bien connues des gens de l'art, et que l'on trouvera présentées par: - Golub G., H.et al (1983); "Résolution numérique des grands systèmes linéaires" (Eyrolles) La minimisation simplifiée de la fonction cot permet de réduire considérablement

le temps nécessaire à la résolution numérique.

Premier exemple de mise en oeuvre appliqué à l'inversion à ID de données PSV contaminées par une onde de tube Les données PSV (Profil Sismique Vertical) font partie des acquisitions classiques de sismique de puits. Elles ont des applications multiples, l'une des principales étant leur

utilisation pour faire le lien entre données sismiques de surface et mesures diagraphiques.

L'inversion de données PSV telle que décrite par exemple par: - Macé, D. et Lailly, P. (1984) A Solution of an Inverse problem with ID Wave

Equation applied to the Inversion of Vertical Seismic Profile; Proc. of the 6th Intern.

conf. in "Analysis and Optimisation of Systems"; Nice, 2, 309-323; a été développée à cette fin. Dans le problème inverse ID, on suppose le modèle de subsurface invariant latéralement et une excitation par ondes planes se propageant verticalement. Les inconnues du problème sont la distribution d'impédance acoustique en fonction de la profondeur (mesurée en temps de parcours vertical et sur un intervalle démarrant à la profondeur du premier récepteur) ainsi que le mode d'excitation sismique (de façon précise la fonction du temps caractérisant la condition aux limites à la profondeur du premier capteur). Il s'agit là d'un problème inverse non linéaire, les données PSV

dépendant de façon non linéaire de la distribution d'impédance.

Les données PSV sont souvent contaminées par un bruit corrélé de très forte

amplitude: l'onde de boue. Il s'agit d'une onde guidée dans la boue envahissant le puits.

Sa vitesse de propagation est lente par rapport à la vitesse de propagation dans les roches à l'entour du puits. De plus la vitesse de propagation de cette onde dépend de la fréquence (propagation dispersive). Enfin cette onde se réfléchit notamment au fond du puits, donnant de plus naissance à un autre mode de propagation. Des données PSV typiques sont données sur la figure ci-dessous: on y voit les principales caractéristiques de l'onde de tube. Son amplitude et l'extension de la zone contaminée rendent difficilement exploitable

le signal contenu dans les nombreux échantillons sismiques contaminés par l'onde de boue.

Nous décrivons ci-après la mise en oeuvre de la méthode pour l'inversion à ID des

données PSV.

Le protocole expérimental Il consiste en la spécification: - des différentes cotes auxquelles les récepteurs seront placés dans le puits (supposé vertical): dans nos expériences les capteurs couvrent l'intervalle de profondeurs mesurées en temps simple [xm,,jn= O.5s; xmna> = 0.2s], avec un récepteur tous les Ax=lms (en temps simple); nous obtenons ainsi de échantillons numérotés de 0 à I; - de la nature physique de la mesure effectuée par ces capteurs: dans nos expériences les capteurs mesurent le déplacement vertical résultant de la propagation des ondes; de la durée des enregistrements: dans nos expériences les capteurs mesurent l'état vibratoire sur l'intervalle de temps [O,T=1.5s], ces données étant échantillonnées tous les At=lms. Nous obtenons ainsi des échantillons numérotés de 0 à N. Nous appelons di' l'échantillon de données enregistré à la profondeur x,,,in+ i Ax et

au temps n At. Nous avons bien sr: T= N At et Xmax= Xmin+I AX.

Acquisition des données conformes à ce protocole expérimental Nos expériences ont été faites à partir de données synthétiques c'est-à-dire obtenues par modélisation numérique, elle même réalisée en deux étapes: Création de données synthétiques non bruitées: celle-ci a été réalisée par résolution numérique de l'équation des ondes acoustiques à 1D, le modèle étant la distribution d'impédance acoustique représentée sur la Fig. 2 (la profondeur est mesurée en temps simple) et le mode d'excitation (condition aux limites de Neumann à la profondeur 0) étant

une classique ondelette de Ricker.

Les données PSV synthétiques issues de cette modélisation sont représentées dans la figure 3: l'axe vertical représentant le temps d'observation est gradué de 0 à 1.5 s et l'axe horizontal représentant la profondeur (en temps de parcours simple) des différents

capteurs est gradué de 0.05 à 0.2 s.

Adjonction de bruit sur les données Nous allons perturber les données PSV ainsi calculées en leur ajoutant d'une part un bruit aléatoire (mais filtré dans la bande des fréquences sismiques) d'amplitude assez forte et d'autre part un ou plusieurs bruits corrélés d'amplitude très forte. Chacun de ces bruits se propagera vers le bas avec une vitesse lente qui dépend de la fréquence (propagation dispersive): ces bruits corrélés sont censés représenter des ondes de boue. Ils ont été modélisés par résolution numérique (utilisation du schéma centré aux différences finies) d'une équation de transport avec une vitesse de propagation constante. La dispersion des ondes a été obtenue en utilisant dans de grands pas de discrétisation dans le schéma aux différences. Dans le cas de la superposition de plusieurs bruits corrélés, les vitesses de propagation associées à chacun de ces bruits sont différentes. Les deux figures 4, 5 montrent les données PSV contaminées par un seul bruit corrélé (Fig.4) et par la

superposition de deux bruits corrélés (Fig.5).

Spécification de l'opérateur de modélisation La modélisation est réalisée par la résolution de l'équation des ondes: CF (_) -0) (X ) - O dans]xr.2, +cc[ x [0, T] avec la condition aulx limites de Neumann: - C7 (X mrain) (x t) = h(t) sur [O, T] et les conditions initiales y (x,) = t(x) = Osr sr v [ Ot La définition de l'opérateur de modélisation demande la définition préalable de l'opérateur d'observation qui formalise le protocole expérimental décrit au 1. L'opérateur d'observation est donc l'opérateur qui, à la fonction y(x,t) solution du système ci-dessus,

associe les échantillons y(xi,t1) pour xi=xmin+i Ax et tn=n At, i=O,...,I; n=O,...,N.

L'opérateur de modélisation que nous appellerons F(m) est donc l'opérateur (non linéaire) qui au couple de fonctions m=(o(x),h(t)) associe les échantillons y(xi,t') pour

xi=xni+i Ax et t1=n At, i=O,...,j; n=O,....,N.

Spécification de l'opérateur de modélisation associé à chaque bruit corrélé Nous précisons dans ce paragraphe la procédure pouvant être utilisée pour modéliser chaque bruit corrélé. Chaque bruit corrélé sera caractérisé par des directions de corrélation qui lui seront spécifiques et que nous supposons connues, au moins approximativement: ces directions de corrélation sont spécifiées par l'intermédiaire d'un champ de vecteurs de corrélation 3-j(x,t) de composantes (cjx(xt) cjt(x,t)) qui devra être défini en tout point du domaine [Xln,,, x]X[O,T]. Nous pouvons de façon équivalente spécifier les directions de corrélation par l'intermédiaire d'un faisceau de lignes de corrélation qui représentera les lignes dechamp associées à 3j (x, t). De telles lignes de corrélation peuvent être obtenues suivant la technique décrite dans les brevets EP 354.112 (US 4 972 383) du demandeur, à partir d'un pointé de quelques phases caractéristiques du bruit l'information étant étendue à tout le domaine XiiXmax j]X[O,T] en utilisant des procédures classiques d'interpolation et/ou d'extrapolation. La spécification des directions de corrélation revient à spécifier la distribution de vitesse de propagation du bruit cj(x,t) C (X, t) qui est liée au vecteur de corrélation par la relation c, (x,t)= c (xt) Nous supposons ici qu'il n'y a pas de variations d'amplitude de l'onde lors de sa propagation (une situation différente sera présentée pour le 2 me exemple d'application). La modélisation d'un bruit corrélé tel que spécifié ci-dessus est basée sur la remarque suivante: une onde se propageant le long des lignes de corrélation sans changer d'amplitude est solution de l'équation de transport: eb(Xt) c(x:t) = O Dans ce contexte, nous introduisons - un espace Bj de fonctions génératrices de bruit qui sera l'espace des fonctions /j (t) à support dans Fçin, t.axlc [o,T] et que nous prolongeons par 0 sur l'intervalle [o, T] entier; - l'opérateur de modélisation du bruit Tj: solution de I 'équation de transport ebj(X, t) = A(Xt) = i dan m.[ Xma1 x [O, T] et vérifiant, les conditions initiales 4 ( tXI) O b(X Zmin t) d (t) Il faut remarquer que la géométrie du faisceau de lignes de corrélation (ou, ce qui revient au même, la distribution de vitesse de propagation) ne peut pas être quelconque: les lignes de corrélation ne sont pas autorisées à se couper: sinon le problème de résolution de l'équation de transport serait un problème mal posé. Pour les mêmes raisons les lignes s de corrélation ne sont pas autorisées à intercepter deux fois les bords sur lesquels les conditions initiales sont spécifiées. Cette dernière considération peut amener (pour modéliser un bruit corrélé correspondant à une onde montante), à spécifier la condition initiale non plus en x=x mais en x=xr,, ou même en une valeur intermédiaire quitte à décomposer le problème de Cauchy en deux problèmes associés aux domaines situés de

part et d'autre de ce bord intermédiaire.

Nous supposons par ailleurs que la géométrie des lignes de corrélation est telle que (des situations plus générales peuvent être envisagées mais les résultats qui suivent prennent alors une forme différente) çj (Xz t) = (X) x Â(t) ou. encore: cJ(Xit) = ^^ (T) et. -(fxJt) =Ai Enfin nous supposons que le faisceau de lignes de corrélation qui interceptent en x=xmin l'intervalle [tr it7 axJ] n'interceptent pas en t=T l'intervalle [x... Xmax] En fait l'équation de transport est résolue par l'intermédiaire d'un' schéma

numérique. On peut par exemple utiliser le classique schéma centré aux différences finies.

Pour ce faire nous choisissons des pas de discrétisation Axj' et Atj' (qui ne sont pas nécessairement égaux à Ax et At mais que nous supposerons, même si des situations plus générales peuvent être envisagées, sous multiples de ces quantités) et introduisons un maillage dont les noeuds sont les points de coordonnées (x1n.in+ i'Axj',n'/\tj') avec i' et n'E N. Le classique schéma centré aux différences finies explicité ci-dessous permet, partant des conditions initiales, le calcul de proche en proche des différentes valeurs prises par la fonction bj(x,t) aux différents noeuds du maillage pour simplifier les notations nous

avons omis l'indice j associé au bruit corrélé considéré).

(t+2 tn +L + ni + - - bn + 2,At+ ' [bii + 1;+1 -h + -il b]= Dans la formule ci-dessus la quantité oj'1 représente l'évaluation de la quantité cx (cx est un terme générique) au point de coordonnées (xmin+ i'Axj', n'Atj') Bien entendu l'utilisation de pas de discrétisation Axa' et Atj' suffisamment petits devant la longueur d'onde est requise si nous voulons que le schéma numérique donne une

approximation précise de la fonction solution de l'équation de transport.

Cependant, l'utilisation de pas de discrétisation plus grands (mais toujours plus petits que la longueur d'onde) permet de modéliser des phénomènes propagatifs plus complexes en rendant la vitesse de propagation de l'onde dépendante de la fréquence: on

peut ainsi modéliser des propagations dispersives telles la propagation des ondes de tube.

Le choix de pas de discrétisation adaptés à la modélisation d'un mode particulier de l'onde de tube peut être réalisé: En déterminant, par un examen des données, la vitesse de propagation de ce mode (qui dépend ici non seulement de l'espace et du temps mais encore de la fréquence); des techniques de tomographie peuvent être utilisées à cette fin, telles que celles présentées dans la publication suivante: - Ernst F. et al (2000); Tomographt of Dispersive Media; J. Acoustic. Soc. Am., 108, -116. En connaissant les propriétés de propagation dispersive associées au schéma numérique, propriétés qui découlent d'une analyse de la relation de dispersion numérique (comme décrit par exemple dans les publications suivantes: - Alford R. M. et al, (1974), Accurency of Finite Difference Modeling of the Acoustic

Wave Equation; Geophysics, 6, 834-842.

Spécification d'une norme ou d'une semi-norme dans l'espace des données Un premier choix consiste à prendre comme norme dans l'espace des données une norme discrète quelconque qui se veut une approximation (via une formule de quadrature) de la norme définie en continu par: III (LITdz dt 1u12 (XI t)) Nous pouvons aussi, ce qui constitue, comme nous le verrons au suivant, un meilleur choix, utiliser la semi-norme: / I-L jN-1 2 1I1V11D a (UX+ t + bey + (2) o u est un vecteur quelconque de l'espace des données (i et n sont les indices représentant

les numéros d'échantillons respectivement en espace et en temps).

Par ailleurs d'autres choix sont possibles.

Spécification, dans chaque espace Bj de fonctions génératrices du bruit, d'une norme pour laquelle chaque opérateur de modélisation du bruit Tj établit une relation isométrique (ou l'approximation d'une relation isométrique) entre l'espace des fonctions génératrices du bruit et l'espace des données (muni de la semi-norme définie en 5) Nous précisons dans ce paragraphe la procédure pouvant être utilisée pour définir la norme dans l'espace des fonctions génératrices associé à chaque bruit corrélé. Cette procédure étant la même pour chaque bruit, nous la décrivons de façon générique en

omettant, dans les formules, l'indice j associé au bruit considéré.

Un premier choix consiste à munir chaque espace B de fonctions génératrices de bruit d'une norme discrète quelconque qui se voudrait une approximation (via une formule de quadrature) de la norme définie en continu par: IlIIIB = (X jT dt rifj(t)() auquel cas, en utilisant le fait que, grâce aux hypothèses énoncées au 4, on a bj(x,T)=O pour xe [xr nXmax], l'opérateur linéaire Tj réalise une isométrie entre Bj et D. Mais, comme le calcul de T1 /j est effectué numériquement le caractère isométrique de l'opérateur T ne sera pas parfait mais seulement approché, l'approximation étant d'autant meilleure que les pas de discrétisation Ax, At, Ax ',At 'sont petits. Le fait que nous n'ayons alors qu'approximativement une isométrie se traduira par des performances moins bonnes

au niveau de la mise en oeuvre algorithmique présentée au paragraphe 8 ciaprès.

Un meilleur choix consiste à munir chaque espace Bj de fonctions génératrices de bruit de la serai norme: IIBIIBs = (AÄLS At I >E '+ 1 + 1) ) (3) Dans ce cas l'opérateur Tj réalise rigoureusement une isométrie entre Bj et D, du moins si u N = O Vi = 0, I -1. Cela résulte de l'égalité d'énergie discrète:

JN-i I-1. fa-.

i Z + L(7 +710] = E E n [( 2]

+ E E aAt [7k+1 + L)]-

i=0 i0= égalité valable dans le cas (simplifié pour les besoins de la présentation) O A.ô = zAx et. At' = At et. si l'on a défini: _= Définition de la fonction cot Nous définissons la fonction cot par la formule: C(m,,,,61.. 16J-bc F (m) +ô Tj,6j - d j=l D Recherche du modèle et des fonctions génératrices du bruit minimisant la fonction cot, cette recherche se faisant par une méthode algorithmique tirant parti, explicitement

ou implicitement, des propriétés d'isométrie des opérateurs de modélisation du bruit (cf.6).

La méthode algorithmique exploite la spécificité du complément de Schur, spécificité associée au caractère isométrique des opérateurs Tj. Ceci peut se faire par exemple en réalisant que minimiser C(m,,/3, /B2,...- 3) revient à minimiser C'(m) = C(m, /3 (m), /2 (m),-...e (m)) o (1 (m), / 2 (m),...A (m)) minimise C(m, 36,/2,-,...) pour m donné. Notons que la hessienne associée à cette forme quadratique est le complément de Schur. La minimisation de C'(m) peut être réalisée par une quelconque méthode d'optimisation telle que la version BFGS d'une méthode de

Quasi Newton.

La détermination de (>S1 (m), 2 (m),... (m)) est très aisée grâce au caractère isométrique des opérateurs Tj., notamment si nous avons effectués les choix (2) et (3) respectivement pour définir 11 |D et i1 il. S'il n'y a qu'un seul type de bruit (J=1) et si l'isométrie est parfaite, la détermination de /31 (m) demande simplement une évaluation du gradient de la forme quadratique. S'il y a superposition de plusieurs bruits (J>1), différents algorithmes peuvent être envisagés pour profiter du caractère isométrique des opérateurs Tj. On peut par exemple déterminer (/3k(m), /2(m),.... J(m)) par une méthode de relaxation (alias Gauss-Seidel) par bloc, un bloc étant associé à un groupe d'inconnues Taj (m), une itération de relaxation ne demandant simplement (toujours dans le cas d'une isométrie parfaite) qu'une seule évaluation du gradient de la forme quadratique associée au bloc considéré. Pour les problèmes rendus difficiles par la proximité des directions de corrélation des différents types de bruits, on peut envisager une méthode de gradient conjugué préconditionné, la matrice de préconditionnement étant de type Gauss-Seidel par bloc (alias relaxation) symétrique. Là encore la mise en oeuvre du préconditionnement est très facile, la résolution du problème associé à un bloc ne demandant qu'une évaluation du

gradient de la forme quadratique associée au bloc considéré.

Les Fig.6 à 11 montrent les résultats obtenus sur trois types d'expériences de mise en oeuvre de la méthode. Dans le cas d'une inversion de données bruitées (Fig.4) à la fois par la superposition d'un bruit corrélé et d'un bruit aléatoire, la Fig.6 donne la distribution d'impédance

retrouvée, tandis que la Fig.7 montre les résidus d'inversion.

Dans le cas d'une inversion de données bruitées (Fig.5) à la fois par la superposition de deux bruits corrélés et d'un bruit aléatoire, la Fig.8 donne la distribution d'impédance

retrouvée, tandis que la Fig.9 montre les résidus d'inversion.

Dans le cas d'une inversion de données bruitées (Fig.4) à la fois par la superposition d'un bruit corrélé et d'un bruit aléatoire o l'on ne connaît pas de façon précise les propriétés de propagation de l'onde de boue, il est naturel de rechercher un modèle de bruit corrélé constitué de la superposition de deux bruits corrélés ayant des vitesse de propagation erronées mais encadrant la vraie vitesse de propagation de l'onde de boue. La Fig. 10 donne la distribution d'impédance retrouvée, alors que la Fig. 11 montre les résidus d'inversion. Deuxième exemple d'application de la méthode à la détermination d'un modèle par inversion linéarisée de données multi-déports (multi-offsets) contaminées par des réflexions multiples Cette application se distingue de la première principalement par le fait qu'elle illustre les possibilités qu'offre la méthode vis à vis de la prise en compte de variations d'amplitude du bruit le long des directions de corrélation. Une autre différence réside dans les propriétés de propagation des bruits corrélés, la propagation des réflexions multiples

n'étant pas dispersive.

L'inversion linéarisée figure parmi les méthodes sophistiquées utilisées par les géophysiciens pour obtenir un modèle quantitatif des hétérogénéités de la subsurface, ce

type d'information étant très précieux notamment pour la caractérisation des réservoirs.

Les données sont ici des données sismiques de surface. Cette méthode est décrite par exemple dans la publication suivante: - Bourgeois A. et al., (1989) "The Linearized Seismic Inverse Problem: an Attractive Method for a Sharp Descriptiopn of Strategic Traps": 59th Ann. Internat. Mtg. Soc. Expl. Geoph. Expanded Abstract, pages 973-976. Par ailleurs il est nécessaire de disposer d'un modèle de référence, constitué (dans le contexte d'un modèle acoustique) d'une distribution de vitesse et d'une distribution d'impédance acoustique: ce sera le modèle autour duquel on effectuera une linéarisation dans l'équation des ondes (approximation de Born). Pour simplifier la présentation nous nous placerons dans le cas ID mais cette hypothèse est d'autant moins essentielle que les applications pratiques concernent pour l'essentiel la sismique 3D. Dans le contexte 1D, toute l'information sismique est contenue dans un seul point de tir. Par ailleurs le modèle de référence est lui aussi iD: il ne dépend que de la profondeur. Les inconnues du problème sont la distribution de perturbation d'impédance acoustique et de vitesse en fonction de la profondeur. En revanche, contrairement au premier exemple d'application le mode d'excitation sismique (et notamment rondelette sismique) seront supposés connus. fl

s'agit là d'un problème inverse rendu linéaire grâce à la linéarisation.

Les données sismiques de surface sont souvent contaminées par des réflexions multiples qui peuvent être d'assez forte amplitude. Leur cinématique, qui, dans le contexte iD, se caractérise par le "moveout", est en général assez différente de celle des réflexions primaires: ces deux types d'ondes ont, en général, traversé des couches géologiques ayant des vitesses de propagation différentes. Voilà pourquoi l'inversion linéarisée constitue en elle-même, un technique efficace pour atténuer les réflexions multiples: la cinématique des réflexions primaires (les événements modélisés par l'approximation de Born) est entièrement définie par la distribution de vitesse dans le modèle de référence: les événements modélisés ne peuvent donc s'ajuster sur des réflexions multiples pour peu que le "moveout" de ces dernières diffère de celui des primaires. Cependant, en raison des stationnarités à l'offset nul, il y a, en pratique un ajustement partiel et les résultats d'inversion linéarisée s'avèrent contaminés par la présence du bruit corrélé que constituent les réflexions multiples. Les données sismiques de la Fig.12A ont été modélisées avec prise en compte des réflexions multiples (il n'y a que 3 primaires qui arrivent à l'offset nul successivement aux temps:500; 1000; 1650 ms, la 2ème arrivée étant d'amplitude

particulièrement forte).

Sur la réponse sismique du modèle obtenu après inversion linéarisée classique (cf. Fig.12B), le multiple de forte amplitude arrivant vers 1500 ms à l'offset nul n'a été que faiblement atténué et le modèle se trouve contaminé par le multiple. Il est alors naturel de tenter d'exploiter les possibilités offertes par la méthode de Nemeth pour mieux discriminer multiples et primaires, notamment avec les améliorations apportées par notre méthode. Nous décrivons ci-après la mise en oeuvre de la méthode pour l'inversion linéarisée

à 1D de données sismiques de surface.

Le protocole expérimental Il consiste en la spécification: - Du mode d'excitation sismique: dans nos expériences la source sismique est un point source modulé en temps par une fonction (ondelette sismique) notée w(t); - des différents emplacements auxquels les récepteurs sont disposés: dans nos expériences les capteurs sont situés à la profondeur lOm et couvrent l'intervalle d'offsets [xmin=O,x,,,a=SOOm], avec un récepteur tous les Ax= 100m; nous obtenons ainsi des traces sismiques numérotées de 0 à 1; - les capteurs mesurent le champ de pression résultant de la propagation des ondes; - de la durée des enregistrements: dans nos expériences les capteurs mesurent l'état vibratoire sur l'intervalle de temps [0,T=1800ms], ces données étant échantillonnées tous les At=lms. Nous obtenons ainsi des échantillons numérotés de 0 à N. Nous appelons di' l'échantillon de données enregistré pour l'offset Xmin+ i Ax et au

temps n At. Nous avons bien sr: T= N At et Xma,= Xmin+I Ax.

Acquisition des données conformes à ce protocole expérimental Nos expériences ont été faites à partir de données synthétiques: ces données ont été obtenues par résolution numérique (méthode des différences finies) de l'équation des ondes 2D: -2ii -V-7 P = ô (., z-zs) w(t) dans e x a*+ x [0, T] avec la condition aulx limites: P(z O) O et les conditions initiales P(t= (}) = (t) =o équation dans laquelle: - x, z, t désignent respectivement la coordonnée latérale, la profondeur et le temps; - et c représentent respectivement les distributions (fonctions de la profondeur) d'impédance acoustique et de vitesse de propagation L'ondelette choisie est une classique ondelette de Ricker de fréquence centrale Hz. et les distributions d'impédance a et de vitesse c sont données sur les figure 13A, 13B. La réponse sismique de ce modèle est donc constituée des données représentées sur la Fig. 12A. Nous appelons FNL l'opérateur qui, au couple de fonctions (Y(z),c(z)), associe les

données sismiques.

Spécification de l'opérateur de modélisation Nous choisissons un modèle de référence décrit par le couple de fonctions (co(z),co(z)). Dans les expériences présentées dans la suite nous avons choisi co(z)=c(z) et Go(z) =cte (=1). Nous choisissons pour opérateur de modélisation l'opérateur jacobien de l'opérateur FNL au point (c0(z),co(z)). En fait, puisque nous avons choisi co(z)=c(z), la seule inconnue est 86(z)= 6(z)- 6o(z) et seule la composante correspondante du jacobien

est importante.

L'opérateur de modélisation que nous appellerons F(Èm) est donc l'opérateur (linéaire) qui au couple de fonctions âm=(&6(z),8c(z)) associe les échantillons Èy(xi,t') qui

sont les composantes du vecteur F (8m) pour xi=xm+i Ax et t'=n At, i=0,... ,I; n=0,..,N.

Spécification de l'opérateur de modélisation associé à chaque bruit corrélé Nous précisons dans ce paragraphe la procédure pouvant être utilisée pour modéliser chaque bruit corrélé. Chaque bruit corrélé sera caractérisé par des directions de corrélation qui lui seront spécifiques et que nous supposons connues, au moins approximativement: ces directions de corrélation sont spécifiées par l'intermédiaire d'un champ de vecteurs de corrélation 3j (x,t) de composantes (cjx(xt), cjt(x,t)) qui devra être défini en tout point du domaine [Xl,,,x.,i]X[O,T]. Nous pouvons de façon équivalente spécifier les directions de corrélation par l'intermédiaire d'un faisceau de lignes de corrélation qui représentera les lignes de champ associées à 3j(x,t). De telles lignes de corrélation peuvent être obtenues à partir d'un pointé de quelques phases caractéristiques du bruit l'information étant étendue à tout le domaine xm,,Xmax.]X[O,T] en utilisant des procédures classiques d'interpolation et/ou d'extrapolation, comme décrit dans le brevet EP 354.112 du demandeur déjà cité. Pour ce deuxième exemple d'application dans lequel nous voulons éliminer le multiple précédemment mentionné, nous pointons ce multiple ( par exemple le pic d'amplitude), ce qui nous permet de définir les variations de temps d'arrivée en fonction de l'offset et nous définissons un continuum de lignes de corrélation par une simple translation verticale de la ligne pointée (dans ces conditions le vecteur 3(x,t) ne dépend pas de t). La spécification des directions de corrélation revient à spécifier la distribution de vitesse de propagation du bruit cj(x,t) qui est liée au vecteur de corrélation par la relation ci (xt) = c (xt) t) Contrairement au ler exemple d'application déjà décrit, nous avons ici de variations d'amplitude du bruit le long de la ligne de corrélation. Nous supposons que nous pouvons estimer ces variations d'amplitude à partir d'une mesure sur les données ou bien à l'aide de

considérations théoriques: cette mesure nous défini une fonction g(x).

La modélisation du bruit se fera donc par la composition de deux opérateurs: Un opérateur de transport, comme dans le premier exemple d'application Une modulation d'amplitude, c'est à dire la multiplication de la solution de

l'équation de transport par la fonction g(x).

En reprenant le schéma donné pour le ler exemple d'application (et notamment en faisant appel aux mêmes hypothèses et notations), nous introduisons: - un espace Bj de fonctions génératrices de bruit qui sera l'espace des fonctions /j (t) à support dans [t, lt.ax]c [0,T] et que nous prolongeons par 0 sur l'intervalle [0,T] entier; l'opérateur de modélisation du bruit Tj: Tj:,?j (t) (E tj - bj (x,t) e D solhtion de l'équation de transport j (3, t) - È} (X, t) = O dans [ min Xrnrj.z] x [P, T] et vérifiant les conditions initiales: b b(3, D O) - (i -j u(z = Xm.i, t) =,fI(t) En fait l'équation de transport est résolue par l'intermédiaire d'un schéma

Pour ce faire nous choisissons des pas de discrétisation Axj' et Atj' (qui ne sont pas nécessairement égaux à Ax et At mais que nous supposerons, même si des situations plus générales peuvent être envisagées, sous multiples de ces quantités) et introduisons un maillage dont les noeuds sont les points de coordonnées (xmjn+ i'Axj',n'Atj') avec i' et n'a N. Le classique schéma centré aux différences finies explicité ci-dessous permet, partant des conditions initiales, le calcul de proche en proche des différentes valeurs prises par la fonction bj(x,t) aux différents noeuds du maillage (pour simplifier les notations nous

avons omis l'indice j associé au bruit corrélé considéré).

42tiI+: [1 bn;?r+ - I n.l]+ 2êA b i'2+L + bet bjié+l -bil + 1 [t7n;+Ll + n+ t _n 1+ 1 _ n]_( Dans la formule ci-dessus la quantité cj,fl représente l'évaluation de la quantité cE

(ca est un terme générique) au point de coordonnées (x1,a i'Axj',n 'Atj') .

Ici les pas de discrétisation Axj' et Atj' doivent être choisis suffisamment petits devant la longueur d'onde car nous voulons que le schéma numérique donne une approximation précise de la fonction solution de l'équation de transport (nous ne voulons

pas de dispersion).

La dernière étape dans la modélisation du bruit consiste à sélectionner les échantillons bi, qui appartiennent aux noeuds communs aux deux maillages et à les multiplier par g(iAx): -nous obtenons ainsi les différentes composantes du vecteur (de

l'espace des données) T1 (/3j).

Spécification d'une norme ou d'une semi-norme dans l'espace des données Un premier choix consisterait à prendre comme norme dans l'espace des données une norme discrète quelconque qui se voudrait une approximation (via une formule de quadrature) de la -norme définie en continu par: IIUIID dx dti u t) Nous pouvons aussi, ce qui constitue, comme nous le verrons au suivant, un meilleur choix, utiliser la semi-norme: Tn+ (4) o u est un vecteur quelconque de l'espace des données (i et n sont les indices représentant

les numéros d'échantillons respectivement en espace et en temps).

Par ailleurs d'autres choix sont possibles.

Spécification, dans chaque espace de fonctions génératrices du bruit, d'une norme pour laquelle chaque opérateur de modélisation du bruit établit une relation isométrique (ou l'approximation d'une relation isométrique) entre l'espace des fonctions génératrices du bruit et l'espace des données (muni de la semi-norme définie au 5) Nous précisons dans ce paragraphe la procédure pouvant être utilisée pour définir la

norme dans l'espace des fonctions génératrices associé à chaque bruit corrélé.

Un premier choix consisterait à munir chaque espace Bj de fonctions génératrices de bruit d'une norme discrète quelconque qui se voudrait une approximation (via une formule de quadrature) de la norme définie en continu par (là encore nous avons omis l'indice j pour alléger les notations): 11,31 = Xx gJz" (3X)), dé rt'2 (t)) auquel cas, en utilisant le fait que bj(x,T)=O pour xe [xni,Xmax], l'opérateur linéaire Tj réalise une isométrie entre Bj et D. Mais, comme le calcul de TS6 est effectué numériquement le caractère isométrique de l'opérateur Tj ne sera pas parfait mais seulement approché, l'approximation étant d'autant meilleure que les pas de discrétisationAx,At,Ax ',At 'sont petits. Le fait que nous n'ayons alors qu'approximativement une isométrie se traduira par des performances moins bonnes au

niveau de la mise en oeuvre algorithmique présentée au paragraphe 8.

Un meilleur choix consiste à munir chaque espace Bj de fonctions génératrices de bruit de la semi norme (là encore nous avons omis l'indice j pour alléger les notations): IIt|IIB = AN /\XL (E 2(jAX)\ (3 + IPIIB N E n+.) (5) Dans ce cas l'opérateur Tj réalise rigoureusement une isométrie entre Bj et D, du moins si la solution du schéma de discrétisation de l'équation de transport est

identiquement nulle pour n=N.

Définition de lafonction cot Nous définissons la fonction cot par la formule: C(mn,f1,fi2,v.fi1) = FNL(mo) + F(8m) + Y, TJj - d j=1 D

Notons que, dans la formule ci-dessus, FNL(mO) est un vecteur constant.

Recherche du modèle et des fonctions génératrices du bruit minimisant la fonction cot, cette recherche se faisant par une méthode algorithmique tirant parti, explicitement ou implicitement, des propriétés d'isométrie des opérateurs de modélisation du bruit

(cf. 6).

La méthode algorithmique exploite la spécificité du complément de Schur, spécificité associée au caractère isométrique des opérateurs Tj. Ceci peut se faire par exemple en réalisant que minimiser C(&m,,l3132,./39 revient à minimiser C'(&n) = C(Sm, l (ân), 1,2 (&n),....1, (Sm)) o (/(3k (e5n), i2 (mt),.../l3, (Sm)) minimise C(&m, 1, '32,-...,3) pour 8m donné. A noter que la hessienne associée à cette forme quadratique est le complément de Schur. La minimisation de C'(Sm) peut être réalisée par une quelconque méthode d'optimisation; la fonction cot C'(Sm) étant quadratique, l'utilisation d'une méthode de gradient conjugué est ici tout à

fait indiquée.

La détermination de (/3k (m), /2 (m),...fi (m)) est très aisée grâce au caractère isométrique des opérateurs Tj., notamment si nous avons effectués les choix (4) et (5) respectivement pour définir i1 ||D et || |lB.S'il n'y a qu'un seul type de bruit (J=1), ce qui est le cas pour notre illustration et si l'isométrie est parfaite (ce qui est également le cas grâce aux choix que nous avons effectués) la détermination de /,3 (m) demande simplement une évaluation du gradient de la forme quadratique. S'il y a superposition de plusieurs bruits (J>1), différents algorithmes peuvent être envisagés pour profiter du caractère isométrique

des opérateurs Ti. comme exposé pour le premier exemple d'application.

La Fig. 14B montrent les résultats obtenus (réponse sismique du modèle solution) par mise en oeuvre de la méthode dans le cas o la fonction g(x)est choisie sous forme très simple (fonction affine). On peut constater l'amélioration par rapport au résultat

d'inversion classique (Fig. 14A).

On a décrit des exemples de mise en oeuvre o le paramètre physique dont on modélise la distribution dans le sous-sol est l'impédance acoustique. Il est bien évident que la méthode dans sa définition la plus générale, peut s'appliquer à la recherche de la distribution de grandeurs physiques affectés par des bruits corrélés dans tout milieu hétérogène.

Claims

REVENDICATIONS

1) Méthode pour estimer, à partir de données obtenues par exploration d'une zone d'un milieu hétérogène, un modèle représentatif de la distribution, dans la zone, d'au moins une grandeur physique, affranchi de la présence de bruits corrélés pouvant entacher les données, caractérisée en ce qu'elle comporte les étapes suivantes: a) acquisition de mesures renseignant sur certaines caractéristiques physiques de la zone en suivant un protocole expérimental défini; b) spécification d'un opérateur de modélisation, qui, à un modèle de chaque grandeur physique associe des données synthétiques constituant la réponse du modèle, les mesures et les données synthétiques appartenant à l'espace de données; c) pour chaque bruit corrélé, repéré par un indice j variant de 1 à J, choix d'un opérateur de modélisation qui associe un bruit corrélé à une fonction génératrice du bruit appartenant à un espace (Bj) préalablement défini de fonctions génératrices du bruit; d) spécification d'une norme ou d'une semi-norme dans l'espace des données; e) spécification d'une semi-norme dans l'espace des fonctions génératrices de bruit, pour laquelle chaque opérateur de modélisation du bruit établit sensiblement une relation isométrique entre l'espace des fonctions génératrices du bruit et l'espace des données; f) définition d'une fonction cot quantifiant l'écart entre les mesures d'une part et la superposition de la réponse du modèle et des bruits corrélés associés aux fonctions génératrices de bruit d'autre part; et g) ajustement du modèle et des fonctions génératrices du bruit en minimisant la fonction cot, par une méthode algorithmique tirant partie, des propriétés d'isométrie des

opérateurs de modélisation du bruit.

2) Méthode selon la revendication 1, dans laquelle on recherche la distribution en fonction de la profondeur de l'impédance acoustique dans le milieu, les bruits corrélés affectant les données sont des ondes de tube identifiées chacune par des paramètres caractérisant leur propagation, les données mesurées sont des données de PSV en utilisant des capteurs adaptés à détecter le déplacement des particules dans le milieu en réponse à une excitation sismique localisée, l'emplacement des capteurs, la durée d'enregistrement et les points d'échantillonnage en temps étant définis, et l'opérateur de modélisation choisi associe les données synthétiques à une distribution d'impédance acoustique en fonction de la profondeur évaluée en temps de parcours et à la contrainte verticale mesurée en fonction du temps à la profondeur du premier capteur.

3) Méthode selon la revendication 2, dans laquelle la fonction cot quantifiant

l'écart est le carré de la semi-norme de cet écart dans l'espace des données.

4) Méthode selon la revendication 2 ou 3, dans laquelle l'ajustement du modèle et des fonctions génératrices du bruit est obtenu par une méthode de relaxation par bloc pour éliminer les inconnues correspondant à chaque fonction génératrice de bruits corrélés, cette méthode de relaxation étant mise en oeuvre à l'intérieur des itérations d'un algorithme de

quasi-Newton pour le calcul du modèle.

) Méthode selon l'une des revendications 2 à 4, dans laquelle le calcul numérique

de l'image d'un modèle par l'opérateur de modélisation est effectué par résolution numérique de l'équation des ondes ID pour le modèle considéré, en sélectionnant des valeurs prises par le déplacement des particules aux emplacements de capteurs et aux points d'échantillonnage en temps préalablement spécifiés et en appliquant un opérateur

propre à compenser les effets de divergence sphérique et d'atténuation.

6) Méthode selon l'une des revendications 2 à 5, dans laquelle l'opérateur

numérique de modélisation du bruit est un schéma numérique centré aux différences finies pour discrétiser l'équation de transport du bruit et la fonction génératrice de bruit intervenant comme condition initiale au bord de la zone d'observation appartient à, un

espace (Bj) de fonctions du temps à support dans un intervalle de temps donné.

7) Méthode selon l'une des revendications 2 à 6, dans laquelle dans laquelle la

semi-norme choisie pour l'espace des données est:

0 I- O N1 ( E 2

117s11 A t t ft, E n+1(s + uj2) 2=0 n=O et la semi-norme choisie pour l'espace des fonctions génératrices de bruit est: I1,31IB (xAxsARI EAt 1 ( n.+ +S n)2) 8) Méthode selon la revendication 1, dans laquelle on recherche la distribution de perturbations, par rapport à un modèle de référence préalablement choisi, de l'impédance et de la vitesse dans la dite zone du milieu, les bruits corrélés affectant les données sont dus à des réflexions multiples dont la cinématique et les variations d'amplitude avec l'offset ont été préalablement estimés, les données mesurées sont captées par des capteurs sismiques de surface, l'emplacement des dits capteurs, le mode d'excitation sismique, la durée d'enregistrement et le pas d'échantillonnage en temps étant définis, et l'opérateur de i0 modélisation étant défini via une linéarisation de l'équation des ondes autour du modèle de référence. 9) Méthode selon la revendication 8, dans laquelle la fonction cot quantifiant

l'écart est le carré de la semni-norme de cet écart dans l'espace des données.

) Méthode selon la revendication 8 ou 9, dans laquelle l'ajustement du modèle et des fonctions génératrices du bruit est obtenu par une méthode de relaxation par bloc pour éliminer les inconnues correspondant à chaque fonction génératrice de bruits corrélés, cette méthode de relaxation étant mise en oeuvre à l'intérieur des itérations d'un algorithme de

gradient conjugué pour le calcul du modèle.