CA2177180C - Procede de prospection sismique avec application d'un filtre d'erreur de prediction auto-deconvolue - Google Patents

Abstract

The invention relates to a seismic prospection method wherein a seismic shock is caused in the earth subsoil and seismic data is collected by means of sensors, said seismic data being processed in order to draw useful information concerning the geology of the subsoil, the seismic data containing a signal y¿o? intended to be isolated, said signal y¿o? being present in the seismic data embedded in an additive signal. The method comprises the following steps: 1) calculating for a given integer p a prediction error filter a having p coefficients and cancelling at most the signal y¿o?; 2) applying to the seismic data the self-deconvoluted prediction error filter in order to obtain the filtered data wherein the signal y¿o? is absent; 3) subtracting from the initial data the filtered data in order to obtain the processed data containing the signal y¿o? without the additive signal.

Description

2 ~ 7 718 0 p~~~p1231 PROCEDE DE PROSPECTION SISMIQUE AVEC APPLICATION D'UN FILTRE
D'ERREUR DE PREDICTION AL1T0-DECONVOLUE
La présente invention concerne la prospection sismique, et notamment un procédé destiné à atténuer les bruits affectant les traces S sismiques.
Le principe général de la prospection sismique consiste à
provoquer, à l'aide d'une source sismique, un ébranlement dans le sous-sol et à~ enregistrer à l'aide de capteurs (géophones ou hydrophones) des données sismiques pour en tirer une information sur la géologie du sous-sol, en particulier détecter la présence d'hydrocarbures.
Les capteurs sont habituellement régulièrement espacés selon tes noeuds d'une grille spatiale et les enregistrements y(t) en fonction du temps t des signaux recueillis par ces capteurs, encore appelés traces sismiques, sont regroupés et juxtaposés selon les valeurs de la coordonnée x (coordonnée spatiale) des capteurs sur la grille, pour former une section sismique t-x constituée de données sismiques y(t, x).
Les données sismiques comprennent une information utile (par exemple une succession d'échos en sismique réflexion) noyée dans un bruit de fond global que l'on cherche à éliminer par un traitement approprié, qui doit en outre tenir compte de la possibilité de traces manquantes sur la section sismique. En effet, il arrive dans la pratique qu'en certains points de la grille spatiale on ne dispose d'aucun enregistrement utilisable, soit parce que la mise en place d'un capteur en cet emplacement n'était pas possible en raison de la nature du terrain, soit à cause d'une défaillance de l'installation d'acquisition des données, ou soit encore parce que le signal délivré par un capteur a été saturé pendant l'enregistrement par un bruit très fort.
Le traitement des données sismiques, pour éliminer les bruits, peut s'effectuer dans un plan f-x sur des données Y(f,x) après transformée de Fourier des données sismiques initiales y(t,x) pour la variable temps t, ou dans un plan f-k sur des données Y(f,k) après transformée de Fourier des données Y(f,x) pour la variable spatiale x, où f désigne une variable fréquentielle et k une variable nombre d'onde.
On distingue habituellement les bruits présentant un caractère cohérent, qui sont engendrés en général par des phénomènes physiques identifiés tels qu'un certain mode de propagation des ondes sismiques (cas des ondes se propageant en surface par exemple) et que l'on souhaite éliminer parce que ne contenant aucune information utile sur le sous-sol, ou encore wo 96/09562 PCT/FR95/01231 parce que l'extraction de l'information utile qu'ils contiennent serait trop complexe (cas des réflexions dites multiples), et les bruits présentant un ~
caractère non cohérent, aléatoire, dont l'origine physique n'est en général pas identifiée (par exemple dûs à des phénomènes naturels tels que la houle ou les microséismes).
On a proposé de nombreux procédés pour traiter les données sismiques en vue d'atténuer les bruits, et ces procédés peuvent se classer dans plusieurs grandes catégories, selon qu'ils sont destinés plus particulièrement à éliminer des bruits ayant un caractère localisé ou non.
Dans un premier type de procédé, encore appelé procédé
d'atténuation "déterministe", on se place classiquement dans le plan f-k et on espère que tes bruits affectant les données sismiques Y(f,k) pourront être identifiés et localisés dans une région particulière de ce plan. A titre indicatif, un bruit se propageant avec une vitesse constante se traduit par une droite dans le plan f-k et pourra être éliminé par un simple filtrage portant sur la variable nombre d'onde k. Plus généralement, on traite les données dans le plan f-k en les multipliant par un filtre dit "projectif", valant en principe 0 dans les régions où le bruit est supposé présent (domaine affaibli du filtre) et 1 ailleurs (domaine passant du filtre). Une transition targe doit toutefois être prévue entre les domaines passant et affaiblis du filtre pour limiter l'extension de ce dernier dans la section sismique t-x, sous peine d'être confronté à des effets de bord lorsque le filtre sort de la section sismique et que les données manquantes sont prises égales à
zéro (ce qui revient à annuler certains coefficients du filtre et donc à
en détériorer les performances). Les effets de bord sont particulièrement gênants dans le domaine spatial pour lequel le nombre d'échantillons, qui correspond au nombre de traces, est relativement faible. Le choix d'une transition large, qui correspond à une faible extension spatiale, conduit alors à un filtre insuffisamment sélectif, c'est-à-dire que pour atténuer un bruit donné dans le plan f-k on affaiblira également une partie non négligeable du signal utile qui se trouve à proximité. On a proposé, dans le brevet américain US-5 067 112, de remédier aux inconvénients du filtrage classiquement opéré dans le plan f-k en effectuant un filtrage par nombre d'onde dans le plan f-x. Le procédé décrit ne concerne cependant que le cas particulier où pour chaque fréquence du plan f-x il n'y a qu'un seul nombre d'onde kp à enlever, et d'autre part le filtre utilisé, non récursif, présente le double inconvénient d'être long à appliquer et de s'accompagner d'effets de WO 96109562 217 718 0 p~~9g~01231 bord lorsque l'on recherche une grande slectivit en nombre d'onde.

Un second type de procd, encore appel procd d' attnuation "statistique", s'applique au cas o l'on ne dispose pas d'informations sur la localisation du bruit dans la section sismique. On a propos, dans l'article Expanded Abstracts of Tbe 54th Annual SEG Meeting, Atlanta, Canales, L. L., 1984, Random Noise Reduction, d'effectuer un filtrage dit "prdictif' dans le plan f-x. Selon l'auteur de l'article prcit, le bruit affectant une section sismique peut tre isol du reste des donnes sismiques en appliquant ces dernires un filtre d'erreur de prdiction pour la variable spatiale x. Les coefficients de ce filtre sont choisis de faon minimiser l'nergie du bruit rsultant de l'application du filtre aux donnes. Le filtrage qui est propos par Canales n'est pas de type projectif, car il n'annule pas les donnes supposes ne contenir que du bruit, et il ne peut valoir 1 sur une bande de frquences donne. Or, tout procd d'attnuation des bruits prservant les signaux utiles devrait tre de type projectif . lorsque les donnes ne contiennent que du bruit, elle doivent tre multiplies par 0 pour tre limines, et dans le cas contraire elles doivent tre multiplies par 1 pour tre conserves intactes. Dans le procd propos par Canales, les donnes sismiques subissent une dgradation importante.

Un troisime type de procd, propos notamment dans US 5,182,729 (Duren et al.) consiste sparer un signal suppos parfaitement prdictible d'un bruit additif en utilisant une matrice de projection.
Ce type de procd est bas sur les considrations suivantes Soit N donnes sismiques y(0),..,y(N), o N est un nombre entier, correspondant la somme d'un signal yp(0),..,yp(N) parfaitement prdictible et d'un bruit et e(0),..,e(N).
La prédictibilité de yp - [yp(0),yp(1),..,yp(N)]T (T désignant l'opération de transposition matricielle) pour un filtre d'erreur de prédiction a à p+1 coefficient a(0), a(I),..,a(p) où p est un nombre entier inférieur à
N, se traduit par le fait que yp appartient au sous-espace Ayo=0, A étant la matrice à N+1 colonnes et N+1-p lignes a(p) ... a(0) 0 ... ... 0 0 a(p) ... a(0) 0 ... 0 3 5 A = ... ... ... ... ... ... ... (I) 0 ... ... 0 a(p) ... a(0) i On peut alors estimer yo en projetant [y(0).y(1),..,y(N)]T sur ce sous-espace. Il est connu en algèbre linéaire que ceci s'écrit yp = (I-A* (AA*)-tA)y (II) e = A* (AA*)-~ Ay (III) où I est la matrice unité; A* étant la matrice transposée conjuguée de A, e étant le bruit additif estimé.
L'inversion de la matrice AA* fait alors appel au calcul matriciel, ce qui limite très sérieusement l'intérêt pratique de ces procédés. Le coût en temps de calcul et en place mémoire d'une inversion de matrice de dimension égale à celle des données est en général prohibitif dans le cas du traitement sismique, qui traite des quantités énormes de données.
Dans le brevet Duren et al., il est appliqué une matrice de projection du type de celle de l'équation III. Ceci est viable dans la mesure où
il considère un cas particulier où le nombre de données N est de 3 à 5.
L'invention vise à remédier aux inconvénients des procédés de l'art antérieur.
Elle n'utilise que des procédés de type filtrage sur les données, bien plus économique que les procédés de type matriciel.
Les procédés de type filtrage se distinguent des procédés de type matriciels par le critère suivant . dans les procédés de type matriciel, le calcul des coefficients de la matrice à appliquer aux données nécessite de connaître le nombre d'échantillons et de calculer et stocl'er pour chaque échantillon des données, un ensemble de coefficients. Les procédés de type filtrage se caractérisent par le fait qu'ils n'ont pas besoin de connaître le nombre d'échantillons de données au moment du calcul des coefficients. Les coefficients sont calculés une fois pour toute, pour toutes les données. .
On comprendra que pour des traitements sismiques, dans lesquels on manipule une grande masse de données, les procédés de type filtrage sont .
particulièrement avantageux par rapport aux méthodes matricielles.
L'invention enseigne comment réaliser une projection approchée par filtrage.
Les filtrages souffrent en outre de l' inconvénient des effets de ~. wo mo9s6z 21 7 71 8 0 bords. C'est un autre enseignement de cette invention que de montrer comment s'en affranchir, sans augmenter le coût de calcul, en utilisant la structure de filtrage non-stationnaire.
Par filtrage on entend par conséquent dans tous le présent texte 5 une opération consistant à calculer des échantillons de sortie à partir de combinaisons linéaires des échantillons d'entrée avec un ensemble de coefficients calculés préalablement et identiques pour les différents échantillons d'entrée.
Dans le cas où un échantillon de sortie est calculé à partir des échantillons d'entrée uniquement, on parle de filtrage non récursif.
Dans le cas où un échantillon de sortie est calculé à partir de combinaisons linéaires des échantillons d'entrëe et des échantillons de sortie précédemment calculés, on parle de filtrage récursif.
Par filtrage non-stationnaire, on entend un filtrage dont les coefficients sont non stationnaires et recalculés à chaque fois pour les premiers échantillons d'entrée, jusqu'à converger sur des valeurs de ceofficients utilisés ensuite de façon stationnaire pour les autres échantillons d'entrée. Le nombre d'ensembles de coefficients à calculer est supérieur à 1 mais reste petit et est de toute façon indépendant du nombre d'échantillons de données.
Par opération matricielle sur les données, on entend une opération consistant à calculer autant d'ensemble de coefficients qu'il y a d'échantillons dans les données.
L'invention propose quant à elle un proocédé de prospection sismique dans lequel - on provoque un ébranlement sismique dans le sous-sol, - on recueille à l'aide de capteurs des données sismiques échantillonnées y= [y(0),...,y(N)]T, où N est un nombre entier, contenant un signal yo que l'on cherche à isoler noyé dans un bruit additif, - on réalise dans le domaine temporel ou fréquentiel un filtrage des données sismiques, pour obtenir des données filtrées dans lesquelles le signal à
isoler est absent, - on soustrait aux données sismiques initiales les données filtrées pour obtenir des données traitées yo(0),...,yo(N) correspondant au signal yp dépourvu de bruit additif, - on traite les données ainsi acquises pour en déduire une information utile wo Wo956z PCT/FR93/01231 2;77180 sur la géologie du sous-sol, le filtre utilisé pour le filtrage des données sismiques étant un filtre fréquentiel d'erreur de prédiction autodéconvolué M(f) tel que M(f) = IA(f)I?/R(f), , ou un filtre équivalent dans le domaine temporel, A(f) étant le spectre d'un filtre a d'erreur de prédiction à p+1 coefficient a(0), ...,a(p), où p est un nombre entier inférieur à N, a étant préalablement choisi pour annuler au mieux par convolution le signal yp que l'on cherche à
isoler, R(f) étant une autocorrélation précolorisée de ce filtre de prédiction a vérifiant R(fj = IA(f)I? + 82 IB(f)I?, où ~ et B(f) sont respectivement un facteur et un filtre de précolorisation préalablement choisis en fonction de la sélectivité souhaitée pour le filtrage.
Comme on le comprendra par la suite, un tel filtrage est approximativement projectif.
En outre, le procédé permet de réaliser une atténuation "déterministe" avec un filtre sélectif sans rencontrer d'effets de bord et une atténuation "statistique" par filtrage de type projectif.
L'invention sera mieux comprise à la lecture de la description détaillée qui va suivre, d'exemples non limitatifs de mise en oeuvre du procédé selon l'invention, et .à ~ l'examen du dessin annexé sur lequel - la figure 1 illustre, sous une forme schématique, le procédé selon l' invention, - la figure 2 représente l'évolution en fonction de la fréquence du module du spectre d'un filtre d'erreur de prédiction, - la figure 3A représente le spectre du filtre d'erreur de prédiction auto-déconvolué correspondant à la figure 2, pour une valeur donnée d'un paramètre de calcul E, - la figure 3B illustre l'évolution du filtre d'erreur de prédiction auto-déconvolué représenté sur la figure 3A, pour une valeur différente du paramètre de calcul ~, - la figure 4 représente une section sismique synthétique, - la figure 5 représente la partie prédictibte yo(t,x) des données sismiques y(t,x) représentées sur la section sismique de la figure 4, obtenue par la wo ~s~o~ss2 2 i 7 718 0 rcr/»sioi2m mise en oeuvre du procédé selon l'invention, - la figure 6 représente la partie non prédictible e(t,x) des données sismiques y(t,x) représentées sur la section sismique de la figure 4, . obtenue par la mise en oeuvre du procédé selon l'invention, - la figure 7 représente la partie prédictible des données sismiques y(t,x) représentées sur la section sismique de la figure 4, obtenue par la mise en oeuvre d'un procédé selon l'art antérieur, - la figure 8 représente la partie non prédictible des données sismiques y(t,x) représentées sur la section sismique de la figure 4, obtenue par la mise en oeuvre d'un procédé selon l'art antérieur, - la figure 9 représente la section sismique de la figure 4 avec 10 traces manquantes, - la figure 10 représente une section sismique restaurée par la mise en oeuvre du procédé selon l'invention, - la figure 11 représente la partie prédictible des données sismiques restaurées de la figure 10, et - la figure 12 représente la partie non prédictible des données sismiques restaurées de la figure 10.
Dans la description qui va suivre, on applique l'invention à des données sismiques y(t) fonction de la variable temps t. Les données sismiques y(t) sont supposées enregistrées sous forme numérique et discrétisées avec un pas d'échantillonnage 0'C=1/fe, où fe est la fréquence d'échantillonnage.
Dans la description qui suit, 0'L est pris égal à 1 pour simplifier les notations.On peut bien entendu appliquer de la même façon l'invention à des données sismiques y(x) fonction de la variable espace x.
L'invention permet, lorsqu'appliquée aux données y(t), d'isoler des signaux de bruit de fréquence telle que 50 Hz ou 60 Hz affectant les données sismiques, en vue de leur élimination, et lorsqu'appliquée aux données y(x), d'enlever des perturbations dues par exemple à la configuration spatiale des capteurs sur la grille.
La variable entière n est utilisée dans la suite de la description - pour désigner les échantillons de données sismiques discrétisées. Lorsque l'invention est appliquée à des données sismiques fonction de la variable temps t, n est compris entre 0 et N, le nombre total d'échantillons enregistrés par un capteur pendant la durée de la mesure étant égal à N+1. Lorsque l'invention est appliquée dans la section sismique t-x ou le plan f-x considéré

wo 96/09562 PCT/FR95/01231 g à des données sismiques fonction de la variable espace x, n est compris entre 0 et le nombre total moins un de traces juxtaposées dans la section sismique t-x considérée.
On considère, de façon générale, dans l'exposé qui suit, que l'on dispose de données sismiques discrétisées y(n) extraites d'une section sismique obtenue de façon connue en soi (étape 1 sur la figure 1), pouvant se décomposer en une composante prédictible yo(n) que l'on cherche à isoler et une composante non prédictible e(n) Y(n) = Yo(n) + e(n) (1) Conformément à l'invention, on calcule un filtre d'erreur de prédiction a de longueur p, c'est-à-dire ayant p+I coefficients, annulant au mieux yo(n) (étape 2 sur la figure 1). On considère ici, pour faciliter la compréhension de l'invention, que le filtre d'erreur de prédiction a annule complètement ya(n), c'est-à-dire que les p+1 coefficients a(0), a(1), ..., a(p) du filtre a sont tels que a(0) yo(n) + a(1) yo(n-1) + ... + a(p) yo(n-p) = 0, avec a(0) = 1 (2) Dans le domaine fréquentiel, obtenu en prenant la transformée de Fourier de y(n) pour la variable temps, par exemple par une méthode de transformée de Fourier rapide réalisant l'opération Y(m') = E y(m)e2Jn mm'/M m'= 0,...,M-1 m=0 où M est un entier choisi supérieur aux nombres d'échantillons temporels du signal, les équations (1) et (2) s'écrivent, en utilisant des lettres majuscules pour désigner e(n) et yo(n) après transformée de Fourier Y(~ = Yo(~ + E(~ (l') A(~ 1'0(~ = 0 (2') où A(f) est le spectre du filtre d'erreur de prédiction n.
On déduit de (l') et (2') A(~ E(~ = A(17 Y(17 (3') Soit D(f) un opérateur de déconvolution pour A(f).
On applique aux données Y(f), conformément à l'invention, le filtre d'erreur de prédiction auto-déconvolué D(f) A(f) (étape 3 sur la figure 1) et l'on obtient E(17 = D(17 A(17 Y(i. (4') wo ~ro~2 217 718 0 On peut ensuite (étape 4 sur la figure I) isoler le signal utile Yo(f) ~ par Yo(f) = Y(f) - E(f) Dans la pratique, le filtre d'erreur de prédiction a n'est jamais rigoureusement inversible et l'on ne peut trouver un opérateur de déconvolution D(f) tel que D(f) A(f) = 1 ; l'opérateur de déconvolution D(f) est donc dans la pratique un opérateur inverse approché, déterminé par exemple en minimisant ID(f)A(f) - 11'- + >r'- IB(f)l'- ID(f)12 Cette expression est minimale pour ' D(f)= A(f)/R(f) (5') avec R(f) = IA(f)IZ + ~2 IB(fji'- (6') Dans cette expression, désigne le complexe conjugué, ~ et B(f) contrôlent le calcul et sont choisis par l'homme du métier en fonction de la marge d'erreur que l'on s'autorise dans le calcul de D(f). ~ est un scalaire appelé
"facteur de précolorisation", B(f) est un filtre appelé "filtre de précolorisation". Dans le domaine temporel, le filtre de préeolorisation b est de longueur finie q, préférentiellement inférieure à p, avec pour coefficients b(0). ..., b(q) et b(0)=1. Lorsque B(f) = 1, >r est appelé
"facteur de préblanchiment".
L'équation (4') s'écrit encore E(~ = M(~ Y(>7 (7') avec M(f) = IA(f)I? / [IA(f)l'- + >r2 IB(f)l'-] (8') soit M(f) = IA(f)l'- / R(fj (9' ) R(f) est un filtre appelé "autocorrélation précolorisée" du filtre d'erreur de prédiction a.
Ainsi, le filtre M(f) obtenu et appliqué conformément à
l'invention aux données sismiques Y(f) est analogue à un filtre de type projectif, car M(f) = 0 si A(f) = 0 M(f) = 1 si IA(f)I » ~ IB(f)I

wo 96ro9562 ~ PCT/FR95101231 lo Le calcul du filtre d'erreur de prédiction auto-déconvolué M(f) à
appliquer aux données sismiques nécessite la connaissance du filtre d'erreur de prédiction A(f) et de l'autocorrélation précolorisée R(f) d'après l'équation (9'). Après avoir calculé le filtre d'erreur de prédiction a et l'autocorrélation précolorisée, le filtre M(f) à appliquer aux données dans le domaine fréquentiel résulte de la même équation (9'). Il sera détaillé plus loin comment l'appliquer dans le domaine temporel.
On rappelle, que de façon connue en soi par l'homme du métier, la convolution de deux séquences x(M1), ..., x(M2) et y(N1), ..., y(N2) est la séquence z(M1+N1), ...., z(M2+N2) telle que z(i) = E xÜ) y(i-j) j Nous utiliserons par la suite la notation:
z(i)=x(i)*y(i) pour désigner cette opération de convolution.
On désignera l'opération de corrélation qui est la convolution de la séquence x(i) avec la séquence retournée conjuguée yr(i) telle que:
yr(i)=y(-i) pour i = -N2, ..., -N1 par la notation z(i)=x(i)*y(-i) 8(i) dénotera par la suite une séquence telle que S(0)=1 et 8(i)=0 si i est différent de 0.
CALCUL DETERMINISTE DU FILTRE D'ERREUR DE PREDICTION ET DE
L'AUTOCORRELATION PRECOLORISEE
Dans une première mise en oeuvre du procédé selon l'in~~ention, les p+1 coefficients du filtre d'erreur de prédiction a sont calculés à partir de paramètres physiques connus du signal yo que l'on cherche à isoler, et de préférence, les p+1 coefficients du filtre d'erreur de prédiction n s o n t calculés à partir du support du spectre fréquentiel Yo (f) du signal yo (ensemble des valeurs de f pour lesquelles Yo(f) ~ 0).
Le filtre d'erreur de prédiction n peut être calculé de plusieurs façons, lorsque l'on connaît le support du spectre du signal prédictible yo que l'on désire isoler des données sismiques y(n), ce qui est le cas par exemple wo 96/09562 ~ ~ ~ ~ ~ ~ ~ PCT/FR95/01231 lorsque l'on désire éliminer des données sismiques le bruit de fréquence a 50 Hz ou 60 Hz dû aux installations électriques. On traduit dans l'invention, en ce qui concerne le filtrage déterministe, l'hypothèse qu'un signal est à
- bande limitée en terme de, quasi-prédictibilité de ce signal.
S SiEnal v" présentant un spectre de raies isolée ~ Dans cet exemple, le signal yo à isoler dans les données sismiques y(n) est supposé présenter un spectre Yo(f) constitué de p raies fq (q = 1, ..., p), On définit, pour chaque raie de fréquence fq, un filtre d'erreur de prédiction élémentaire aq de longueur l, à deux coefficients a9 (0) et aq(1) égaux à
aq(0) = 1 aq(1) _ _ e-'-jafq Les coefficients a(0), a(1), ..., a(p) du filtre d'erreur de prédiction global a à appliquer aux données sismiques y(n) sont alors définis par a(i) = at(i) * a~(i) * ... * ap(i) De même, on calcule pour chaque entier q l'autocorrélation précolorisée rq avec préblanchiment, de longueur 2, à trois coefficients i = -1, 0, 1 rq(i) = aq(i)* aq(-i) + E2 S(i) puis on calcule r(i) = rt(i)*r2(i)*...*rp(i) a(i) et r(i) permettent, après transformée de Fourier, de construire un filtre M(f) d'après l'équation (9') qui rejette les p fréquentés fq. ~ contrôle la sélectivité du filtre ; plus E est petit, plus le filtre est sélectif.
Si pal y~, présentant un spectre Y~,(fl non nul pour au moins un intervalle de fréguences ~ Dans cet exemple, le signal yo à isoler est supposé présenter un spectre Yo(f) non nul pour au moins un intervalle de fréquences.
Mthode par moindrescarrs Dans cette premiremthode on dtermine les coefficientsa(0), a(1), ..., a(p) du filtred'erreurde prdiction a, de spectre A(f),une par mthode des moindres carrs, c'est--dire qu'on les choisit que tels l'intégrale J~p W(f) IA(f)1? df soit minimum, sachant que a(0) = 1 et que W(f) wo 9sio~2 2 ~ 7 ~ 18 0 PCT/FR~S/01231 est une fonction qui vaut 1 lorsque la fréquence f est dans un intervalle de fréquences où le spectre Yo(f) du signal yo à isoler est non nul, et 0 ailleurs.
Minimiser l'intégrale définie plus haut revient à résoudre, dans le domaine temporel, après transformée de Fourier inverse de W(f) pour obtenir w(i) dans le domaine temporel, de façon connue de l'homme du métier, un système de Toeplitz formé à l'aide de w(-p),...,w(p), par exemple au moyen d'un algorithme de Levinson sans second membre.
La matrice de Toeplitz M de coefficients m(i,j) se calcule par m(i,j) = w(i-j) pour i,j = 0, ..., p et on résout, de façon connue en soi P
E m(i~j) aU) = a2 s(i) j=0 avec a(0) = 1, a un scalaire.
I1 est préférable, pour le calcul de l'autocorrélation précolorisée par l'équation (6'), d'utiliser B(f)=1 et un facteur de préblanchiment ~ qui soit grand devant dans lesdits intervalles les valeurs de du module IA(f)) frquences.

Mth ode avec lvnmes de Chebychev ro En variante, prfrentiellement le cas o le procd et dans selon l'invention appliqu traitement de sismiques y(t) ou est au donnes y(x) ne dpendant que de la seulevariable temps t ou espace x, on dcompose dans une deuxime mthode support du spectreYo(f) en Q intervalles le de 2 S fréquence [fq - ~fq, fq + Ofq]
0 fq étant la demi-largeur du sous-intervalle centré sur la fréquence fq.
On définit les Q fonctions Sq(f) = e-JnPf TPq[sin(7t(f-fq)) / sin(JC ~fq)]
où TPq est le polynôme de Chebychev de première espèce et d'ordre pq, défini récursivement pour une variable a par To(a) = 1, T1(a) = a, TP(a) = 2aTP_1 (a) - TP_~(a) On peut vérifier que ~''~ wo ~'°~2 21 7 71 8 0 pq Sy(f) = E sq(m) e-J2xmf m=0 C'est-à-dire qu'après transformée de Fourier inverse de Sq(f), on obtient dans le domaine temporel, le signal sq(i) dont les seuls coefficients non nuls sont sq(0), ..., sq(pq).
On calcule ensuite des coefficients de filtres d'erreur de prédiction élémentaires aq (après avoir posé ~,q = 1/ sq(0)) par : aq(i) = sq(i) ~,q pour i = 0, ..., pq et le filtre d'erreur de prédiction a par a(i) = at(i) * a~(i) * ... * aQ(i) On calcule les autocorrélations précolorisées élémentaires rq (i), avec Eq » %1, q par rq(i) = aq(i) * aq(-i) + ~'- s(i), puis l'autocorrélation précolorisée r(i) par r(i) = rl(i) * r2(~) * ... * rQ(i).
On a tracé sur la figure 2 le module IA(f)I du spectre A(f) en fonction de la fréquence, pour Yo(f) de support centré sur fl - 1/10 avec D f 1 - 1/20. On remarquera que le filtre d'erreur de prédiction a obtenu présente un spectre de module IA(f)I petit dans l'intervalle de fréquence du support du spectre Yo(f).
On a respectivement tracé sur les figures 3A et 3B le spectre M(f) pour deux valeurs du facteur de préblanchiment ~ respectivement égales à
10-2 et 10-3. On remarquera que le paramètre E permet de contrôler la pente du filtre M(f).
On décrira plus loin comment calculer de façon "statistique", conformément à l'invention, le filtre d'erreur de prédiction a lorsque le signal prédictible yo à isoler n'est pas connu ou identifiable a priori parmi les données sismiques.
On va exposer dans ce qui suit différentes manières d' appliquer, conformément à l'invention, le filtre d'erreur de prédiction auto-déconvolué

wo 96/09562 PCT/FR95/01231 aux données sismiques.
APPLICATION AUX DONNEES SISMIQUES DU FILTRE D'F_RREUR DE PREDICTION
AUTO-DECONVOLUE
Processus de filtrage fréauentiel ~ Une première façon d'appliquer le filtre d'erreur de prédiction auto-déconvolué aux données sismiques consiste à appliquer directement aux données sismiques, dans le domaine fréquentiel, le filtre M(f) tel que défini précédemment par l'équation (9').
Processus de filtrage temnorel non récursif ~ Une deuxième façon d'appliquer le filtre d'erreur de prédiction auto-déconvolué consiste à effectuer un filtrage temporel non récursif.
L'autocorrélation précolorisée R(f) - IA(f)I= + ~'- IB(f)l'- dans le domaine fréquentiel a pour équivalent dans le domaine temporel, en écriture symbolique r(i) = a(i) * a(-i) + E2 b(i) * b(-i) (10) De préférence, pour ne pas alourdir le temps de calcul, la longueur du filtre de précolorisation b(i) est choisie inférieure à la longueur de a, de sorte que l'ensemble des indices pour lesquels r(i) est non nul est limité à -p, ..., 0, ..., p.
Après avoir calculé l'autocorrélation précolorisée r(i), on calcule l'équivalent g(i) dans le domaine temporel de G(f) qui vérifie R(f) = 1/
LG(f)l'-.
par résolution d'un système de Toeplitz formé à partir de l'autocorrélation précolorisée r(i), à l'aide d'un algorithme de Levinson par exemple . on calcule la matrice M présentant une forme dite de Toeplitz, définie par ses coefficients m(i,j), égaux à
m(i,j) = r(i j) pour i,j = 0, ..., Ls 3 0 Lg vérifiant E m(i,j) g(j) = 8(i) j=0 avec 8(0) = 1 et 8(i) = 0 pour i = 1, ..., Lg et L~ choisi pour obtenir la convergence souhaitée de 1 / IG(f)12 vers R(f).
On calcule ensuite e(n) = a(-i) * g(-i) . g(i) * a(i) * y(n) (11) ~. wo ~io~2 2 i 7 7 i $ 0 On peut résumer le processus de filtrage temporel non récursif de la façon suivante (processus N° 11 1 ) Calcul de r(i), autocorrélation précolorisée de a(i) 5 r(i) = a(i) * a(-i) + E2 b(i) * b(-i) 2 ) Calcul de g(i) à partir de r(i), g(i) étant tel que r(i) * g(i) * g(-i) =
8(i), avec 8(i) = 1 si i = 0 et 8(i) = 0 sinon 3 ) Calcul de e(n) = a(-i) * g(-i) . g(i) * a(i) * y(n) On calcule ensuite le signal à isoler yo(n) = y(n) - e(n).
10 g(i) a une longueur L~ +1 qui peut être grande si E est petit. Dans ce cas, l'application de g(i) aux données sismiques est coûteuse en temps de calcul.
Processus de filtrag-e temporel récursif Pour remédier à cela, on peut factoriser l'autocorrelation 15 précolorisée R(f), qui est de longueur finie, en écrivant l'autocorrélation r(i) de longueur temporelle p comme l'autocorrélation d'un signal c(i) de longueur p et à phase minimale.
~ Une troisième façon d'appliquer le filtre d'erreur de prédiction auto-déconvolué consiste alors à effectuer un filtrage temporel récursif, qui présente l'avantage, comme expliqué plus haut, de nécessiter un temps de calcul moindre que le filtrage temporel non récursif.
On cherche le filtre, de coefficients c(i) (i = 0. .. , p) et appelé pour la suite "filtre d'erreur de prédiction amorti", tel que, en écriture symbolique on ait r(i) = a(i) * a(-i) + ~'- b(i) * b(-i) = c(i) * c(-i) ( 12) qui est la traduction dans le domaine temporel de l'équation (6') dans le domaine fréquentiel R(f) = IA(f)l'- + E2 IB(f)l'- = IC(f)l'-On détermine de préférence les coefficients c(i) par le processus itératif de filtrage temporel récursif suivant (rrocessus N° 2) (on pourra se reporter à l'article de Kunetz, G. "Généralisation des opérateurs d'anti-résonance à un nombre quelconque de réflecteurs" Geophysical Prospectory, vol. 12, pp. 283-289) 1 ) initialisation de variables de calcul en et em ep(0, i) = r(i) pour i = 0, ..., p WO 96109562 - PC"f/FR95/01231 2°. î7180 ep(0, i) = r(i) pour i = 0, ..., p em(0, i) = r(i+1 ) pour i = 0, ..., p-1 em(0, p) = 0 2) itérations Pour j allant de 0 à L, où L est le nombre d'itérations que l'on se fixe K(j) _ - em(j, 0) / ep(j, 0) ep(j+1, i) = ep(j, i) - x(j) em(j, i) pour i = 0, ..., p em(j+1, i) = em(j, i+1) - tc(j) ep(j, i+1) pour i = 0, ..., p-1 em(j+1, p) = 0 Fin de boucle en j.
On remarquera que les tc(j) sont calculés de sorte que em (j+1, -1) = 0 3) normalisation On calcule ensuite c(i) par c(i) = ep(L, i) / J ep(L, 0) pour i = 0, ..., p On choisit L pour avoir la précision recherchée dans la relation asymptotique r(i) = c(i) * c(-i).
On calcule, connaissant c(i) et a(i), u(n), pour n variant de 0 à N, par P P
u(n) _ [~ a(i) y(n-i) - E c(i) u(n-i)] / c(0) (13) i~ i=1 puis, connaissant u(n), on calcule pour n variant de N à 0, e(n) par P P
e(n) _ [~ a(i) u(n+i) - ~ c(i) e(n+i)] / c(0) (1.~) i~ i=1 u(n) s'obtient ci-dessus par un filtrage dit "récursif ascendant"
que l'on écrira par la suite en écriture symbolique a(i) u(n) _ [_____] * y(n) c(i) e(n) s'obtient ci-dessus par un filtrage dit "récursif descendant"
que l'on écrira par la suite en écriture symbolique v.. wo mo9s62 21 7 l 1 8 0 . a(-i) e(n) _ [---] * u(n) c(-i) S
On utilisera aussi par la suite la notation y(n) u(n) _ [____-]
c(i) pour désigner l'opération b(i) u(n) _ [____] * y(n) c(i) L'homme du métier notera que c(i) est à phase minimale et c(-i) à phase maximale, ce qui permet d'obtenir des filtrages récursifs respectivement ascendant et descendant stables.
On peut résumer le processus de filtrage temporel non récursif de la façon suivante (processus N° 3 ) 1) Calcul de r(i), autocorrélation précolorisée de a(i), par r(i) = a(i) * a(-i) + sz b(i) * b(-i) 2) Calcul du filtre d'erreur de prédiction amorti c(i) (i = 0, ..., p) vérifiant c(i) * c(-i) p r(i) 2~ 3) Calcul de u(n) par filtrage récursif ascendant P P
u(n) _ [~ a(i) y(n-i) - ~ c(i) u(n-i)] / c(0) i=0 i=1 4) Calcul de e(n) par filtrage récursif descendant P P
e(n) _ [~ a(i) u(n+i) - ~ c(i) e(n+i)] / c(0) i=0 i=1 . 35 Le signal à isoler yo(n) est ensuite calculé par yo(n) = y(n) - e(n).
Processus de filtrage récursif non stationnaire ~ Une quatrième façon pour appliquer le filtre d'erreur de prédiction auto-déconvolué, conformément à l'invention, consiste à
effectuer un filtrage que nous appellerons "récursif non stationnaire", qui WU 96109562 PC"TIFR95101Z31 ~17718~

présente l'avantage de réduire les effets de bord.
La quasi-prédictibilité de yo peut s'écrire, sous forme matricielle Ae+~Bp=Ay, (15) où A et B sont des matrices construites à partir des filtres a et b et p est un vecteur de calcul dit "vecteur d'innovation", précolorisé par B et multiplié
par le facteur de précolorisation ~ . Si ~ est nul, yo est parfaitement prédictible.
En explicitant l'équation matricielle (IS) on a a(p)...a(0)0 0...0 e(0) b(q)...b(0) 0 ... 0 p(0) a(p)...a(0) 0 0 y(0 0 a(p)...a(0) 0...0 e(1) 0 b(q)...b(0) 0 p(1) y(1 ... +~ ... ... .. ...
0 0 O...a(p)...a(0) 0 0 b(q)...b(0) 0 0 a(p)...a(0) le(N) p(N) y(N
On cherche e et p minimisant =e*e+p*p où * désigne le transposé conjugué d'une matrice ou d'.un vecteur.
est minimum pour e pris égal à emin = A' (A A* + ~'- B B*)-1 Ay e t p pris égal à purin = ~ B* (A A* + ~'- B B*)-1 Ay vaut alors min = emin* emin + pmin* purin soit, en explicitant les valeurs de emin et purin ~ min= Y* A* (A A* + ~2 B B*)-t Ay Posons s = Ay qui est le résultat de la convolution de y par A
et v=R-1 s, avec R=AA*+~'-BB*
Pour trouver v tel ~ que Rv = s, l'homme du métier remarquera qu'il s'agit de la résolution d'un système linéaire, R étant une matrice de Toeplitz complexe p-bande formée à partir de l'autocorrélation de a(i) plus ~'- fois wo mo9s62 217 7 î 8 0 celle de b(i). Une façon pour trouver v consiste à adopter la décomposition de Cholesky pour R
R = C C* (C étant une matrice triangulaire p-bande inférieure).
Mais ceci est une opération matricielle sur les données. En effet, la dimension de R est N+1-p, donc de l'ordre de N qui est la dimension des données. Mais une des caractéristiques de l'invention est de constater que l'addition de la précolorisation permet d'approximer une opération matricielle sur les données par un filtrage non-stationnaire pour les premiers points d'applications. Cette structure de filtrage non-stationnaire permet de concilier le faible coût du filtrage et l'absence d'effets de bords.
L'implémentation de cette structure de filtrage non stationnaire dans le cadre de la décomposition de Choleski consiste à arrêter la décomposition au bout du calcul des L premières colonnes de la matrice C, celle-ci étant fictivement complétée par c;+j , j = ci+t. . L pour j > L et i=O,..p.
Nous appellerons décomposition de Choleski partielle ce procédé.
Ce nombre L d'étapes au bout duquel le calcul peut être arrêté avec une erreur acceptable dépend du facteur de précolorisation et en aucun cas du nombre d'échantillons de données N. L peut être choisi d'autant plus petit que le facteur de précolorisation est grand.
On peut noter à ce propos qu'il est l'usage en calcul numérique de rajouter un petit facteur sur la diagonale d'une matrice positive de l'ordre de la précision numérique du calculateur employé (10-6 en _énéral) pour garantir la stabilité numérique de la décomposition de Choleski.
l'autocorrélation de précolorisation utilisée dans la méthode proposée est d'un ordre de grandeur différent, de l'ordre de 10-'- à 10-t, et vise à rendre possible l'utilisation de la structure de filtrage non-stationnaire pour approximer l'inverse de la matrice à l'aide d'un petit nombre d'ensembles de coefficients.
En variante et préférentiellement, on ne calcule pas les coefficients du filtre non-stationnaire par la décomposition de Choleski partielle, mais on calcule des coefficients ep(i, j) par une méthode que l'on appelera par la suite calcul du filtre d'erreur de prédiction amorti non stationnaire (processus N° 5), les coefficients ep(i,j) étant reliés aux c;,j par la relation ci,j=ep(j,i-j)/~ep(j,0), et les coefficients em(i.j) étant des variables internes du processus.

WO 96/09562 PG"T/FR95/01231 Initialisation ep(0, i) = r(i) pour i = 0. ..., p em(0, i) = r(i+1) pour i = 0, ..., p-1 em(0, p) = 0 5 Itérations Pour j allant de 0 à L, où L inférieur ou égal à N est le nombre d'itérations que l'on se fixe pour obtenir la convergence des ep(j, i) quel que soit i K(j) _ _ emV, 0) / epU~ 0) 10 ep(j+1, i) = ep(j, i) - K(j) em(j, i) pour i = 0, ..., p em(j+1, i) = em(j, i+1) - K(j) ep(j, i+1)pour i = 0. .. , p-1 em(j+1, p) = 0 Fin de boucle en j On remarquera que les K(j) sont calculés de sorte que 15 em (j+1, -1) = 0 Normalisation Pour j allant de 0 à L, = I/ '~ ep(j, 0) ep(j, i) _ ~ . ep(j, i) i = 0, .... p 20 Fin de boucle en j.
On complète ensuite fictivement les ep par ep(j.i) - ep(L,i) pour j=L+ 1, .., N et i=0,..,p.
Par compléter fictivement, nous entendons que, lorsque par la suite la notation ep(j,i) apparaîtra, il faudra la remplacer par ep(L.i) si j est supérieur à L. I1 est bien entendu que nous n'avons pas, lors de l'implantation du procédé sur un calculateur, à prévoir un espace mémoire pour stocker les ep(j,i) pour j>L, ni d'ailleurs les c;,j fictivement complété
dans le processus précédent. Ceci montre bien que les Processus N° 4 et N° S
que nous venons de décrire sont bien des filtrages non stationnaires et non des calculs matriciels sur les données.
Le processus N°2 de calcul du prédicteur amorti c(i) revient à
prendre:
c(i)=ep(L,i) En conséquence, ep(j,i) peut être considéré comme un filtre d'erreur de prédiction amorti non-stationnaire.

2171180v Le filtrage récursif non stationnaire préféré comporte en résumé
les étapes suivantes (rzrocessus N° 6) 1 ) Calcul de l' autocorrélation précolorisée r(i) = a(i) * a(-i) + E2 b(i) * b(-i) 2 ) Calcul d'après le processus N° 5 à partir de r(i) du filtre d'erreur de prédiction amorti non-stationnaire ep(0, ..., L ; 0, ..., p) complété fictivement à ep(0, ..., N ; 0. ..., p) 3 ) Calcul de P
s(n-p) = E a(i) y(n-i) n = p, ..., N
i=0 4 ) Calcul de i=min(p,n) u(n) _ (s(n) - E ep(n-i, i) u(n-i)~ / ep(n, 0) i=1 pour n=0,...,N-p 5 ) Calcul de i=min(p,N-p-n) v(n) _ (u(n) - E ep(n, i) v(n+i)~ / ep(n, 0) i=0 pour n=N-p,...,0 6 ) Calcul de i=min(p,N-p-n) e(n+p) = E . a(i) v(n+i) i=max(0,-n) p o u r n = -p, .. , N-p Dans la suite, on utilisera symboliquement la même écriture u = C-~s pour désigner le calcul de u à partir de s par une méthode filtrage non-stationnaire, indépendamment de la méthode utilisée, que ce soit la décomposition de Choleski partielle (processus N° 4) ou la méthode préférentielle du processus N° 5 et de l'étape 4) du processus N° 6 ci-dessus.
ou tout autre méthode comprenant l'étape de choix d'un nombre L
d'itérations au-delà de laquelle l'approximation stationnaire est faite. En particulier, le procédé de filtrage non récursif (procédé N° 1) peut être lui aussi généralisé en un procédé de filtrage non récursif non stationnaire.En effet, l'algorithme de Levinson appliqué pour un nombre d'itérations égal au nombre d'échantillons des données permet de calculer une matrice G telle que : R-~ = G*G, ce qui permet d'écrire 21~7~80 e = A*G* GAy (4 ) Ceci est une opération matricielle sur les données. L'enseignement de l'invention permet d'approximer cette opération par un filtrage non stationnaire en arrêtant l'algorithme de Levinson au bout de L itérations et en complétant fictivement la matrice G.
L'écriture calcul du filtre récursif non stationnaire C tel que OC' = R et application de ce filtre par u = C- t s pourra donc aussi être comprise comme le calcul du filtre non récursif non stationnaire G tel que G'G ~ R-t et application de ce filtre par u = Gs.
De même on utilisera symboliquement l'écriture v = C*-1 u, pour désigner le calcul de v à partir de u, indépendamment de la méthode utilisée (par exemple étape 5) du processus N° 4, mais de façon préférentielle par l'étape 5) du processus N° 6).
Dans le cas où le filtre d'erreur de prédiction a n'est pas calculé de façon déterministe, à partir du spectre Yo(f) du signal yo, on calcule a par itérations, de façon dite "statistique".
CALCUL STATISTIQUE DU FILTRE D'ERREUR DE PREDICT10N
Tout comme l'application du filtre d'erreur de prédiction auto déconvolué, le calcul statistique du filtre d'erreur de prédiction n peut s'effectuer de différentes manières. Dans tous les cas, la partie non prédictible des données se calcule à l'aide d'un opérateur linéaire M : e = My M peut prendre différentes expressions - en fréquentiel M(f) est défini par l'équation (S') ;
- en temporel, filtrage non récursif : m(i) = a(-i) * g(-i) * g(i) * a(i) filtrage récursif stationnaire a(-i) a(i) m(i) _ (_______j * (_____ c(-i) c(i) filtrage récursif non stationnaire : M = A* C*-1 C-1 A.
filtrage non récursif non stationnaire : M = A* G' GA
Filtrage récursif non stationnaire a peut être calculé, conformément à l'invention, par minimisation d'une norme d'erreur dépendant de a et de y ; ainsi, il est possible pour calculer a de minimiser la norme de la partie non prédictible e des données wo 96ro9562 217 718 0 p~~~g/01231 sismiques : Ilell2 = e* e.
e s'exprime à partir des données par e=My M étant un opérateur linéaire dont l'expression dépend de la méthode de S déconvolution utilisée. Cet opérateur dépend du filtre d'erreur de prédiction.
Ona:
Ilell2 = y* M* M Y
Il est toutefois préférable de chercher à minimiser en a(0), ..., a(p) min = e* e + p * p = y* I~Iy ( 15' ) soit compte tenu de ce qui précède, de chercher à minimiser IIC-lAyil'-.
On va maintenant décrire plusieurs façons de calculer le filtre d'erreur de prédiction minimisant IIC-lAYll'-~ conformément à l'invention.
Pour calculer le filtre d'erreur de prédiction minimisant IIC-lAyll'-une première façon consiste à poser d = C-1 Ay, ( 16) qui s'écrit encore d = C-1 Y~
où ~ est un vecteur dont les composantes sont a(0), ..., a(p) et Y la matrice de Toeplitz formée à partir de y, à p+1 colonnes yo. .. . YP
I y(p-i) I i = 0, ..., p avec y; = I ... I
I y(N-i) I.
Minimiser IIC-tAy112 revient à minimiser d* d.
Calculons pour ce faire la matrice U = C-tY, puis U* U, pour résoudre finalement, en négligeant la dépendance de U en (U* U)h = (1, 0, ..., 0)*.
Aprs avoir calcul h, de ..., on calcule composantes h(p) h(0), ensuite a(i) en normalisant de tellefaon que a(0)=1par h(i) a(i) = h(i) / h(0) i = , ..., p.

De faon prendre en compte la dpendance U en on reprend de a_, ensuite les calculs prcdentsavec nouvelles valeursde a(i)jusqu' les la convergence souhaite.

Le processus dcrit plus pour calculer de faonstatistique haut n peut se rsumer de la faon cessus N 7) suivante (pro 1 ) Choix de p+1 coefficients a(i) initiaux wo 9~s62 21 7 7 i ô ~

2) Calcul du filtre non stationnaire C à partir de l'autocorrélation procolorisée r(i), C étant tel que CC* = R, R étant la matrice de Toeplitz formée à partir de r(i) 3 ) Calcul par filtrage récursif non-stationnaire des p+1 vecteurs ui : ui = C-tYi 4 ) Calcul de la matrice de covariance des ui : COV"(i, j) = ui~ u~
5 ) Résolution de P
COV"(i, j) h(j) = 8(i) i = 0, ..., p avec b(1) = 0 et 8(i) = 0 pour i = 1, ..., p j=0 6 ) Calcul de a(i) = h(i) / h(0) i = 0, ..., p 7 ) Reprendre à l'étape 2) jusqu'à la convergence souhaitée des valeurs a(i) 8 ) Fin.
~ En variante, pour calculer le filtre d'erreur de prédiction, on peut effectuer le filtrage non-récursif non stationnaire en remplaçant les étapes 2) et 3) du processus N° 7 par 2) Calcul du filtre non stationnaire G à partir de l'autocorrélation précolorisée r(i), G étant tel que G*G $ R-1, R étant la matrice de Toeplitz formée à partir de r(i) 3 ) Calcul par filtrage non-récursif non-stationnaire des p+1 vecteurs u; tels que ui = Gy;
~ En variante, pour calculer le filtre d'erreur de prédiction minimisant IIC-lAYll'-. on peut effectuer le filtrage dit récursif stationnaire.
L'équation précédente d - C-tAy s'écrit, en écriture symbolique dans le domaine temporel d(n) = a(i) * u(n), (17) y(n) avec u(n) _ [-------J (18) c(i) Après avoir calculé le filtre d'erreur de prédiction amorti c(i), on calcule u(n) par l'équation (18) puis on cherche les valeurs a(i) minimisant Id(n)l'-.
I1 existe de nombreuses manières connues de l'homme du métier pour trouver les a(i) minimisant E Id(n12.

w0 2 PCT/FR95/01231 Si on utilise une méthode dite de corrélation, on considère que u(n) - est connu pour n = 0. ..., N et est nul en dehors de ces valeurs de n, donc d(n) est non nul pour n = 0, ..., N+p et nul ailleurs. On minimise donc en fait 5 N+p E Id(n)l'-n=0 De prfrence, on utilise une mthode dite de covariance qui a 10 l' avantage de mieux prendre en compte les effets de bord.
On considre que u(n) est connu pour n = 0, ..., N et est inconnu ailleurs, donc d(n) n'est connu que pour n = p, ..., N et on minimise N

E Id(n)I? .

15 n=p Le processus de calcul statistique du filtre d'erreur de prdiction par filtrage rcursif stationnaire peut se rsumer de la faon suivante (pro cessus N 8) 20 1 Choix d'un filtre d'erreur de prdiction initial n, de ) coefficients a(0), a( 1 ), .. . , a(p) 2) Calcul de l'autocorrlation prcolorise r(i) = a(i) * a(-i) + '- b(i) * b(-i) 3 Calcul du filtre d'erreur de prdiction amorti c(i) partir ) de r(i) par le processus N 2.

25 4 Calcul par filtrage rcursif de u(0), ..., u(N) par ) y(n) u(n) _ [____-___]

c(i) 5 Calcul de la matrice de covariance COV" de u(n) ) n = N - min(i, j) COV" (i, j)= ~ u(n-i)u(n-j) n = max(i, j) pour i=0....,p et j=0,...,p 6) Calcul du vecteur h de composantes h(0), ..., h(p) vrifiant COV"h = [1. 0, ..., 0]* (par la mthode de Cholesky par exemple) 7 Calcul de a(i) par normalisation de h(i) ) a(i) = h(i) l h(0) 8) Reprise l'tape 2) jusqu' la convergence souhaite des valeurs a(i).

En variante, pour calculer le filtre d'erreur de prdiction minimisant IIC-1 A y I l'-, on peut effectuer un filtrage non récursif, qui comporte les étapes suivantes (processus N° 9) : ' 1 ) Choix d'un filtre d'erreur de prédiction initial a, de coefficients a(0) =
1.
a(1), ..., a(p) ' 2 ) Calcul de l'autocorrélation précolorisée r(i) = a(i) * a(-i) + ~'- b(i) *
b(-i) 3 ) Calcul de g(i) à partir de r(i) par le processus N° 1 (étape 2) 4 ) Calcul de. u(n) = g(i) * y(n) 5 ) Calcul de la matrice de covariance COV" de u(n) n = N - min(i, j) COV" (i, j) = E u * (n-i) u(n-j) n = max(i, j) avec i=0....,p et j=0....,p 6 ) Calcul de h(0), ..., h(p) vérifiant COV"h = [1, 0, ..., 0]* (par la méthode de Cholesl:y par exemple) 7 ) Calcul de a(i) par normalisation de h(i) a(i) = h(i) / h(0) 8 ) Reprise à l'étape 2) jusqu'à convergence souhaitée des valeurs a(i).
~ Dans une troisième variante, pour calculer le filtre d'erreur de prédiction minimisant IIC-lAyll'-, en oeuvre les tapes suivantes (processus N 10) on met 1 ) Choix d'un re d'erreur de prdiction initial a, dc coefficients filt a(0) = 1, a(1), ..., a(p) 2 ) Calcul de par transforme de Fourier de a(i) A(f) 3 Calcul de = IY(f)l'- ! [IA(f)l'- + >r'- IB(f)l'-]
) UU(f) 4 ) Calcul de par transforme de Fourier inverse de UU(f) uu(i) 5 ) Calcul de en rsolvant le systme d'quations de Toeplitz form a(i) avec uu(i) par l'algorithme de Levinson sans second membre par exempte 6) Reprise l'tape 2) jusqu' la convergence souhaite des valeurs a(i).

