Brevet FR2691812A1 — publié le 19931203

DOMAINE TECHNIQUE DE L'INVENTION La présente invention concerne le domaine de la prospection géophysique ou sismique. ETAT GENERAL DE LA TECHNIQUE On sait que la prospection géophysique ou sismique, destinée à permettre l'étude de la structure interne de la terre, consiste fréquemment à enregistrer des ondes sismiques issues d'une source artificielle (explosif, générateur de vibrations, etc) après leur réflexion et/ou leur réfraction dans le sous-sol. Cette méthode de prospection est fondée sur la propriété qu'ont les ondes sonores provoquées par une source artificielle de subir des réfractions et des réflexions aux surfaces de contact de couches ayant des vitesses de transmission et des densités différentes, suivant des lois analogues à celles de l'optique. Les systèmes utilisés pour la prospection géophysique comprennent donc généralement: une source artificielle d'une onde sonore ou ébranlement à base d'explosif ou équivalent, un réseau de capteurs, tels que des géophones ou équivalents, disposés sur la surface du sol et sensibles aux ondes réfléchies vers le haut par les couches internes du sol, un enregistreur relié à ces capteurs pour enregistrer les signaux électriques issus de ceux-ci, sous forme de sismogrammes, des moyens de traitement des signaux enregistrés et des moyens de visualisation d'un résultat exploitable après traitement de ces signaux. Les moyens de traitement doivent effectuer un traitement assez complexe des signaux enregistrés, car ces signaux ne peuvent être utilisés directement pour visualiser les caractéristiques des couches analysées du sous-sol, du fait qu'ils sont pertubés. Les perturbations de ces signaux sont assez bien connues de l'homme de l'art Elles résultent en particulier du fait que chaque point lieu d'une variation d'impédance acoustique renvoie les ondes sonores incidentes dans une pluralité de directions, ce qui conduit donc non pas à un point unique, mais à une trace hyperbolique dans l'enregistrement, puisque chaque point diffractant est situé à des distances différentes des divers capteurs utilisés. Pour éliminer cette perturbation, on utilise généralement une technique dite de migration qui focalise l'énergie diffractée et permet de donner une image claire du sous-sol, considéré comme un ensemble de sources secondaires De nombreuses techniques de migration sont actuellement utilisées Parmi les diverses techniques ainsi connues, la présente invention concerne plus précisément la catégorie des méthodes récursives basées sur l'extrapolation en profondeur (z) du champ d'onde enregistré en surface (z=O) Dans cette catégorie, les méthodes les plus couramment utilisées travaillent dans le domaine fréquence temporelle () et espace soit (x,z) dans le cas à deux dimensions ( 2 D), et (x,y, z) dans le cas à trois dimensions ( 3 D) On distingue encore les méthodes d'extrapolation implicites pour lesquelles le champ extrapolé s'obtient par résolution d'un système linéaire et les méthodes explicites pour lesquelles le champ extrapolé s'obtient par application d'un opérateur de convolution spatiale Fo(w,c,x,y) suivant l'équation: Wo,z+dz(x Y) = Fo(C Oecxy) * W w,z(Xy) ( 1) o le signe * dénote la convolution en (x,y). W,z(xy) s'obtient à partir de Wt,z(x,y), qui est le champ d'onde à la profondeur z, par une transformation de Fourier sur la variable t L'opérateur d'extrapolation FO(cec,x,y) dépend de c qui est la vitesse de propagation locale des ondes dans le milieu. Cette équation permet de façon récursive de calculer le champ d'onde W,, z(x,y), pour toutes les profondeurs z, à partir d'enregistrements de surface pour la profondeur z=O fournissant W,,z=o(x,y). Les méthodes implicites ont l'avantage d'être inconditionellement stables, mais elles conduisent à des imprécisions numériques croissantes avec les pendages des réflecteurs géologiques Le "pendage" correspond à l'inclinaison des surfaces sensibles des capteurs par rapport à l'horizontale, soit à l'inclinaison de la propagation par rapport à la verticale De plus les méthodes implicites se généralisent difficilement au cas 3 D pour lequel on est conduit à utiliser une approximation anisotrope dite de "splitting" qui consiste à décomposer l'extrapolation en 2 étapes 2 D suivant les directions x et y Il en résulte une mauvaise image des réflecteurs dont le pendage n'est ni dans la direction x ni dans la direction y. Les méthodes explicites ont l'avantage de se généraliser facilement au cas 3 D En revanche, le calcul des coefficients des opérateurs de convolution est délicat pour garantir stabilité et précision Dans leur mise en oeuvre en général, les opérateurs sont calculés à l'avance dans une table en fonction d'un échantillonnage des paramètres de vitesse du sous-sol. NOTATIONS UTILISEES Les données dépendant de (x,y) ne sont connues dans la pratique que sur une grille discrète de valeurs, correspondant en surface aux positions des capteurs Les pas d'échantillonnages spatiaux en x et y sont Dx et Dy. (wk Xky) = (fréquence, nombre d'onde en x, nombre d'onde en y) sont les variables correspondant à t,x,y par transformation de Fourier knyqx*=/Dx dénote le nombre d'onde de Nyquist pour le pas Dx, et ky Ty=n/Dy de même pour le pas Dy Nous désignerons désormais le champ d'onde par W 0,z(mx,my), mx étant l'indice sur l'axe des x et my l'indice sur l'axe des y, Ww,z(mx,my) désigne le champ échantillonné en (x,y) sur une grille de maillage Dx en coordonnée x et Dy en coordonnée y Les indices (mx,my) correspondent donc à des coordonnées (mx Dx,my Dy) Le champ des vitesses de propagation à la profondeur z est noté cz(mx,my) L'opérateur d'extrapolation Fo(w,c,x,y), d'un champ d'ondes montantes, d'une profondeur Dz pour une vitesse de propagation c et une fréquence w, s'écrit dans le domaine (kx,ky): Fo(oe,c,kx,ky) = expj Dz(w 2/c 2-kx 2-ky 2)1/2 ( 2) L'équation ( 1) s'écrivant dans le domaine (kx,ky): WW,z+dz(kx,ky)=Fo(o,c,kx,ky)WW,Z(kx,ky) ( 3) Pour appliquer la méthode à des ondes se propageant vers le bas, il suffit de prendre le conjugué de Fo. Les méthodes explicites d'extrapolation en profondeur sont basées sur le fait que F O ne dépend pas de et de c séparément mais du rapport w/c seulement On peut donc écrire Fo(w,c,kx,ky) = Fo(w/c,kx,ky). ETAT DE LA TECHNIOUE DANS LE DOMAINE DE L'EXTRAPOLATION 2 D Dans le cas de l'extrapolation à deux dimensions, l'opérateur d'extrapolation s'écrit: Fo(co/c,k) = exp j Dz(o 02/c 2 k 2)1/2 Les méthodes explicites consistent à faire la synthèse spectrale de F 0, c'est-à-dire calculer une banque de coefficients a(w/c,n) pour w/c parcourant un échantillonnage de 1 ' intervalle l in/cmax,oax/cminlr, o lmin,%maxl représente le contenu fréquentiel des données et lCmin, Cmaxl les bornes des vitesses de propagation, et n El 0,Nl, tels que l'opérateur synthétisé: n=N F(o/c,k)= j a(w/c,n) cos(n Dx k) n= O approxime Fo(w/c,k) pour k El 0,kc=wsin O/cl, bande passante, et soit de module inférieur à 1 pour k Elkc,knyq=T/DXl, bande d'arrêt O désigne le pendage géologique maximum correspondant à une propagation précise kc=wsin O/c est le nombre d'onde correspondant au w/c donné et à un angle de propagation O par rapport à la verticale La contrainte sur le module est rendue nécessaire par le fait que ces opérateurs étant destinés à être appliqués récursivement, un module supérieur à 1 générerait une instabilité. Il faut donc assurer la précision sur la bande passante et la stabilité sur la bande d'arrêt. Holberg l 1 l utilise pour celà un procédé classique de moindres carrés non linéaire Cette méthode entraine cependant un coût de calcul élevé. Hale l 2 l propose un procédé plus rapide en calculant les coefficients a(n) qui égalent les premiers coefficients des développements en série de Taylor de F(k) et Fo(k) en k La stabilité est assurée en imposant à F(k) d'être nul pour certaines valeurs de k dans la bande d'arrêt Cependant, une série de Taylor tronquée n'étant une bonne approximation que pour des faibles valeurs de k, donc pour de faibles pendages O, cette méthode n'est pas optimale. ETAT DE LA TECHNIQUE DANS LE DOMAINE DE L'EXTRAPOLATION 3 D Dans le cas à 3 dimensions, l'opérateur d'extrapolation d'un champ d'ondes montantes, d'une profondeur Dz pour une vitesse de propagation c et une fréquence w, est: Fo(w/c,kx,ky) = exp j Dz(w 2/c 2 kx 2 ky 2)1/2 Dans la référence l 3 l Blaquière propose de réaliser l'extrapolation explicite 3 D en généralisant le cas 2 D, c'est-à-dire en faisant une synthèse spectrale 2 D: une banque de coefficients a(w/c, nx,ny) est calculée, pour w/c parcourant un échantillonnage de l'intervalle l%in/Cmax,max/cminl, o lwmin,maxl représente le contenu fréquentiel des données et lcmin,cmaxl les bornes des vitesses de propagation, et nx El-Nx,Nxl, ny El-Ny,Nyl, tels que F(w/c,kx,ky) définit ci-dessous: nx= Nx ny= Ny F(w/C,kx,ky)= E 2 a(w/c,nx, ny)explj(nx Dx kx+ny Dy ky)l n=-Nx ny=-Ny approxime Fo(w/c,kx, ky) sur le disque (kx 2 +ky 2)l/2 s w sin 0/c, et soit de module inférieur à 1 ailleurs O est le pendage géologique maximum qui correspond à une extrapolation précise Ce paramètre est choisit par l'utilisateur. Les a(w/c,nx,ny) sont appliqués de la manière suivante W, tz+dz (Ilim) n,= Nx ny= Ny z E aloe/cz+dz(m X,my),n X, nyl Ww,z(mx-nx,my-ny) nx-Nx ny=-Ny qui est la version en (mxmy) de l'équation: Wco,z+dz(kx,ky) = F(w/c,kx,ky) W ,z(kx,ky) Hale, dans la référence l 4 l, reprend l'idée de la tabulation précédente en tirant partie de surcroit de la symétrie circulaire en (kx,ky) de Fo(o/C,kx, ky) Il utilise pour celà les coefficients a(w/c,n) à une dimension correspondant à une extrapolation 2-D et le processus de transformation proposé par McClellan dans la référence l 5 l Hale écrit avec Dx=Dy: F O (o/c,kx,ky)=expj Dz(o 2/c 2-k 2)l/2 avec k = (kx 2 +ky 2)1/2 et fait la synthèse spectrale: n=N F(a/c,kx,ky)= j a(o/c,n) cos(n Dx k) n= O O Les a(o/c,n) sont calculés de telle sorte que F(w/c,k) approxime Fo(o/c,k) pour les petites valeurs de k et soit partout de module inférieur à 1. Les cos(n Dx k) peuvent se déduire de cos(Dx k) par la formule récursive des polynômes de Chebychev: cos(nh) = 2 cos(h) cosl(n-1)hl cosl(n-2)hl o h=k Dx = longueur d'onde normalisée. De plus cos(Dx k) peut être approximé par un filtre à deux dimensions pour les pas d'échantillonnage Dx et Dy: nx= N Mx ny= NMY cos(Dx k) = 2 2 M(nx,ny) exp j(nx Dx kx + ny Dy ky) nx=-N Mx ny=-N My Le développement le plus simple de cos(Dx k) est basé sur la transformée de McClellan à 9 termes, correspondant à N Mx=Nmy=l: cos(Dx k) = -1 + 0,5 l 1 + cos(Dx kx)l l 1 + cos(Dx ky)l Hale propose également d'exploiter une transformée de McClellan améliorée à 17 termes, correspondant à NMX=N My= 2 avec certains termes nuls Mais, comme on le verra par la suite, même avec cette transformée améliorée, le résultat final comporte des pertubations. Plus précisément, la méthode consiste donc à calculer les a(w/c,n), w/c parcourant un échantillonnage de l'intervalle loin/cmaxfmax/cmnl et les M(nx,ny). Définissons l'opération W' = M*W par: nx=N Mx ny=N My W,,z+dz(mx,my) = E M(nx,ny) Wo,z(mx-nx,my-ny) nx=-Nmx ny=-Nmy Avec cette notation, le calcul du champ extrapolé Ww,z+dz(mx,my) à partir de Ww,z(mx,my) se fait, pour chaque w, par le pseudo-code suivant: 1) Initialisation: To =Wz T 1 = M*Wz Wz+dz(mx,my)=bl(e/cz+dz(m X,my), OlTO(mx,my)+ b lE/cz+dz(mx, my) 1 l Tl(mx, my) 2) Itération, pour n= 2 jusqu'à N faire: T = 2 * M * T 1 To WZ+d Z(m X,m Y) = WZ+d Z(m Xmy) + blw/c Z+d Z(m Xmy),nl T(mx,my) To = T 1 T 1 = T 3) Fin de la boucle en n. BUT DE L'INVENTION La présente invention a pour but de perfectionner les techniques existantes de traitement de signaux pour prospection géophysique exploitant un opérateur d'extrapolation d'un champ d'onde. L'objet de l'invention est en particulier de faciliter le calcul des coefficients des opérateurs de convolution et de simplifier la construction et l'application des opérateurs en 3 D. APPLICATIONS INDUSTRIELLES DE L'INVENTION L'extrapolation du champ d'onde est nécessaire pour réaliser la migration, qui est une étape importante du traitement des données géophysiques La présente invention peut être appliquée à la migration à 2 ou 3 dimensions, avant ou après sommation, à la migration en profondeur ainsi qu'à la migration en temps Elle peut s'appliquer aussi à la modélisation, qui est l'opération inverse de la migration. EXPOSE GENERAL DE L'INVENTION Pour atteindre le but précité il est proposé dans le cadre de la présente invention, un procédé d'extrapolation d'un ensemble de traces enregistrées à l'aide d'un réseau de capteurs, dans un processus de prospection géophysique, caractérisé par le fait qu'il comprend les étapes consistant à: i) approximer un moins Laplacien Lo(k,ky) = k 2 +ky 2 par la somme de deux filtres à une dimension, et ii) approximer un opérateur d'extrapolation Fo par un polynome du moins Laplacien L. Selon une caractéristique avantageuse de la présente invention, l'étape i) consiste à approximer le moins Laplacien par la somme de deux filtres à une dimension de la forme: n=NLX n= N Ly L = j dx(n) cos(n Dx k) + J dy(n) cos(n Dy ky) n= O n=O Selon une autre caractéristique avantageuse de la présente invention, à l'étape ii) l'opérateur d'extrapolation est approximé par un polynome du moins Laplacien de la forme: n=N F(oe/c,k,, ky)= 1 p(w/c,n) LOI n=O Dans le cadre de la présente invention il est également proposé un système pour prospection géophysique du type comprenant: une source artificielle d'une onde sonore, un réseau de capteurs disposés sur la surface du sol et sensibles aux ondes réfléchies vers le haut par les couches internes du sol, un enregistreur relié à ces capteurs pour enregistrer les signaux électriques issus de ceux-ci, des moyens de traitement des signaux enregistrés et des moyens de visualisation d'un résultat exploitable après traitement de ces signaux, caractérisé par le fait que: les moyens de traitement sont adaptés pour décomposer un opérateur d'extrapolation en partie réelle et imagnaire, puis effectuer une synthèse optimale au sens de la norme L-w sur chacune des deux parties. D'autres buts, caractéristiques et avantages de la présente invention apparaîtront à la lecture de description détaillée qui va suivre et en regard des il dessins annexés donnés à titre d'exemple non limitatif et sur lesquels: la figure 1 représente une vue schématique, sous forme de blocs fonctionnels, d'un système conforme à la présente invention, la figure 2 a représente la partie réelle d'un filtre d'extrapolation conforme à la présente invention, la figure 2 b représente le module maximal d'un tel filtre, la figure 2 c représente la partie imaginaire du filtre d'extrapolation, la figure 2 d représente le module effectif du filtre d'extrapolation conforme à la présente invention, la figure 3 représente la réponse impulsionnelle obtenue en 2 D dans le cadre de la présente invention, la figure 4 représente une réponse impulsionnelle 3-D avec utilisation d'un Laplacien à 17 termes conforme à la présente invention, pour une tranche verticale y=O, la figure 5 représente une réponse impulsionnelle 3-D avec utilisation d'un Laplacien à 17 termes conforme à la présente invention, pour une tranche horizontale z= 450 m, la figure 6 représente une réponse impulsionnelle 3-D avec transformée de McClellan à 17 termes conforme à l'état de la technique, pour une tranche horizontale z= 450 m. STRUCTURE GENERALE DE L'INVENTION La structure générale du système utilisé dans la présente invention est sensiblement conforme à la structure connue rappelée dans le préambule. Ce système comprend une source artificielle 10 d'une onde sonore ou ébranlement à base d'explosif ou équivalent, un réseau de capteurs 20, tels que des géophones ou équivalents, disposés sur la surface S du sol et sensibles aux ondes réfléchies vers le haut par les couches internes du sol, un enregistreur 30 relié à ces capteurs 20 pour enregistrer les signaux électriques issus de ceux-ci, sous forme de sismogrammes, des moyens de traitement 40 des signaux enregistrés et des moyens 50 de visualisation d'un résultat exploitable après traitement de ces signaux. La présente invention se distingue cependant de l'état de la technique par la structure et la fonction des moyens de traitement 40. On va par la suite décrire la fonction de ces moyens de traitement 40 dans le cadre de la présente invention. INVENTION DANS LE DOMAINE DE L'EXTRAPOLATION 2 D Comme indiqué précédemment selon l'invention, on propose d'effectuer la synthèse spectrale en décomposant Fo(k) en partie réelle et imaginaire après une rotation de phase, puis d'effectuer une synthèse optimale au sens de la norme L-co sur chacune des deux parties. Ces 2 synthèses optimales L-0 x peuvent être effectuées par exemple en utilisant le "Remez exchange algorithm" ou "algorithme de Remez", qui est décrit dans la référence l 6 l, et qui résoud le problème suivant: Soit So(h) définit pour h El 0,nl, un spectre réel symétrique donné. Soit W(h) définit pour h El 0,Trl, une fonction poids réelle positive donné. Soit N la demi-longueur du filtre donnée. n=N S(h)= X a(n) cos(nh) est le spectre synthétisé. n= O E(h)=W(h)lS(h)-So(h)l est la fonction d'erreur pondérée entre le spectre synthétisé et le spectre idéal donné. Trouver les a(n), n ElO,Nl tels que: Norme L-m de E(h) = maxh Elo 0,,rl IE(h) l soit minimum. Plus précisément, la procédure proposée par la présente inventon comprend les étapes qui consistent à: i) Décomposer Fo(k) exp(jc) en parties réelles et imaginaires: F O (k) exp(j 4 c) = R O o(k) + j Io(k) : est calculé de telle façon que R O (k) présente un maximum pour k=kc En l'occurence, Ro(k,)= 1 et Io(kc)= 0. ii) Par l'algorithme de Remez ou tout autre algorithme, calculer ar(n), n El 0,Nl, tels que: n=N R(k) = 2 ar(n) cos(n Dx k) n= O corresponde au minimum de maxkelO,kcl IR(k)-RO(k)l, tout en imposant aux extrema de R(k) sur lkc,knyql d'avoir des ordonnées a, o a est un nombre plus petit que 1 La figure 2 a annexée représente la partie réelle d'un tel filtre d'extrapolation. iii) Puis calculer la fonction de pondération W(k) pour k Elkc,knyql: W(k)=lM 2 (k)-R 2 (k) l-1/2 o M(k) est le module maximum dans la bande lkc,knyql, définit par l'utilisateur Le minimum de M(k) doit être supérieur à " et son maximum inférieur à 1 un exemple d'une telle fonction est montré sur la figure 2 b annexée. iv) Par l'algorithme de Remez ou tout autre algorithme, calculer ai(n), n El 0,Nl, tels que: n=N I(k) = 2 ai(n) cos(n Dx k) n= O corresponde au minimum de maxkel O o,kcl I(k)-IO O (k)l, tout en imposant que maxk Elk C,knyql W(k)II(k)l soit égal à 1 La figure 2 c annexée représente la partie imaginaire d'un tel filtre d'extrapolation. F(k)=lR(k)+j I(k)l exp(-j>c) répond alors au problème posé Les coefficients recherchés sont: a(n) =lar(n)+jai(n)l exp(-j 4 c) On a représenté sur la figure 2 d annexée le module du filtre d'extrapolation. La méthode de synthèse spectrale qui vient d'être décrite présente des avantages par rapport aux deux méthodes existantes. L'une (référence l 1 l) est relativement optimale mais le calcul des coefficients est long, et l'optimalité l'est au sens de la norme L-2 Or, il est préférable que l'approximation en bande passante soit faite pour la norme L-^, comme dans la méthode proposée En effet ces opérateurs doivent être utilisés de façon récursive. Minimiser l'erreur sur l'opérateur global FP(k), o p est le nombre de pas d'extrapolation en z, pour n'importe quelle norme, revient en effet à minimiser l'erreur sur l'opérateur unitaire pour la norme L-. L'autre (référence l 2 l) est relativement rapide mais sous-optimale Elle revient à minimiser l'erreur sur un intervalle local kl 0,El, avec k infiniment petit, alors que la méthode proposée minimise l'erreur sur tout l'intervalle k El 0,k,== sin S/cl, O étant l'angle maximum de propagation exacte choisi par l'utilisateur, c'est-à-dire que l'erreur est minimisée pour tous les pendages inférieurs à un seuil fixé par l'utilisateur. On a représenté sur la figure 3 la réponse impulsionnelle obtenue pour une plage de fréquence comprise entre fmin = 2 Hz et fmax = 48 Hz, une vitesse c = 2000 m/s et une cible à une profondeur z = 900 m. INVENTION DANS LE DOMAINE DE L'EXTRAPOLATION 3 D Comme indiqué précédemment, l'un des buts importants de la présente invention est de simplifier le processus de calcul mis en oeuvre pour extrapoler le champ d'onde 3-D représenté par les traces enregistrées, et en améliorer la précision. Ce but est atteint selon la présente invention grâce à un procédé d'extrapolation en profondeur d'un champ d'onde enregistré à l'aide d'un réseau de capteurs dans un processus de prospection géophysique comprenant les étapes qui consistent à mettre l'opérateur d'extrapolation sous la forme: Fo(w/c,kx,ky) = exp j Dz(m 2/c 2-Lo)1/2 dans lequel L O désigne le moins Laplacien: L O O(kx,ky) = kx 2 +ky 2. Plus précisément selon l'invention le procédé d'extrapolation comprend les étapes qui consistent à: i) approximer un moins Laplacien L O (kx, ky) = kx 2 +ky 2 par la somme de deux filtres à une dimension, et ii) approximer un opérateur d'extrapolation Fo par un polynome du moins Laplacien Lo. Le moins Laplacien est de préférence approximé par la somme de deux filtres à une dimension de la forme: n=N Lx n= N Ly L = dx(n) cos(n Dx kx) + X dy(n) cos(n Dy ky) n= O n=O Par ailleurs à l'étape ii) l'opérateur d'extrapolation est de préférence approximé par un polynome du moins Laplacien de la forme: n=N F(o/c,kx,ky)= p(o/c,n) L On' n= O L'opérateur moins Laplacien L O peut être scindé, sans approximation, sous forme de la somme de deux opérateurs de moins dérivée seconde, chacun d'eux pouvant être approximé par un filtre à une dimension avec N Lx et NLY termes en cosinus. La présente invention exploite la symétrie circulaire de l'opérateur d'extrapolation, comme le propose Hale dans le cadre d'une mise en ouvre de transformée de McClellan. Cependant la présente invention permet d'éviter le calcul du filtre 2-D approximant coslDx(kx 2 +ky 2)1/2 l grâce à l'utilisation d'un opérateur moins Laplacien, L O = kx 2 + ky 2, à l'aide duquel l'opérateur d'extrapolation F O est développé sous la forme de polynome et dont la synthèse ne nécessite que des filtres 1-D. Plus précisément, la procédure proposée dans le cadre de la présente invention comprend les étapes qui consistent à: 1) Faire une synthèse spectrale de l'opérateur moins dérivée seconde en x, c'est-à- dire calculer un filtre symétrique de longueur 2 N Lx+l, pour le pas d'échantillonnage Dx, dont le spectre Dx(k), définit par: n= N Lx Dx(k) = dx(n)exp(j n Dx k) avec dx(-n)=dx(n) n=-N Lx approxime le spectre de la moins dérivée seconde Do(k)=k 2, sur k El 0,axknyqxl, avec O certaine valeur). 2) Faire la même chose pour y: n= N Ly Dy(k) = Z dy(n) exp(j n Dy k) avec dy(-n)=dy(n) n=-N Ly 3) Définir L(kx,ky) = Dx(kx) + Dy(ky) qui est la synthèse spectrale du moins Laplacien Lo(kx,ky)=kx 2 +ky 2 sur le rectangle k 5 l -oxknyx, xknycxl, ky El kny, k nyyl l. En pratique, on calcule d'abord les dx(n), dy(n) (qui sont a priori indépendants de w/c mais qui peuvent aussi en dépendre) Ensuite les bornes lLmin Lmaxl de L(kx,ky) sont calculées en tant que somme des bornes de Dx(k) et de Dy(k) On choisit O qui est l'angle de propagation maximum considéré par rapport à la verticale. Pour une pulsation et une vitesse donnée, celà correspond à un nombre d'onde kc, sin O/c Puis pour chaque w/c parcourant un échantillonnage de l'intervalle l%min/Cmaxpax/cminl, on calcule les p( /c,n), en tant que développement polynomial de F O (w/c,L) F(w/c,L) doit approximer Fo(w/c,L) pour LElLmin,Lc=kc 2 =W 2 sin 28/c 2 l et doit être de module inférieur à 1 pour LElLc,t Lmxl. Définissons l'opération W'=L*W par: n= N Lx n= NL W' (nx, ny) = dx(n) W(nx-n,ny) + S dy(n) W(nx,ny-n) n=-N Lx n=-N Ly W' approxime le moins Laplacien de W. Avec cette notation, le calcul du champ extrapolé WW,z+dz(mx,my) à partir de Wo,z(mx,my) se fait, pour chaque a, par le pseudo-code comprenant les étapes qui consistent à: 1) Initialiser: Wz+dz(mxmy) = plC/c(mx,my),Nl Wz(mx,my) puis 2) Pour n=N-1 jusqu'à O faire: Wz+dz(mx,my) = (L*Wz+ dz)(mx,my) + plw/c(Im C,my),nl Wz(mx,my) On va maintenant préciser la méthode proposée dans le cadre de la présente invention pour calculer les coefficients en extrapolation 3 D. Nous avons vu comment ramener la synthèse de Fo(w/c,kx,ky) à celle de coefficients 1-D: une synthèse polynomiale 1-D, p(w/c,n), et deux synthèses spectrales 1-D de spectres réels symétriques, dx(n) et dy(n) Nous allons décrire une procédure de calcul pour ces coefficients. Les synthèses de dx(n) et dy(n) peuvent se faire simplement en utilisant l'algorithme de Remez, avec (h)=h 2, W(h)= 1 pour h El 0,axnl,W(h)= pour h Elaxn,nl, o s est un petit nombre de l'ordre de 0,05. On montre que la synthèse polynomiale peut se ramener à une synthèse spectrale De plus celle-ci peut être opérée avec le procédé décrit précédemment pour l'extrapolation 2 D. La synthèse polynomiale à effectuer est la suivante: Soit une fonction complexe Fo(L) donnée, ainsi qu'un intervalle lLmin Lmaxl et une valeur Lc Il faut trouver des coefficients p(n), n El 0,Nl, tels que: n=N F(L)= p(n) Ln soit n= O une approximation de Fo(L) sur lLmin,Lcl et que le module de F(L) soit inférieur à 1 sur lLó,Llxl Faisons le changement de variable: L= 0,5 (Lmin-Lmax) c OS (h) + 0,5 (Lmin+Lmax)=g(h) avec LElL ^n,Lmxl, h El 0,nl. F O lg(h)l = So(h) est un spectre symétrique connu. Posons h O = g-l(Lc). On résoud le problème de synthèse spectrale suivant, en utilisant par exemple la méthode proposée précédemment: Trouver les t(n), n El 0,Nl, tels que n=N S(h) = j t(n) cos(nh) n=O approxime So(h) pour h ElO,h,l et soit de module inférieur à 1 sur lh,1 rl. Les t(n) définissants S(h) étant calculés, la méthode décrite dans la suite, F(L)=S(g-l(L)) répond au problème posé. En effet, le changement de variable L=g(h) transforme un polynôme de degré N en L en un polynôme de degré N en cos(h) et donc, d'après les polynômes de Chebychev, en une somme de cos(nh) pour N variant de O à N Donc réciproquement, h=g-1 (L), transforme un spectre de cosinus en un polynôme De plus, la norme L- est invariante par changement de variable Donc si la synthèse spectrale est optimale pour la norme L-o, la synthèse polynomiale le sera aussi. Plus précisément, on peut écrire F(L) sous la forme: n=N F(L) = 2 t(n) Tnl( 2 L-Lmin-Lmax)/(Ljin-Lmax)l ( 4) n= O T(x) étant le polynôme de Chebychev de degré n définit par Tn(cos h) = cos nh. On se ramène ensuite à la forme: n=N F(L)= E p(n)(AL/I Lyx)n ( 5) n= O o A est un facteur d'échelle qui vaut typiquement 4 Le fait de développer en (AL/Lmax)n plutôt qu'en Li permet d'avoir des coefficients p(n) de moindre dynamique et donc d'avoir moins de problèmes numériques. Le passage de la forme ( 4) à la forme ( 5) pour F(L) se fait en écrivant: n=m p(m) = 5 a(m,n) t(n) n= O o les a(m,n) sont définis par: n=m Tml ( 2 L-Lmtn-Lmx)/(Lmin-Lmax)l = a(m,n) (AL/T Lax)a n=O Les a(m,n) se calculent récursivement en écrivant: Tm(x) = 2 x Tm-l(X) Tm-2 (X) avec x= ( 2 L-Lmin-mx) / (lin-lmax) Lorsque N est grand, il peut être impossible de trouver un A qui permette de ne pas avoir de problèmes de précision machine Dans ce cas, il faut en rester à la forme en polynômes de Chebychev ( 4) Le pseudo- code de la méthode devient alors similaire à celui de la méthode de Hale, mais l'opérateur M, qui a une forme rectangulaire dans le domaine (x,y) est remplacé par l'opérateur ( 2 L-Lmin- Lmax) / (Lmin-Lmax) qui a une forme de croix dans le domaine (x,y). Avec la forme ( 4) de F(L), le pseudo-code utilisé est le suivant - 1) initialiser: To=Wz Tl=L*Wz Wz+dz(mx,my) = tlO/cz+dz(m X,my),Ol TO(mx,my) + tlCO/cz+dz(m X,my),1 l Tl(mx,my) 2) Pour n= 2 jusqu'à N faire: T= 2 *L*T 1 l-TO Wz+dz(m X,my) = Wz+dz(mx,my) + tlC/Cz+dz(mx,my), nl T(mx,my) T O =T 1 T 1 =T Fin de la boucle en n avec W'=L*W définit par n= N Lx W'(nxny) = l 2 dx(n)/(Lmin-Lmaxl W(nx-nny) n=Nny nf N Ly + l 2 dy(n)/(Lmin-Lmaxl W(nx, ny-n) n 7 =-N Ly (Lmin+Ix)W(nx, ny)/ (I Rn-I*ax). Avec la forme ( 5) de F(L), le pseudo-code utilisé est le suivant: 1) initialiser Wz+dz(mx,my) = p EO/c(mx,my),Nl Wz(mx,my) puis 2) Pour n=N- 1 jusqu'à O faire: Wz+dz(mx,my) = (L*Wz+dz)(mx,my) + plo/c(n%,m%),nl Wz(mx,my) avec W'=L*W définit par n N Lx W' (nx,ny) = lAdx(n)/Lmaxl w(nx-n,ny) n=-N Lx n= N Ly + lAdy(n)/Lmaxl W(nx,ny-n) n=-N Ly. EXTENSION DE LA METHODE A L'EXTRAPOLATION 3 D EN TEMPS L'extrapolation en temps d'un champ d'onde permet de faire la migration en temps, avant ou après sommation, une étape souvent utilisée en traitement géophysique La migration en temps est moins précise que la migration en profondeur, mais elle peut se contenter d'un modèle de vitesse cz(mx,my) beaucoup moins précis, ce qui en fait l'utilité. Contrairement aux méthodes de l'état de la technique précédemment décrites, la méthode proposée par la présente invention peut s'étendre au cas de l'extrapolation en temps La fonction d'extrapolation en kxky est dans ce cas, pour un pas d'extrapolation temporel Dt: G O (,,c, kx,ky)≥expj Dt(C 2-c 2 L O O)1/2 avec Lo (kx, ky)=kx 2 +ky 2 Go(t,c) dépend cette fois-ci séparément de W et de c et ne peut se mettre sous la forme G O (/c) Une méthode basée sur la tabulation nécessiterait un filtre pour chaque et chaque c, ce qui serait prohibitif en temps de calcul et en stockage. La méthode proposée, utilisant un développement en monômes de la forme Ln, permet de ne tabuler les coefficients que pour la vitesse minimale cmizn, et de déduire les coefficients correspondant à une vitesse c quelconque par une opération très simple qui peut être faite au moment de l'application. On écrit: Go ( Cmin kx,ky)=expj Dt( 2-Cmin 2 L) 1/2 = n=N P((,Cmin,n)(AL/Lmax)n=G(,cmin,kx,ky) n=O On précalcule et tabule les p(c,cmin,n) pour un échantillonnage de WEl(min Omaxl Pour une vitesse c quelconque,écrivons: Go(co,c,kx,ky) = exp j Dt(oe 2-c 2 L)l/2 n=N p(o, c,n) (AL/-max)n = G(o,c,kx,ky) n=O Go(w,c,kx,ky) = exp j Dtl 02-Cmin 2 (LC 2/Cmin 2) l 1/2 n=N E p(o,cmin,n) (C/Cmin)2 N (AL/Imax)n = G(co,c,kx,ky) n= O Donc p(a,c,n) = (c/cmin)Znp(o, cmin,n) fournit une approximation de G O (o,c,kx,ky). Les p(o,cmin,n) doivent être calculés sur lLin Cmax 2/cmint Lmx Cmax 2/Cmin 2 l au lieu de lLmin,Lmaxl De cette façon, la stabilité (le fait que le module ne soit jamais supérieur à 1), est assurée pour toutes les vitesses c et non pas seulement pour la vitesse tabulée La tabulation se fait avec Lc=o 2 sin 2 o/cmin 2, de sorte que cette valeur de Lc est plus grande que les valeurs de Lc nécessaires pour les autres vitesses C'est la raison pour laquelle il faut tabuler pour cmin plutôt que pour un co arbitraire. Si le rapport Cmax/cmin est trop grand, on peut tabuler pour 2 ou 3 valeurs de c: Cmint Cint:,Cint 2, et se baser pour calculer p(wc,n) sur la valeur tabulée juste inférieure à c. AVANTAGES DE L'INVENTION Le procédé conforme à la présente invention offre des avantages notables par rapport à la technique antérieure, en particulier par rapport au procédé proposé par Hale dans la référence l 4 l. Tout d'abord le procédé conforme à la présente invention, basé sur un développement polynomial et une synthèse spectrale de Laplacien, requiert seulement le calcul et l'application d'une banque de coefficients polynomiaux à une dimension p(w/c,n) et de deux filtres symétriques à une dimension d (n) et dy(n); alors que le procédé proposé par Hale, requiert le calcul et l'application d'une banque de filtres à une dimension b(co/c,n) et d'un filtre à deux dimensions M(n,,n,) (transformée de McClellan à 9 ou 17 termes par exemple). La méthode de la référence l 3 l requiert pour sa part le calcul et l'application d'une banque de filtres à deux dimensions a(oe/cn Xny) et est donc beaucoup plus lourde. Par ailleurs, nous avons montré que le calcul d'un développement polynomial se ramène à celui d'un spectre symétrique On peut donc utiliser pour le développement polynomial le même algorithme que pour le développement spectral Les approximations polynomiales et spectrales ont alors une précision similaire, pour un nombre de termes identiques. Enfin avec un nombre de termes identique, la synthèse spectrale du Laplacien donne une meilleure précision que la synthèse spectrale de la fonction coslDx(k 2 +ky 2)1/2 l dont la transformée de McClellan est un cas particulier Avec la méthode proposée, n'ayant que des spectres symétriques à une dimension à calculer, nous pouvons utiliser des algorithmes optimaux et rapides De plus, l'application d'un Laplacien, qui a une forme de croix dans le domaine (x,y) est plus rapide que l'application d'un coslDx(kx 2 +ky 2)1/2 l, qui a une forme rectangulaire. D'autre part, le cas Dxg Dy, qui est le plus courant en pratique, complique la méthode de Hale alors qu'elle ne pose aucun problème avec la méthode proposée. Un autre avantage notable de la méthode proposée est que l'extrapolation en temps est possible avec une table à une dimension, alors qu'elle nécessiterait une table à deux dimensions avec les méthodes existantes (ce qui la rend non pratiquable). La figure 4 représente une tranche verticale pour y=O d'une réponse impulsionnelle correspondant à: Dx=Dy=Dz= 10 m, c=l OO Om/s, 0 = 60 degrés, fmin= 2 Hz, fmax= 48 Hz, z= 900 m, N Lx=N Ly= 4 de sorte que le Laplacien comprend 17 coefficients Le développement polynomial est de degré 15. La figure 5 représente une tranche horizontale obtenue dans le cadre de la présente invention à l'aide d'un procédé mettant en oeuvre un développement polynomial à 15 termes et une synthèse de Laplacien à 17 termes, pour z= 450 m, soit un angle de 60 degrés. La figure 6 représente une tranche horizontale obtenue à l'aide d'un procédé classique comprenant un développement spectral à 12 termes et une transformée de McClellan à 17 termes, pour la même profondeur z= 450 m Le développement spectral à 12 termes est calculé par la méthode de synthèse spectrale par 2 optimisations L-w décrite précédemment. On peut noter que les perturbations visibles sur la figure 6, ne possèdent pas de symétrie circulaire Il est donc clair qu'elles ne résultent pas d'une erreur dans le développement à 12 termes, mais nécessairement dans la transformée de McClellan à 17 termes. Bien entendu la présente invention n'est pas limitée au processus qui vient d'être décrit mais s'étend à toutes variantes conformes à son esprit. Références bibliographiques: l 1 l O Holberg: Towards Optimum one-way Wave Propagation, Geophysical Prospecting, vol 36 n 2, pp 99- 114, Février 1988. l 2 l DHale: Stable, Explicit Depth Extrapolation of Seismic Wavefields, Geophysics, vol 56 n 11, pp 1770- 1777, Novembre 1991. l 3 l G Blaquiere, H W J Debeye, C P A Wapenaar, A J Berkhout: 3-D table-driven Migration, Geophysical Prospecting, vol 37 n 8, pp 925-958, Novembre 1989. l 4 l D Hale: 3-D Depth Migration via McClellan Transformations, Geophysics, vol 56 n 11, pp 1778-1785, Novembre 1991. l 5 l J H McClellan: The Design of Two-Dimensional Filters by Transformations, Proc 7th Annual Princeton Conf on Inform Sci and Syst, pp 247-251,1973. l 6 l J H McClellan, T W Parks: Equiripple Approximation of Fan Filters, Geophysics, vol 37 n 4, pp. 573-583, Août 1972. REVENDICATIONS 1 Procédé d'extrapolation d'un ensemble de traces enregistrées à l'aide d'un réseau de capteurs, dans un processus de prospection géophysique, caractérisé par le fait qu'il comprend les étapes consistant à: i) approximer un moins Laplacien Lo(kx, ky) = k X 2 +ky 2 par la somme de deux filtres à une dimension, et ii) approximer un opérateur d'extrapolation Fo par un polynome du moins Laplacien Lo. 2 Procédé selon la revendication 1, caractérisé par le fait que l'étape i) consiste à approximer le moins Laplacien par la somme de deux filtres à une dimension de la forme: n=N Lx n= N Ly L = j dx(n) cos(n Dx kx) + E dy(n) cos(n Dy ky) n=O n=O. 3 Procédé, selon l'une des revendications 1 ou 2, caractérisé par le fait qu'à l'étape ii) l'opérateur d'extrapolation est approximé par un polynome du moins Laplacien de la forme: n=N F(o/c, kx,ky)= E p(o/c,n) Lon n= 0. 4 Procédé selon l'une des revendications 1 ou 3, caractérisé par le fait que la synthèse du moins Laplacien consiste à: 1) Faire une synthèse spectrale de l'opérateur moins dérivée seconde en x, en calculant un filtre symétrique de longueur 2 N Lx+l, pour le pas d'échantillonnage Dx, dont le spectre Dx(k), définit par: n= N Lx Dx(k) = J dx(n)exp(j n Dx k) avec dx(-n)=dx(n) n=-NL, approxime le spectre de la moins dérivée seconde Do(k)=k 2, sur k El 0,axknyqxl, avec O rectangle kx El -axknyqx, xknyqxl, k El-ayknyqy, aykny l- Procédé selon l'une des revendications 1 à 4, caractérisé par le fait que les synthèses des coefficients dx(n) et dy(n) du moins Laplacien sont effectuées en utilisant l'algorithme de Remez, avec 50 (h)=h 2, W(h)=l pour h ElO,axrrl,W(h)=e pour h Elax Tr,l, o S est un petit nombre de l'ordre de 0,05. 6 Procédé selon l'une des revendications 1 à 5, caractérisé par le fait que la synthèse polynomiale de l'opérateur d'extrapolation est opérée sous forme d'une synthèse spectrale. 7 Procédé selon l'une des revendications 1 à 6, caractérisé par le fait que l'extrapolation en profondeur consiste à: calculer les coefficients dx(n), dy(n) du moins Laplacien, définir les bornes lLmn, Lmxl du moins Laplacien, pour chaque w/c parcourant un échantillonnage de l'intervalle loien/cmx, max/cminl calculer les coefficients p(w/c, n) du développement polynomial de l'opérateur d'extrapolation F O (w/c, L). 8 Procédé selon l'une des revendications 1 à 7, caractérisé par le fait que l'extrapolation en profondeur consiste à: 1) initialiser Wz+dz(mx,my) = plU/c(mx,my),Nl Wz(mx,m) puis 2) Pour n=N-1 jusqu'à O faire: Wz+dz(mxmy) = (L*Wz+dz)(m Xrmy) + plE/c(mrmy),nl Wz(mx, my) avec W'=L*W définit par n= N Lx N NL W'(nx,ny) = dx(n) W(nx-n,ny) + dy(n) W(nx,ny-n) n=-N Lx n=-N Ly. 9 Procédé selon l'une des revendications 1 à 8, caractérisé par le fait que la synthèse polynomiale de l'opérateur d'extrapolation est opérée sur la base des étapes consistant à: faire le changement de variable: L= 0,5 (Lmin-Lmax) cos(h) + 0,5 (Lmin+ Lmax)=g(h) avec LElL Qn,I Lxl, h El 0,l. Folg(h)l = So(h) étant un spectre symétrique connu. trouver les t(n), n ElO,Nl, tels que n=N S(h) = O t(n) cos(nh) n=O approxime So(h) pour h ElO,hcl et soit de module inférieur à 1 sur lhe, nl, avec he = g-l(Lc), puis définir F(L)=S(g-l(L)). Procédé selon l'une des revendications 1 à 9, caractérisé par le fait qu'il comprend l'étape consistant à écrire l'opérateur d'extrapolaion sous la forme: n=N F(L) = E t(n) Tnl( 2 LLmin-Lmax)/(min-lmax)l ( 4) n= O Tn(x) étant le polynôme de Chebychev de degré n définit par Tn(cos h) = cos nh. 11 Procédé selon l'une des revendications 1 à , caractérisé par le fait qu'il comprend l'étape consistant à écrire l'opérateur d'extrapolation sous la forme: n=N F(L)= 2 p(n)(AL/Lmx)n ( 5) n= O o A est un facteur d'échelle qui vaut typiquement 4. 12.Procédé selon l'une des revendications 1 à 11, caractérisé par le fait que l'extrapolation en profondeur consiste à: 1) initialiser W 2 +dz(mx,my) = plI/c(mx,my),Nl Wz(mx,mn) puis 2) Pour n=N-1 jusqu'à O faire: Wz+dz(m Xm Y) = (L*Wz+dz)(mx,my) + pl(/c(n,%my), nl Wz(mx,mn) avec W'=L*W définit par n= N Lx W' (nx,ny) = lAdx(n)/Lmaxl W(nx-nny) n=-N Lx n= N Ly + lAdy(n)/Lmaxl W(nx,ny-n) n=-N Ly. 13.Procédé selon l'une des revendications 1 à 12, caractérisé par le fait que l'extrapolation en profondeur consiste à: 1) initialiser: T O =Wz Ti=L*Wz Wz+dz(mx,my) = tluo/cz+dz(m X,my),Ol T O (mx,my) + tloe/CZ+d Z(m X,m Y),1 l T 1 (m Xm Y) 2) Pour n= 2 jusqu'à N faire: T= 2 *L*T 1-TO Wz+dz(mx,my) = Wz+dz(mx,my) + tl)/Cz+dz(mx,my),nl T(mx,my) To=Tl Ti=T Fin de la boucle en n avec W'=L*W définit par n= NLX W' (nx,ny) = l 2 dx(n)/(L=min-Lmaxl W(nx-nny) n=-Nta n= N Ly + l 2 dy(n)/(Lmin-Lmxl W(nx,ny-n) n=-N Ly (in+x)W(nxny)/ (I-L"ax). 14 Procédé selon l'une des revendications 1 à 13, caractérisé par le fait qu'il comprend en outre une phase d'extrapolation 3 D en temps sur la base de la fonction d'extrapolation qui s'écrit pour un pas d'extrapolation temporel Dt: Go(w,c,kx,ky)=expj Dt(W 2-c 2 Lo) 1/2 avec Lo(kx,ky)=kx 2 +ky 2 laquelle phase d'extrapolation 3 D en temps consiste à: tabuler les cofficients p(o, cmin, n) du développement polynomial pour au moins une valeur de vitesse déterminée (cmin), et calculer les autres coefficients requis p(w, c, n) par multiplications des coefficients tabulés par (c/Cmin)2 n, soit p(o,c,n) = (c/Cin)2 np(o,cmin,n). Procédé selon l'une des revendications 1 à 14, caractérisé par le fait qu'il comprend l'étape consistant à décomposer un opérateur d'extrapolation ( 20) en partie réelle et imaginaire, puis effectuer une synthèse optimale au sens de la norme L-o sur chacune des deux parties. 16 Procédé selon la revendication 15, caractérisé par le fait que la synthèse optimale L-w sur la partie réelle et la partie imaginaire de l'opérateur d'extrapolation est effectuée à l'aide de l'algorithme de Remez. 17 Procédé selon l'une des revendications 15 ou 16, caractérisé par le fait que la synthèse optimale de la partie réelle et de la partie imaginaire de l'opérateur d'extrapolation consiste à rechercher les a(n), n El 0,Nl du spectre synthétisé: n=N S(h)= a(n) cos(nh) n= O tels que Norme L-o de E(h) = maxhe Ol 0,,l IE(h) soit minimum, avec: So définissant pour h ECl 0,l, un spectre réel symétrique donné, W(h) définissant pour h El 0,,l, une fonction poids réelle positive donné, N étant la demi-longueur du filtre donnée, et E(h)=W(h)lS(h)-So(h)l étant la fonction d'erreur pondérée entre le spectre synthétisé et le spectre idéal donné. 18 Procédé selon l'une des revendications 15 à 17, caractérisé par le fait que l'opérateur d'extrapolation est synthétisé à l'aide des étapes consistant à: i) Décomposer l'opérateur d'extrapolation F O (k) exp(j 4 c) en parties réelles et imaginaires: Fo(k) exp(jc) = Ro(k) + j Io(k), : étant calculé de telle façon que RO(k) présente un maximum pour k=kc, ii) calculer ar(n), n ElO,Nl, tels que: n=N R(k) = Z ar(n) cos(n Dx k) n= O corresponde au minimum de ma Xk El 0,kcl IR(k)-R Po(k), tout en imposant aux extrema de R(k) sur lkc,knyql d'avoir des ordonnées +a, o a est un nombre plus petit que 1, iii) calculer la fonction de pondération W(k) pour k EEl kc, knyql: W(k)=lM 2 (k)-R 2 (k) l-1/2 o M(k) est le module maximum dans la bande lkc,knyql, définit par l'utilisateur, iv) calculer ai(n), n ElO,Nl, tels que: n=N I(k) = ai(n) cos(n Dx k) n= O O corresponde au minimum de maxkel 0, k Cl I(k)-Io O (k), tout en imposant que maxk Elkc,knyql W(k)II(k)l soit égal à 1, v) définir l'opérateur d'extrapolation par sommation des synthèses obtenues aux étapes ii) et iv), soit, F(k)=lR(k)+j I(k)l exp(-j c) les coefficients recherchés de l'opérateur étant: a(n) =lar(n)+jai(n)l exp(-j 4 c) 19 Procédé selon l'une des revendications 15 à 18, caractérisé par le fait qu'à l'étape i) on impose: R O (kc)= 1 et 10 (kc)= 0. Procédé selon l'une des revendications 15 à 19, caractérisé par le fait qu'à l'étape iii) on impose que le minimum de M(k) soit supérieur à a et que son maximum soit inférieur à 1. 21 Système pour prospection géophysique du type comprenant: une source artificielle ( 10) d'une onde sonore, un réseau de capteurs ( 20) disposés sur la surface du sol et sensibles aux ondes réfléchies vers le haut par les couches internes du sol, un enregistreur ( 30) relié à ces capteurs pour enregistrer les signaux électriques issus de ceux-ci, des moyens ( 40) de traitement des signaux enregistrés et des moyens ( 50) de visualisation d'un résultat exploitable après traitement de ces signaux, caractérisé par le fait que les moyens de traitement ( 40) sont adaptés pour décomposer un opérateur d'extrapolation ( 20) en partie réelle et imaginaire, puis effectuer une synthèse optimale au sens de la norme L-x sur chacune des deux parties. 22 Système selon la revendication 21, caractérisé par le fait que les moyens de traitement ( 40) effectuent la synthèse optimale L-x sur la partie réelle et la partie imaginaire de l'opérateur d'extrapolation à l'aide de l'algorithme de Remez. 23 Système selon l'une des revendications 21 ou 22, caractérisé par le fait que les moyens de traitement ( 40) effectuent la synthèse optimale de la partie réelle et de la partie imaginaire de l'opérateur d'extrapolation en recherchant les a(n), n ElO, Nl du spectre synthétisé: n=N S(h)= a(n) cos(nh) n= O tels que Norme L-o de E(h) = maxh Elo 0,,l IE(h) I soit minimum, avec: définissant pour h El 0,nl, un spectre réel symétrique donné, W(h) définissant pour h El 0,xl, une fonction poids réelle positive donné, N étant la demi-longueur du filtre donnée, et E(h)=W(h)lS(h)-So(h)l étant la fonction d'erreur pondérée entre le spectre synthétisé et le spectre idéal donné. 24 Système selon l'une des revendications 21 à 23, caractérisé par le fait que les moyens de traitement ( 40) synthétisent l'opérateur d'extrapolation à l'aide des étapes consistant à: i) Décomposer l'opérateur d'extrapolation Fo(k) exp(j 4 c) en parties réelles et imaginaires: F O (k) exp(j 4 ó) = Ro(k) + j Io(k), : étant calculé de telle façon que R O (k) présente un maximum pour k=kc, ii) calculer ar(n), n El 0,Nl, tels que: n=N R(k) = 2 ar(n) cos(n Dx k) n= O corresponde au minimum de maxkel 0,kcl IR(k)-R O (k)I, tout en imposant aux extrema de R(k) sur lkc,knyql d'avoir des ordonnées +a, o a est un nombre plus petit que 1, iii) calculer la fonction de pondération W(k) pour k E l kc, knyq l: W(k)=lM 2 (k)-R 2 (k) l-1/2 o M(k) est le module maximum dans la bande lkc, knyql, définit par l'utilisateur, iv) calculer ai(n), n ElO, Nl,,tels que: n=N I(k) = 2 ai(n) cos(n Dx k) n= O corresponde au minimum de maxk_lO,kcl II(k)-Io O (k)l, tout en imposant que maxkrlkc, knyql W(k)II(k)I soit égal à 1, v) définir l'opérateur d'extrapolation par sommation des synthèses obtenues aux étapes ii) et iv), soit, F(k)=lR(k)+j I(k)l exp(-j 4 c) les coefficients recherchés de l'opérateur étant: a(n) =lar(n)+jai(n)l exp(-j 4 ó) Système selon l'une des revendications 21 à 24, caractérisé par le fait qu'à l'étape i) les moyens de traitement ( 40) imposent: Ro(kc)=l et Io(kc)=O. 26 Système selon l'une des revendications 21 à , caractérisé par le fait qu'à l'étape iii) les moyens de traitement ( 40) imposent que le minimum de M(k) soit supérieur à a et que son maximum soit inférieur à 1. 27.Système selon l'une des revendications 21 à 26 caractérisé par le fait que les moyens de traitement ( 40) sont adaptés pour: i) approximer un moins Laplacien L O o(kx, ky) = kx 2 +ky 2 par la somme de deux filtres à une dimension, et ii) approximer un opérateur d'extrapolation Fo par un polynome du moins Laplacien Lo. 28 Système selon la revendication 27, caractérisé par le fait que à l'étape i) les moyens de traitement ( 40) approximent le moins Laplacien par la somme de deux filtres à une dimension de la forme: n=N Lx n= NLY L = dx(n) cos(n Dx kx)+ dy(n) cos(n Dy ky) n= O n=O. 29 Système selon l'une des revendications 27 ou 28, caractérisé par le fait que à l'étape ii) les moyens de traitement ( 40) approximent l'opérateur d'extrapolation par un polynome du moins Laplacien de la forme: n=N F(o/ckx,ky)= S p(o/c,n) L O n n= 0. Système selon l'une des revendications 21 à 29, caractérisé par le fait que pour la synthèse du moins Laplacien les moyens de traitement sont adaptés pour: 1) Faire une synthèse spectrale de l'opérateur moins dérivée seconde en x, en calculant un filtre symétrique de longueur 2 N Lx+l, pour le pas d'échantillonnage Dx, dont le spectre Dx(k), définit par: n= N Lx Dx(k) = dx(n)exp(j n Dx k) avec dx(-n)=dx(n) n=-N Lx approxime le spectre de la moins dérivée seconde Do(k)=k 2, sur k(l 0,xknyqxl, avec O rectangle kx El-xknyqx knyqxl, ky El -yknyayknyqyl- 31 Système selon l'une des revendications 21 à , caractérisé par le fait que les synthèses des coefficients dx(n) et dy(n) du moins Laplacien sont effectuées par les moyens de traitement ( 40) en utilisant l'algorithme de Remez, avec 50 (h)=h 2, W(h)= 1 pour h ElO,axel,W(h)=E pour h Elaxir,rl, o S est un petit nombre de l'ordre de 0,05. 32 Système selon l'une des revendications 21 à 31, caractérisé par le fait que la synthèse polynomiale de l'opérateur d'extrapolation est opérée par les moyens de traitement ( 40) sous forme d'une synthèse spectrale. 33 Système selon l'une des revendications 21 à 32, caractérisé par le fait que les moyens de traitement ( 40) opèrent l'extrapolation en profondeur à l'aide des étapes qui consistent à: calculer les coefficients dx(n), dy(n) du moins Laplacien, définir les bornes lL Min, Lmxl du moins Laplacien, pour chaque w/c parcourant un échantillonnage de l'intervalle loin/cmax, oemax/cmnl calculer les coefficients p( d'extrapolation F O (w/c, L). 34 Système selon l'une des revendications 21 à 33, caractérisé par le fait que pour l'extrapolation en profondeur les moyens de traitement ( 40) sont adaptés pour 1) initialiser Wz+dz(mx, my) = plw/c(mx,my),Nl Wz(mx,my) puis 2) Pour n=N-1 jusqu'à O faire: Wz+ dz(mx,my) = (L*Wz+dz)(mx,my) + pl)/c(n,my),nl Wz(mx,my) avec W'=L*W définit par n= N Lx n= NL W' (nx,ny) = dx(n) W(nx-n,ny) + dy(n) W(nx,ny-n) n=-N Lx n=-N Ly. Système selon l'une des revendications 21 à 34, caractérisé par le fait que la synthèse polynomiale de l'opérateur d'extrapolation est opérée par les moyens de traitement ( 40) sur la base des étapes consistant à: faire le changement de variable: L= 0,5 (Lmin-Lmax) cos(h) + 0,5 (Lmin+Lmax)=g(h) avec LElLmin,Lmaxl, h El 0,Tl- F O lg(h)l = So(h) étant un spectre symétrique connu. trouver les t(n), n ElO,Nl, tels que n=N S(h) = X t(n) cos(nh) n=O approxime So(h) pour h El 0,hcl et soit de module inférieur à 1 sur lhc,,l, avec h = g-l(Lc), puis définir F(L)=S(g-l(L)). 36 Système selon l'une des revendications 21 à , caractérisé par le fait que les moyens de traitement ( 40) écrivent l'opérateur d'extrapolaion sous la forme: n=N F(L) = 1 t(n) Tnl ( 2 L-Lmin-Lax)/(Lmn-Lmax)l ( 4) n= O Ta(x) étant le polynôme de Chebychev de degré n définit par T,(cos h) = cos nh. 37 Système selon l'une des revendications 21 à 36, caractérisé par le fait que les moyens de traitement ( 40) écrivent l'opérateur d'extrapolation sous la forme: n=N F(L)= X p(n)(AL/I Lax)n ( 5) n= O o A est un facteur d'échelle qui vaut typiquement 4. 38.Système selon l'une des revendications 21 à 37, caractérisé par le fait que pour l'extrapolation en profondeur les moyens de traitement ( 40) sont adaptés pour: 1) initialiser Wz+dz(mx, my) = ploe/c(mx,my),Nl Wz(mx,my) puis 2) Pour n=N-1 jusqu'à O faire: Wz+ dz(mx,my) = (L*Wz+dz)(mx,my) + pl)/c(mx,my),nl Wz(mx,my) avec W'=L*W définit par n= NLX W'(nx,ny) =S lAdx(n)/L mxl W(nx-n,ny) n=-N Lx n NLY + X lAdy(n)/Lmaxl W(nx,ny-n) n=-N Ly. 39 Système selon l'une des revendications 21 à 38, caractérisé par le fait que pour l'extrapolation en profondeur les moyens de traitement ( 40) sont adaptés pour: 1) initialiser: To=Wz Ti=L*Wz Wz+dz(mx,my) = tlCe/cz+d Z(mx,my),Ol T O (mx,my) + tlC/cz+dz(mx, my),1 l Tl(mx,my) 2) Pour n= 2 jusqu'à N faire: T= 2 *L*T 1-To Wz+ dz(mx,my) = Wz+dz(mx,my) + tlco/cz+dz(mx,my),nl T(mx,my) T O T 1 T 1 T Fin de la boucle en n avec W'=L*W définit par n= N Lx W (nx,ny) = l 2 dx(n)/(L Imin-Lmaxl W(nx-nny) n N Lx n= N Ly + l 2 dy(n)/(Lmin-Lmaxl W(nx,ny-n) n=-N Ly (Lijn+Inax)W(nx, ny)/ (Imnin 4 znax). 40.Système selon l'une des revendications 21 à 39, caractérisé par le fait que les moyens de traitement ( 40) sont adaptés pour opérer une phase d'extrapolation 3 D en temps sur la base de la fonction d'extrapolation qui s'écrit pour un pas d'extrapolation temporel Dt: Go(oe,c,kx,ky) = exp j Dt(w 2-c 2 LO)1/2 avec L O (k X,ky)=kx 2 +k Y 2 grâce aux étapes qui consistent à: tabuler les cofficients p(w, cmin, n) du développement polynomial pur au moins une valeur de vitesse déterminée (Cmin), et calculer les autres coefficients requis p(w, c, n) par multiplications des coefficients tabulés par (c/cin)2 n, soit p(o,c,n) = (c/cmin)2 np(oecminn).