Calcul st atistique du filtre de ",nrcolorisation On peut, sans sortir du cadre de l'invention, utiliser des techniques de minimisationfonctions, connues de l'homme du mtier, pour minimiser de une norme donne en fonction des coefficients b(0), .. , b(q) du filtre de prcolorisation en fixant ou non la valeur . On peut en particulier b, dans un calcul itératif estimer ~ et/ou b en fonction des statistiques du bruit wo Wo9s6Z 217 718 Q PCT/FR95~01231 estimé e et celles du vecteur d'innovation p à l'itération précédente.
L'expression du vecteur d'innovation p, p=~ B* R-1 Ay, montre qu'il peut être calculé en tant que sous produit du calcul de e.
EXEMPLE-S D'APPLICATION D'UN PROCEDE SELnN L'INVI=NTION DANS LE PLAN
L'invention permet, lorsque le filtrage récursif est utilisé.
d'effectuer un filtrage déterministe dans le plan f-x, quelle que soit la forme du signal yo à isoler dans le plan f-k, sans que la sélectivité du filtre ne se traduise par un temps de calcul élevé, et, dans le cas où le filtrage récursif non stationnaire décrit plus haut est employé, l'invention permet d'éviter les inconvénients dûs aux effets de bord rencontrés dans les procédés de l'art antérieur.
~ Un exemple préféré de mise en oeuvre du procédé selon l'invention est le suivant, lorsque le signal à isoler est défini par sa localisation dans le plan f k (filtrage déterministe) 1 ) Détermination par l'utilisateur d'une zo n e du plan f-k contenant le signal yo à isoler, 2 ) Pour toute valeur de la variable espace x, effectuer la transformée de Fourier des données y(t, x) de la variable temporelle t pour obtenir des donnes transformes Y(f, x) dans le plan f-x, 3 ) Pour chaque frquence f du plan f-x, 3a) prendre y(n) = Y(f,n), 3b) dterminer les intervalles pour la variable k, le signal contenant prdictible yo que l'on cherche isoler, en effectuantcoupe une f constant de ladite zone, 3c) calculer le filtre d'erreur de prdiction spatial lmentaire (variable x) correspondant chacun de ces intervalles en k donn (en pour un E

utilisant par exemple la mthode de polynmes de Chebychev dcrite plus haut), 3d) calculer le filtre d'erreur de prdiction global et l'auto-a(i) corrlation prcolorise correspondante r(i), 3e) calculer le filtre d'erreur de prdiction non-stationnaireep(j,i) partir de l' autocorrlation prcolorise r(i) (processusN 5 de prfrence), 3f) calculer la partie non prdictible e(n) des donnestapes par les 3), WO 96109562 PCT/FIt95/01231 21 ?7180 4), 5), 6) du processus de filtrage récursif non stationnaire (processus N° 6), 3g) soustraire la partie non prédictible e(n) aux données y(n) pour obtenir la partie prédictible yo(n), 3h) Reprendre à l'étape 3a) pour la fréquence f suivante sauf si la dernière fréquence du plan f-x est atteinte, 4 ) Pour tout x, effectuer la transformée de Fourier inverse de la variable f de e(n) ou de yo(n) pour revenir dans le plan t-x.
~ L'application de l'invention aux données sismiques y(f, x) dans le plan f-x et dans le cas où le filtre d'erreur de prédiction a est calculé de manière statistique, permet de pallier aux inconvénients du filtrage prédictif rencontrés dans l'art antérieur, notamment tel que proposé par Cavales.
Un exemple préféré de mise en oeuvre du procédé selon l'invention, dans le cas d'un filtrage statistiaue, est le suivant 1 ) Pour tout x, transformée de Fourier des données sismiques y(t, x) pour la variable t pour se placer dans le plan f-x 2 ) Pour chaque fréquence f du plan f-x 2a) prendre y(n) = Y(f, n) 2b) choix d'une longueur p pour le filtre d'erreur de prédiction, d'un facteur ~ et d'un filtre de précoloriement b(0). ..., b(q), ainsi que d'un nombre d'itérations pour la convergence de a(i) 2c) calcul par filtrage récursif ascendant de y(n) u(n) _ (_______~ .
c(i) c(i) étant le filtre d'erreur de prédiction amorti de l'itération précédente (celui de la dernière itération de la fréquence précédente si c'est la première itération pour la fréquence en cours d'itération) 2d) calcul de la matrice de covariance de u(n) (étape 5) du processus N° 8) ; en variante, calcul de la matrice de corrélation de u(n) 2e) résolution du système pour obtenir les coefficients a(i) du filtre d'erreur de prédiction (étapes 6) et 7) du processus N° 8 de préférence), 2fj calcul de l'autocorrélation précolorisée r(i) (étape I) du processus N° 6) et du filtre d'erreur de prédiction amorti c(i) (par processus N° 2 de préférence) 2g) retour en 2c) jusqu'à la convergence souhaitée des valeurs de a(i) 2h) obtention de la partie non prédictible e(n) des données par le processus N° 6 ,. WO 96/09562 217 718 0 PCT/FR95I01231 2i) soustraction de la partie non prédictible e(n) aux données y(n) pour obtenir la partie prédictible yo(n) 2j) retour en 2a) sauf si dernière fréquence 3 ) Pour tout x, transformée de Fourier inverse de yo(n) de la variable t pour obtenir yo(t, x).
Le filtre d'erreur de prédiction est avantageusement calculé d'une façon récursive stationnaire, l'application du filtre d'erreur de prédiction autodéconvolué s'effectuant de façon récursive non stationnaire.
La figure 4 représente une section sismique t-x synthétique, l'axe de temps étant vertical et l'axe des x horizontal. Les données sismiques y(t, x) contiennent trois évènements prédictibles (matérialisés par des droites) que l'on cherche à isoler du bruit aléatoire non prédictible.
La figure 5 représente la partie prédictible yo(t, x) après mise en oeuvre du procédé selon l'invention. On remarquera la forte atténuation du I 5 bruit aléatoire.
La figure 6 représente la partie non prédictible e(t, x) et l'on peut remarquer l'absence de signal prédictible même sur les bords (pas d'effets de bord).
Les figures 7 et 8 illustrent respectivement les parties prédictibles et non prédictibles obtenues par la mise en oeuvre du procédé décrit par Canales, amélioré en prenant un filtre d'erreur de prédiction bilatéral a(-p).
.. , a(0) = 1, ..., a(p) et en considérant que la partie non prédictible s'obtient en appliquant aux données ce filtre d'erreur de prédiction bilatéral.
On remarquera sur la figure 7 la moins bonne atténuation du bruit aléatoire et sur la figure 8 un résidu de signal prédictible.
Ainsi, appliquée à l'atténuation du bruit aléatoire, l'invention réalise une meilleure atténuation du bruit et une meilleure préservation du signal que celle permise par l'état de la technique.
L'invention n'est bien entendu pas limitée aux exemples qui viennent d'être décrits.
On peut, par exemple, imposer une forme particulière au filtre d'erreur de prédiction a, en choisissant les pn,in premiers coefficients du filtre d'erreur de prédiction nuls. Le filtre d'erreur de prédiction est alors de la forme (1, 0, ..., 0, a(pmin + 1), ..., a(p)]. Il faut bien entendu tenir compte de cette forme particulière du filtre d'erreur de prédiction lors de la résolution du système d'équations formé à partir de la matrice de covariance COV" dans les processus de calcul statistique du filtre d'erreur de prédiction. La forme particulière du filtre peut avantageusement servir à modéliser le bruit prédictible dû aux réflexions multiples, en vue de son élimination.
Les procédés conformes à l'invention, précédemment décrits, 5 peuvent être mis oeuvre dans le domaine temporel ou spatial, et se généraliser aisément à plusieurs dimensions. On peut facilement traduire les équations à une dimension en équations à deux dimensions, etc..., à l'aide d'une notation appropriée.
Ainsi, par exemple, la prédictabilité de yo(nxt, nx~) de support nxl =0, ..., 10 Nxt, nx2 = 0, ..., Nx2 par le filtre d'erreur de prédiction a(i, j) de support i = 0, .. .
Pxt, j = 0. .~.. Px2 Peut s'écrire a(i. .l) ** Yo(nxt, nx2) = 0 (19) où * * est l'écriture symbolique d'un produit de convolution à deux dimensions.
15 Soit Yonxt le vecteur de composantes que Yonxt (nx2) = Yo(nxt, nx2) et Ai la matrice de composantes Ai (n, j) = a(i, n j).
L'équation (19) s'écrit P
Ai Yonxl-i = 0 20 i=0 c'est-à-dire A; (*) Yonx1 = 0 où (*) est l'écriture symbolique pour le produit de convolution matriciel.
Avec cette notation l'équation à une dimension 25 r(i) = a(i) * a(-i) + ~2 b(i) * b(-i) s'écrit Rn = An A*n + E2 Bn B*n où pour tout n, Rn est une matrice de Toeplitz par blocs.
Dans le cas du filtrage non récursif, la matrice du système à
30 résoudre pour trouver les matrices Gn, a la structure d'une matrice de Toeplitz par blocs, le système étant donc résoluble par un algorithme de Levinson par blocs. Accessoirement, chacun des blocs est lui-même une matrice de Toeplitz.
Dans le cas du filtrage récursif, le calcul des matrices Cn ou EPji à
partir de Rn s'effectue par analogie avec les processus N° 2 et 5 de filtrage récursif précédemment décrits. Dans le cas déterministe, la méthode par wo 9b/o9s62 217 718 0 moindres carrés se généralise aisément à plusieurs dimensions.
L'invention peut encore s'appliquer dans le domaine f-kxt-xt en présence de données tridimensionnelles fonction des variables t, xt et x~. On utilise avantageusement, dans le domaine f-x t - x ~ un filtre d'erreur de prédiction à deux dimensions tel que décrit plus haut, et dans le domaine f-kx2-xt on effectue pour chaque kx~ un traitement analogue à celui des données sismiques dans le plan f-x, tel que décrit plus haut.
L'invention peut encore s'appliquer pour effectuer un traitement dit "multicanal", pour lequel on dispose de m jeux de données yj Yj = ~Yj(0). ..., yj(~1 j = 1, ..., m comprenant une partie prédictible yoj pour un filtre d'erreur de prédiction commun.
On a yj(n) = yoj(n) + ej(n) et a(i) * yoj(n) ~ 0.
On effectue conformément à l'invention un calcul statistique multicanal du filtre d'erreur de prédiction en minimisant m Yj * MYj.
j=1 ' par calcul des matrices de covariance COVj correspondant à chaque jeu de données yj, en calculant EjCOVj = COV et en résolvant le système avec COV.
Dans le processus N° 8 cela correspond aux étapes 4) et 5) qu'il faut faire pour tout j, l'étape 6) se faisant avec la matrice COV.
L'invention s'applique avantageusement à la restauration de données manguantes.
I1 arrive souvent, dans la pratique, qu'il manque dans la section sismique t-x la trace correspondant à un capteur placé à une abscisse x sur la grille d'acquisition, ou bien que la trace correspondant à ce capteur soit inutilisable à cause d'un défaut de l'installation d'acquisition des données ou encore parce que le signal délivré par ce capteur a été saturé pendant l'enregistrement par un bruit très fort.
L'invention permet de restaurer la composante prédictible des données manquantes.
Les données sismiques sont alors supposées contenir un signal prédictible pour un filtre d'erreur de prédiction de coefficients a(i), un bruit impulsif pour les valeurs manquantes et un signal additif non prédictible.
On écrit les données sous la forme m y = z - ~ w(j)pj j=1 où z désigne le vecteur contenant les données sismiques existantes et connues, pj le vecteur contenant la réponse impulsionnelle du jième bruit impulsif correspondant à la ji~me valeur manquante d'indice I(j).
Le vecteur pj a des composantes pj(0), ..., pj(N) telles que p~ (n)=1 si n=I(j) et pj(n)=0 si n ~ I(j) w(j) est un scalaire désignant l'amplitude inconnue du ji~me bruit impulsif.
Les données manquantes dans le vecteur z peuvent être prises quelconques, par exemple égales à 0.
Dans ce cas, .
z(I(j)) = 0 et w(j) _ -y(I(j)).
La valeur reconstituée pour l'emplacement I(j) est donc -w(j).
Soit P la matrice dont les colonnes sont les vecteur p~.
Soit w le vecteur dont les composantes sont les scalaires w(j).
On a y=z-Pw On doit effectuer sur ces données reconstituées y la séparation entre la partie prédictible yo et non prédictible e.
Y=Yo+e Ona:
z=yo+e+Pw e se calcule à partir de y à l'aide d'un opérateur linéaire M, dont l'expression varie en fonction du mode de déconvolution utilisé, mais qui se met dans tous les cas sous la forme M=H*Hete=My Dans le cas récursif non stationnaire, on a par exemple : H = C-1 A.
On détermine le vecteur w en minimisant min + ä,2 ~'~'* w = y* My + iv,2 w* w d'après l'équation (15'), ~, étant un facteur de préblanchiment.
On cherche à minimiser la quantité
(z - Pw)* M (z - Pw) + ~.2 w* w (20) On trouve w minimisant (20) en résolvant (P* M P + ~,2 I) w = p* M z, I étant la matrice identité.

wo mo~2 rc~r~siom En posant D = HP et s = Hz on a (D* D + i~,2 I)w = D* s On en déduit l'expression des données reconstituées y . y = z _ p(D* D ~ ~,2 I)-1 D* s L'expression de la partie non prédictible est e = H* H [z - P(D* D + 7i.2 I)-1 D* s]
soit encore e = H* [s - D (D* D + 7~2 I)-I D* s]
et la partie prédictible yo se calcule par Yo=Y-e yo = z - P(D* D + ~,2 I)-1 D* s - H* [s - D (D* D + ~,2 I)-1 D* s]
Si l'entrée est z et la sortie les données reconstituées y, on obtient grâce à l'invention une restauration des données manquantes.
Si l'entrée est z et la sortie la partie prédictible yo, on obtient grâce à l'invention une restauration des données manquantes et une atténuation du bruit.
On peut bien entendu, sans sortir du cadre de la présente invention, généraliser au cas où les pj sont constitués par n'importe quelle réponse impulsionnelie d'un bruit de forme pj connue et d'amplitude w(j) inconnue.
Si le filtre d'erreur de prédiction a est inconnu, on a deux vecteurs w et ~ trouver.

On part de valeurs initiales pour (par exempleavec y(n) les valeurs manquantes 0). On calcule a(i) avec constant y(n) par une mthode de calcul statistique ou dterministe de filtre d'erreur de prdiction. Puis le procd prcdent de calcul de y partir de z tre utilis actualiser peut pour y(n) avec a(i) constant. Puis on revient au calculde a(i) danscas d'un le calcul statistique.

Les donnes d'entres tant les valeurs et chaque connues de y(n) vecteur pj pour j = 1, ..., m, la rponse pulsionnelledu j~~me le procd im bruit, mis en oeuvre comporte les étapes suivantes dans lc cas statistique (processus N° 111 1 ) Initialisation des données y(n).
2) Calcul du filtre d'erreur de prédiction a(i) pour le y(n) courant.
3 ) Actualisation de y(n) pour le a(i) courant 3a) calcul de m vecteurs dj = H pj j = 1. ..., m 3b) calcul du vecteur s = H y 3c) calcul de la matrice x(i, j) = d;* dj i, j = 1, ..., m 3d) addition du facteur de préblanchiment x(i, i) E- x(i, i) + 71.2 i =1, ..., m 3e) calcul du vecteur k(i) = d;* s 3f) résolution du système linéaire d'ordre m m r(i, j) w(j) = k(i) i = 1, ..., m j=1 3g) actualisation de y(n) par m Y E- Y - ~ wV) Pj.
j=1 4 ) Retour en 2) ou 5) jusqu'à la convergence souhaitée pour a(i) 5 ) Calcul de u m u=s-~ w(j)dj.
j=1 6 ) . Calcul de e = H* u 7 ) Calcul de yo = y - e.
Le procd ci-dessus selon l'inventionest avantageusement appliqu aux cas suivants Bruit c ohrent En prenant pour colonnes de P les bruit composantes d'un cohrent d'amplitude inconnue, on peut modliserprsence d'un bruit la cohrent dans les donnes sismiques. Par exemple,le plan f-x, bruit dans un de relation de dispersion connue k;(f) pi(n) = ejnexx;(ti avec n = 0, ..., N

Traces manauantes On a reprsent sur la figure 9 une dans section sismique t-x laquelledix traces sont manquantes.

De faon gnrale, supposons qu'il y m traces manquantes ait d'indices I(1)...I(m) sur la section sismique.

On construit la matrice P m colonnesla ji~me colonneest dont le vecteur pj de composantes pj (0), ..., pj (N) telles que pj (n)=1 si n=I(j) et pj(n)=0 si n x I(j).
On modélise ainsi l'existence d'un bruit impulsif pour les valeurs données, ce qui revient à dire qu'elles sont inconnues.
On réduit avantageusement le temps de calcul en ne calculant 5 qu'un seul vecteur d = Hpo avec po(p) = 1 et po(n) = 0 si n x p.
Dans ce cas Apo = [a(0), ..., a(p), 0, ..., O~T.
10 Si on utilise le filtrage récursif stationnaire pour appliquer H, alors a(i) d(i) _ [__ c(i) Ceci est une division puissances croissantes de deux polynômes qui peut être arrêtée dès que les valeurs de d(i) deviennent négligeables.
Le vecteur dj s'obtient en décalant le vecteur d de I(j)-p échantillons vers le bas et en mettant à 0 les premiers I(j)-p échantillons.
dj(n) = 0 pour n < I(j) - p dj(n) = d(n-I(j)+p) pour n > I(j) - p Dans le cas récursif non stationnaire les opérations de la forme s =
Hy se font par les étapes 3) et 4) du processus N° 6, et les opérations de la forme e = H* u par les étapes 5) et 6) du même processus.
Une variante simplifiée consiste à remplacer l'étape 3) d'actualisation de y(n) à a(i) constant par un procédé itératif ne nécessitant pas de résolution de système linéaire. L'étape 3' remplaçant l'étape 3) est 3') Actualisation de y(n) pour le a(i) courant 3a) initialisation de y;t(n) = y(n) 3b) calcul de la partie non prédictible e(n) de y;t(n) 3c) actualisation de y;i(n) pour les seuls échantillons n correspondant aux valeurs manquantes (n = I(j) j = 1, ..., m) par Yic(n) f- (1 - ~2) Yic(n) + a.2 Y(n) - e(n) 3d) retour en 3b) ou 3e) si convergence souhaitée de y;i(n) 3e) actualisation de y(n) pour les seuls échantillons n correspondant aux valeurs manquantes par Y(n) = Yic(n).

wo s6i rc°r~siom On a représenté sur la figure 10 les données y restaurées conformément à l'invention, et sur les figures 11 et 12 les parties respectivement prédictible et non prédictible des données restaurées. La comparaison des figures 11 et 5 montre que les dix traces manquantes ne dégradent que peu le résultat final.
Portions de traces manguantes Il s'agit du cas où dans la section sismique t-x, on a des valeurs manquantes pour certains couples t, x connus (par exemple pour certaines traces, les valeurs y(t) correspondant à certains intervalles de temps sont manquantes).
Dans ce cas, on ne traite pas indépendamment chaque fréquence.
On a dans la section sismique t-x y=z-Pw y est ici l'ensemble des valeurs en t et x d'une section sismique, z contient l'ensemble des valeurs connues, w l'ensemble des données inconnues, et P est un opérateur linéaire qui positionne ces données inconnues dans l'ensemble des valeurs en t et x de la section sismique.
On utilise avantageusement l'algorithme du gradient conjugué
pour résoudre le système (D* D + 71,2 I) w = D* s . (21 ) Pour cela écrivons l'opérateur linéaire D sous la forme D=HFP (22) F représente l'opérateur linéaire qui effectue la transformation du plan t-x dans le plan f-x par transformée de Fourier de la variable t en F pour chaque x.
Pour appliquer l'algorithme du gradient conjugué au système linéaire (21), il n'est pas nécessaire d'écrire la matrice D, il suffit de savoir réaliser les opérations u' = D v' et r' = D* u' D'après l'expression (22), le procédé pour effectuer u' - D v' consiste à
1 ) prendre les valeurs de v' (v' est un vecteur quelconque qui a la dimension des valeurs manquantes) et les placer sur la section t-x avec des valeurs nulles pour ies valeurs connues (opérateur P') 2) effectuer la transformée de Fourier de t vers f pour tout x (appliquer l'opérateur F) 3 ) appliquer l'opérateur H pour tout f.

'~. WO 96/09562 217 718 0 P~~S/01231 Le procédé pour effectuer r' = D* u' consiste à
4) Appliquer l'opérateur H* pour tout f sur u' (u' étant une section dans le plan f-x) ) effectuer la transformée de Fourier inverse de f vers t pour tout x (sans 5 normaliser les valeurs en les divisant par le nombre d'échantillons (opérateur F*) 6 ) prendre sur la section sismique t-x obtenue les valeurs figurant aux emplacements manquants (opérateur P*).
On a d'autre part s=CFz s peut donc se calculer en appliquant aux données connues les étapes 2) et 3) ci-dessus, y-D* s se calculant à l'aide des étapes 4), 5), 6).
A partir de g et du procédé effectuant en cascade les étapes 1) à 6), l'algorithme du gradient conjugué résout le système (21).
Une fois w calculé, on obtient la section sismique t-x reconstituée par y=z-Pw On peut ensuite sparer y en une partie prdictible et une partie non prdictible dans le plan f-x pour raliser l'attnuation du bruit.

On peut galement, dans le cas o le filtre d'erreur de prdiction est inconnu, intgrer ce procd dans le processus N I1 en remplacement de l'tape 3) de calcul de y(n) dudit processus.

On peut aussi traduire l'tape 3' dans le cas de portions de traces manquantes.
On obtient ainsi l'organigramme suivant.
(processus N
12) 1 Initialisation des donnes y(t, n) ) 2 Passage de y(t, n) Y(f, n) par transforme de Fourier ) 3 Actualisation de a(f, i) : ' ) 3a) pour tout f i) Y(n) = Y(f~ n) ) calcul statistique du filtre d'erreur de prdiction a(f, i) pour y(n) et i =

0, ..., P

iii) initialisation de Y;i(f, n) = Y(f, n) pour tout n iv) retour en 3a)i) ou 4) si dernire frquence.

4 Actualisation de y(t, n) ) 4a) initialisation de y;t(t, n) = y(t, n) , 4b) passage de y;i(t, n) Y;t(f, n) par transforme de Fourier wo 96/09562 PCT/FR95/01231 4c) pour tout f i) Yit(n) = Yic(f~ n) ü) calcul de la partie non prédictible e(n) de yic(n) au moyen de l'opérateur M(f) calculé à partir des a(f, i) déjà calculés en 3c) iii) E(f, n) = e(n) iv) retour en i) ou 4d) si dernière fréquence 4d) passage de E(f, n) en e(t, n) par transformée de Fourier inverse 4e) actualisation de yic(t, n) pour les seuls échantillons (t, n) correspondant aux valeurs manquantes par Yic(t. n) _ (1-~2) Yic(t~ n) + ~2 Y(t~ n) - e(t~ n) 4f) retour en 4b) ou 4g) si convergence souhaitée de yic(t, n) 4g) actualisation de y(t, n) Y(>~ n) = Yic(t~ n) 5 ) Retour en 2) ou 6) si convergence de a(f, i) 6) Obtention des données reconstituées y(t,n) que l'on peut séparer en yo(t,n) et e(t,n) dans le domaine (f x).
Interpolation On peut remplacer une section sismique donnée t-x par une section dont la grille spatiale présente un pas plus petit, et calculer les traces manquantes en mettant en oeuvre le procédé selon l'invention.
On suppose dans la suite que seule une trace sur m est connue.
On applique avantageusement la méthode suivante Soit A(f) le filtre d'erreur de prédiction pour les données interpolées calculé de façon déterministe ou par une, autre méthode.
Soit Z(f) le spectre des données obtenues à partir des données connues supposées régulièrement espacées, en mettant à zéro les données manquantes situées sur une grille m fois plus fine que la grille des données connues, où m est la facteur d'interpolation. Soit Y(f) les données interpolées non débruitées que l'on cherche à estimer dans un premier temps.
Z(f) peut s'écrire Z(f) = Y(f) + Y(f + 1/m) + Y(f + 2/m) + .... + Y(f + (m-1)/m).
Si Y(f) contient p raies de fréquences f~, pour lesquelles A(f) s'annule, le signal Z(f) contiendra des fréquences f~ + i/m, i = 0, ..., m-1.
On construit le filtre d'erreur de prédiction AZ(f) = A(f + 1/m) A(f + 2/m) ... A(f + (m-1)/m) permettant d'enlever les fréquences fj + i/m pour i = 1, ..., m-1 et de wo Wo9562 21 l 718 0 p~~~p1231 récupérer Y(f) Dans le domaine temporel, aZ(i) se calcule par convolution des m filtres d'erreur de prédiction obtenus à partir de a(i) par modulation.
A(f + k/m), k = 1, ..., m - 1 est de longueur p dans le domaine temporel et a pour expression (a(0), a(1) e2in rem). a(2)e2in 2k/m), ..., a(p)e'-in p~i~n]
On calcule d'autre part l'autocorrélation de chaque filtre d'erreur de prédiction modulé, on ajoute un facteur de préblanchiment à chacun, et on convolue.
On obtient l'autocorrélation précolorisée rZ(i).
On peut donc estimer les données interpolées non débruitées Y(fj en appliquant à Z(f) le filtre d'erreur M(fj calculé à partir du filtre d'erreur de prédiction AZ(f). Ceci s'écrit dans le cas fréquentiel IAZ(f)l'-Y(f) = Z(f) R z(~
mais peut aussi se calculer dans le domaine temporel par une des structures de filtrage temporel décrites ci-dessus.
A partir des données interpolées non débruitées Y(f), on peut calculer le signal interpolé débruité Yo(f) à l'aide du filtre d'erreur de prédiction A(f), de son autocorrélation précolorisée et par une des méthodes décrites ci-dessous. Ceci s'écrit dans le cas fréquentiel jA(f)12 Yo(~ = Y(17 - Y(~
R(f) Une autre possibilité plus rapide est de faire Y(f) - Z(f) et de débruiter après, ce qui revient à écrire lA(t712 Yo(~ = Z(17 ~- Z(~
3 0 R(f) Cette méthode peut donner de bons résultats si le filtre jA( f)12 1 - présente une zone affaiblie large.
R(fj

Claims (28)

REVENDICATIONS

1/ Procédé de prospection sismique dans lequel:
- on provoque un ébranlement sismique dans le sous-sol, - on recueille à l'aide de capteurs des données sismiques échantillonnées y= [y(0),...,y(N)]T, où N est un nombre entier, contenant un signal y o que l'on cherche à isoler noyé dans un bruit additif, - on réalise dans le domaine temporel ou fréquentiel un filtrage des données sismiques, pour obtenir des données filtrées dans lesquelles le signal à
isoler est absent, - on soustrait aux données sismiques initiales les données filtrées pour obtenir des données traitées y o(0),...,y o(N) correspondant au signal y 0 dépourvu de bruit additif, - on traite les données ainsi acquises pour en déduire une information utile sur la géologie du sous-sol, le filtre utilisé pour le filtrage des données sismiques étant un filtre fréquentiel d'erreur de prédiction autodéconvolué M(f) tel que:
M(f) = ¦A(f)¦2 /R(f), ou un filtre équivalent dans le domaine temporel, A(f) étant le spectre d'un filtre a d'erreur de prédiction à p+1 coefficient a(0), ...,a(p), où p est un nombre entier inférieur à N, a étant préalablement choisi pour annuler au mieux par convolution le signal y 0 que l'on cherche à
isoler, R(f) étant une autocorrélation précolorisée de ce filtre de prédiction a vérifiant :
R(f) = ¦A(f)¦2 + .epsilon.2 ¦B(f)¦2, où .epsilon. et B(f) sont respectivement un facteur et un filtre de précolorisation préalablement choisis en fonction de la sélectivité souhaitée pour le filtrage.

2/ Procédé selon la revendication 1, caractérisé en ce que le filtre d'erreur de prédiction autodéconvolué est non stationnaire pour des données sismiques proches des bords.

3/ Procédé selon la revendication 1, caractérisé en ce que les p+1 coefficients a(0), a(1),..., a(p) du filtre d'erreur de prédiction a sont calculés de façon déterministe à partir de paramètres physiques connus du signal y o que l'on cherche à isoler.

4/ Procédé selon la revendication 3, le signal y o présentant un spectre Y o(f) constitué de p raies de fréquences f q (q = 1, ..., p), caractérisé en ce que les coefficients a(0), a(1), ..., a(p) du filtre d'erreur de prédiction .alpha.
sont calculés par a(i) = a1(i) * a2(i) * ... * a p(i) a(0) = 1 avec a q(0) = 1 a q(1) = - e-2j.pi.f q

5/ Procédé selon la revendication 3, le signal y o présentant un spectre Y o(f) dont le support comporte au moins un intervalle de fréquence [f q - .DELTA.f q, f q + .DELTA.f q]

.DELTA. f q étant la demi-largeur de l'intervalle centré sur la fréquence f q, caractérisé en ce que les coefficients a(0), a(1), ... , a(p) du filtre d'erreur de prédiction .alpha., de spectre A(f), sont déterminés par une méthode des moindres carrés, en minimisant l'intégrale .intg. W(f) ¦A(f)¦2 df, avec a(0) = 1 et W(f) une fonction qui vaut 1 lorsque la fréquence f est dans un intervalle de fréquences où le spectre Y o(f) est non nul, et 0 ailleurs.

6/ Procédé selon la revendication 3, le signal y o présentant un spectre Y o(f) dont le support comporte Q intervalles de fréquences [f q - .DELTA.f , f q + .DELTA.f q) q = 1 ,..., Q
.DELTA. f9 étant la demi-largeur d'un intervalle centré sur la fréquence f q, caractérisé en ce que les coefficients du filtre d'erreur de prédiction .alpha. sont déterminés par a(i) = a1(i) * a2(i) * ... * a Q(i) avec a q(i) = s q(i) / s q(0) pour q = 1, ..., Q
les coefficients sq(i), i = 0,..., p q étant déterminés par la relation S q(f) = e-J.pi.fp T pq[sin(.pi.(f-f q)) / sin(.pi. .DELTA.f q)] = .SIGMA. s q(m) e-j2.pi.mf où T pq est le polynôme de Chebychev de première espèce et d'ordre p q et en ce que l'autocorrélation précolorisée est calculée par :

et

7/ Procédé selon la revendication 1, caractérisé en ce que le filtrage consiste dans le domaine fréquentiel à multiplier les données sismiques y(0),..., y(N) par le filtre M(f).

8/ Prodécé selon la revendication 1, caractérisé en ce que, dans le cas d'un filtrage dans le domaine temporel, on calcule l'autocorrélation précolorisée r(i) correspondant dans le domaine temporel à R(f) et telle que :

où b est le filtre de précolorisation équivalent à B(f) dans le domaine temporel.

9/ Procédé selon la revendication 8, caractérisé en ce que le filtrage par le filtre d'erreur de prédiction auto-déconvolué comporte les étapes suivantes :
1) Calcul d'un filtre 2) Calcul par filtrage non-récursif pour n variant de 0 à N de e(n), tel que :
e(n) étant le bruit additif.

10/ Procédé selon la revendication 8, caractérisé en ce que le filtrage par le filtre d'erreur de prédiction auto-déconvolué comporte les étapes suivantes :
1) Calcul du filtre c(i) tel que 2) Calcul par filtrage récursif ascendant pour n variant de 0 à N de u(n) tel que :
3) Calcul par filtrage récursif descendant pour n variant de 0 à N de e(n) tel que :

e(n) étant le bruit additif.

11/ Procédé selon la revendication 8, caractérisé en ce que l'application aux données du filtre d'erreur de prédiction auto-déconvolué
comporte les étapes suivantes :

1) Calcul d'un filtre non stationnaire G partir de l'autocorrélation précolorisée r(i), G étant tel que G* G = R-1, R étant la matrice de Toeplitz formée à partir de r(i) 2 Calcul par filtrage non-récursif stationnaire de s - Ay correspondant à
s(n) = a(i) * y(n) 3) Calcul par filtrage non-récursif non stationnaire de u = Gs 4) Calcul par filtrage non-récursif non stationnaire de v = G* u 5) Calcul par filtrage non-récursif stationnaire de e = A* v correspondant à

e(n)tant le bruit additif.

12/ Procédé selon la revendication 8, caractérisé en ce que l'application aux données sismiques du filtre d'erreur de prédiction auto-déconvolué comporte les étapes suivantes :

1) Calcul du filtre non stationnaire C à partir de l'autocorrélation précolorisée r(i), C étant tel que CC* = R, R étant la matrice de Toeplitz formée à partir de r(i) 2) Calculé par filtrage non-récursif stationnaire de s - Ay correspondant à

s(n) = a(i) * y(n) 3) Calcul par filtrage récursif non stationnaire de u = C-1s 4) Calcul par filtrage récursif non stationnaire de v = C-1*u 5) Calcul par filtrage non-récursif stationnaire de e = A*v correspondant à
e(n)tant le bruit additif.

13/ Procédé selon la revendication 1, caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) du filtre d'erreur de prédiction a sont calculés de façon statistique par itérations successives, par minimisation d'une norme dépendant de a et des donnes sismiques y.

14/ Procédé selon la revendication 13, caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) sont calculés par le processus suivant :
1) Calcul de Y(f) par transforme de Fourier des donnes y(n) 2) Choix d'un filtre d'erreur de prédiction initial a, de coefficients a(0) =
1, a(1),., a(p) 3) Calcul de A(f) par transformée de Fourier de a(i) 4) Calcul de UU(f) = ¦Y(f)¦2 / [¦A(f)¦2 + .EPSILON.2 ¦B(f)¦2-]
5) Calcul de uu(i) par transformée de Fourier inverse de UU(f) 6) Calcul de a(i) en résolvant le système d'équations de Toeplitz transformé
avec uu(i) 7) Reprise à l'étape 2) jusqu'à la convergence souhaitée des valeurs a(i).

15/ Procédé selon les revendications 8 et 13, caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) sont calculés par le processus suivant 1) Choix d'un filtre d'erreur de prédiction initial .alpha., de coefficients a(0), a(1), ..., a(p) 2) Calcul de l'autocorrélation précolorisée 3) Calcul de g(i) à partir de r(i), g(i) étant tel que avec .delta.(i) = 1 si i = 0 et .delta.(i) = 0 sinon 4) Calcul de u(n) = g(i) * y(n) 5) Calcul de la matrice de covariance COV u de u(n) avec i = 0....,p et j = 0....,p 6) Calcul de h(0), .. , h(p) vérifiant COV u h = [1, 0, ..., 0]*
7) Calcul de a(i) par normalisation de h(i):
a(i) = h(i) / h(0) 8) Reprise à l'étape 2) jusqu'à convergence souhaitée des valeurs a(i).

16/ Procédé selon les revendications 8 et 13, caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) sont calculés par le processus suivant:
1) Choix d'un filtre d'erreur de prédiction initial .alpha., de coefficients a(0), a(1), ..., a(p) 2) Calcul de l'autocorrélation précolorisée où b(i) est un filtre de précolorisation 3) Calcul du filtre d'erreur de prédiction amorti c(i) à partir de r(i) par un processus itératif de filtrage temporel récursif 4) Calcul par filtrage récursif de u(0), ..., u(N) par 5) Calcul de la matrice de covariance COV u de u(n) pour i = 0....,p et j = 0.....,p 6) Calcul du vecteur h de composantes h(0), ..., h(p) vérifiant COV u h = [1, 0, ..., 0]*
7) Calcul de a(i) par normalisation de h(i):
a(i) = h(i) / h(0) 8) Reprise à l'étape 2) jusqu'à la convergence souhaitée des valeurs a(i).

17/ Procédé selon les revendications 8 et 13, caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) sont calculés par le processus suivant 1) Choix de p+1 coefficients a(i) initiaux 2) Calcul du filtre non stationnaire G à partir de l'autocorrélation précolorisée r(i), G étant tel que G*G = R-1, R étant la matrice de Toeplitz formée à partir de r(i) 3) Calcul par filtrage non-récursif non-stationnaire des p+1 vecteurs u i : u i = G Yi 4) Calcul de la matrice de covariance des u i : COV u(i, j) = u i * u j 5) Résolution de 6) Calcul de a(i) = h(i) / h(0) i = 0, ..., p 7) Reprendre à l'étape 2) jusqu'à la convergence souhaitée des valeurs a(i) 8 ) Fin.

18/ Procédé selon les revendications 8 et 13 caractérisé en ce que les p+1 coefficients a(0), a(1), ..., a(p) sont claculés par le processus suivant:

1) Choix de p+1 coefficients a(i) initiaux 2)Calcul du filtre non stationnaire C à partir de l'autocorélation précolorisé r(i), C étant tel que CC* = R, R étant la matrice de formée à partir de r(i) 3)Calcul par filtrage récursif non-stationnaire des p+1 vecteurs u i tels que:

u i = C-1yi avec yi = [y(p-i),..., y(N-i)] T pour i = 0 à p 4) Calcul de la matrice de covariance des u i : COV u(i, j) = u i * u j 5) Résolution de 6) Calcul de a(i) = h(i) / h(0) i = 0, . . ., p 7) Reprendre à l'étape 2) jusqu'à la convergence souhaitée des valeurs a(i) 8) Fin.

19/ Procédé selon l'une quelconque des revendications 2, 3, 6, 8 et 12, caractérisé en ce qu'il comprend les étapes suivantes:
1) Détermination par l'utilisateur d'une zone du plan f-k contenant le signal y o à isoler, 2) Pour toute valeur de la variable espace x, effectuer la transformée de Fourier des données y(t, x) de la variable temporelle t pour obtenir des données transformées Y(f, x) dans le plan f-x, 3) Pour chaque fréquence f du plan f-x, 3 a) prendre y(n) = Y(f,n), 3b) déterminer les intervalles pour la variable k, contenant le signal prédictible y o que l'on cherche à isoler, en effectuant une coupe à f constant de ladite zone, 3c) calculer le filtre d'erreur de prédiction spatial élémentaire (variable x) correspondant à chacun de ces intervalles en k pour un .epsilon. donné, 3d) calculer le filtre d'erreur de prédiction global a(i) et l'auto-corrélation précolorisée correspondante r(i), 3e) calculer le filtre d'erreur de prédiction non-stationnaire ep(j,i) à
partir de l'autocorrélation précolorisée r(i), 3f) calculer la partie non prédictible e(n) des données par filtrage récursif non stationnaire à l'aide des coefficients ep(j.i), 3g) soustraire la partie non prédictible e(n) aux données y(n) pour obtenir la partie prédictible y o(n), 3h) reprendre à l'étape 3a) pour la fréquence f suivante sauf si la dernière fréquence du plan f-x est atteinte, 4) Pour tout x, effectuer la transformée de Fourier inverse de la variable f de e(n) ou de y o(n) pour revenir dans le plan t-x.

20/ Procédé selon l'une quelconque des revendications 2, 8, 12, 13 et 16, caractérisé en ce qu'il comprend les étapes suivantes:
1) Pour tout x, transformée de Fourier des données sismiques y(t, x) pour la variable t pour se placer dans le plan f-x 2) Pour chaque fréquence f du plan f-x:
2a) prendre y(n) = Y(f, n) 2b) choix d'une longueur p pour le filtre d'erreur de prédiction, d'un facteur .epsilon. et d'un filtre de précoloriement b(0), ..., b(q), ainsi que d'un nombre d'itérations pour la convergence de a(i) 2c) calcul par filtrage récursif ascendant de:
c(i) étant le filtre d'erreur de prédiction amorti de l'itération précédente (celui de la dernière itération de la fréquence précédente si c'est la première itération pour la fréquence en cours d'itération) 2d) calcul de la matrice de covariance ou de corrélation de u(n) 2e) résolution du système pour obtenir les coefficients a(i) du filtre d'erreur de prédiction, 2f) calcul de l'autocorrélation précolorisée r(i) et du filtre d'erreur de prédiction amorti c(i) 2g) retour en 2c) jusqu'à la convergence souhaitée des valeurs de a(i) 2h) obtention de la partie non prédictible e(n) des données par filtrage récursif non stationnaire 2i) soustraction de la partie non prédictible e(n) aux données y(n) pour obtenir la partie prédictible y o(n) 2j) retour en 2a) sauf si dernière fréquence 3) Pour tout x, transformée de Fourier inverse de y o(n) de la variable t pour obtenir y o(t, x).

21/ Procédé selon la revendication 1, appliqué à la restauration de données présentant des valeurs manquantes d'indice I(j), j étant compris entre 1 et m, caractérisé en ce qu'il comprend les étapes suivantes:
1) Initialisation des données à restaurer y(n) 2) Calcul du filtre d'erreur de prédiction a(i) pour le y(n) courant 3) Actualisation de y(n) pour le a(i) courant:
3a) calcul de m vecteurs d j = H p j j = 1, ..., m p j étant un vecteur contenant la réponse impulsionnelle du j ième bruit impulsif correspondant à la j ième valeur manquante d'indice I(j), H étant l'opérateur linéaire tel que M - H* H. M étant l'opérateur permettant de calculer le bruit e à partir des données y 3b) calcul du vecteur s = H y, 3c) calcul de la matrice .chi.(i, j) = d i * d j i,j = 1, ..., m 3d) addition d'un facteur de préblanchiment .lambda.2:

.chi. (i, i) .rarw. .chi. (i,i) + .lambda.2 i = 1, ..., m 3e) calcul du vecteur de composantes k(i) = d i * s 3f) résolution du système linéaire d'ordre m:
m .SIGMA.r(i, j) w(j) = k(i) i = 1, ..., m j=1 3g) actualisation de y(n) par:
m y~.rarw. y -.SIGMA. w(j) p j.
j=1 4) Retour en 2) ou 5) jusqu'à la convergence souhaitée pour a(i) 5) Calcul de u:
m u=s-.SIGMA. w(j)d j.
j=1 6) Calcul de e = H* u 7) Calcul de y o = y - e.

22/ Procédé selon la revendication 21, caractérisé en ce que l'on prend pour p j les composantes d'un bruit cohérent d'amplitude inconnue.

23/ Procédé selon la revendication 22, caractérisé en ce que l'on prend:
p i(n) = e jn.DELTA.xk i(f) avec n = 0, ..., N, k i(f) étant la relation de dispersion vérifiée par le bruit.

24/ Procédé selon la revendication 21, caractérisé en ce que l'on remplace l'étape 3) par l'étape 3') suivante:
3') Actualisation de y(n) pour le a(i) courant:
3a) initialisation de y it(n) = y(n) 3b) calcul de la partie non prédictible e(n) de y it(n) 3c) actualisation de y it(n) pour les seuls échantillons n correspondant aux valeurs manquantes (n = I(j) j = 1, ..., m) par:
y it(n) .rarw. (1 - .lambda.2) y it(n) + .lambda.2 y(n) - e(n) 3d) retour en 3b) ou 3e) si convergence souhaitée de y it(n) 3e) actualisation de y(n) pour les seuls échantillons n correspondant aux valeurs manquantes par:
Y(n) = Y it(n).

25/ Procédé selon la revendication 21, appliqué à la restauration de portions de traces manquantes, caractérisé en ce qu'il comprend les étapes suivantes :
1) Initialisation des données y(t, n) 2 ) Passage de y(t, n) à Y(f, n) par transformée de Fourier 3 ) Actualisation de a(f, i) 3a) pour tout f i) y(n) = Y(f, n) ü) calcul statistique du filtre d'erreur de prédiction a(f, i) pour y(n) et i =
0, ..., P
iii) initialisation de Y it(f, n) = Y(f, n) pour tout n iv) retour en 3a)i) ou 4) si dernière fréquence.

4) Actualisation de y(t, n) 4a) initialisation de y it(t, n) = y(t, n) 4b) passage de y it(t, n) Y it(f, n) par transformée de Fourier 4c) pour tout f :
i) Y it(n) = Y it(f, n) ii) calcul de la partie non prédictible e(n) de y it(n) au moyen de l'opérateur M(f) calcul partir des a(f, i) déjà calculés en 3c) iii) E(f, n) = e(n) iv) retour en i) ou 4d) si dernière fréquence 4d) passage de E(f, n) en e(t, n) par transformée de Fourier inverse 4e) actualisation de y it(t, n) pour les seuls échantillons (t, n) correspondant aux valeurs manquantes par :
Y it(t, n) = (1-.lambda.2) y it(t, n) + .lambda.2 y(t, n) - e(t, n) 4f) retour en 4b) ou 4g) si convergence souhaitée de y it(t, n) 4g) actualisation de y(t, n) y(t, n) = y it(t, n) 5) Retour en 2) ou 6) si convergence de a(f, i) 6) Obtention des données reconstituées y(t,n) que l'on peut séparer en y o(t,n) et e(t,n) dans le domaine (f-x).

26/ Procédé selon la revendication 1, appliqué à l'interpolation de données sismiques, caractérisé en ce qu'il comporte les étapes suivantes, implémentées dans le domaine fréquentiel :
1 ) Calculer un filtre d'erreur de prédiction A(f) pour les données interpolées, 2) Calculer Z(f) par transformée de Fourier des données connues, supposées régulièrement espacées, en mettant à 0 les données manquantes situées sur une grille m fois plus fine que la grille des données connues, où m est le facteur d'interpolation 3) Construire le filtre d'erreur de prédiction A z(f) = A(f + 1/m) A(f + 2/m) ... A(f + (m-1)/m) 4) Calculer les données interpolées non débruitées par R z(f) étant une autocorrélation précolorisée de A z(f).

27/ Procédé selon la revendication 1, appliqué à l'interpolation de données sismiques, caractérisé en ce qu'il comporte les étapes suivantes, implémentées dans le domaine fréquentiel :
1) Calculer un filtre d'erreur de prédiction A(f) pour les données interpolées, 2) Calculer Z(f) par transformée de Fourier des données connues, supposées régulièrement espacées, en mettant à 0 les données manquantes situées sur une grille m fois plus fine que la grille des données connues, où m est le facteur d'interpolation 3) Calculer les données interpolées débruitées par

28/ Procédé selon la revendication 26 ou 27, caractérisé en ce que les étapes du procédé sont mises en oeuvre dans un domaine temporel.