PROCEDE DE PROSPECTION SISMIQUE AVEC APPLICATION D'UN FILTRE D'ERREUR DE PREDICTION AUTO-DECONVOLUE

SEISMIC PROSPECTION METHOD WITH APPLICATION OF A SELF-DECONVOLUTED PREDICTION ERROR FILTER

English 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.

French Abstract

L'invention concerne un procédé de prospection sismique dans lequel on
provoque un ébranlement sismique dans le sous-sol et on recueille à l'aide de
capteurs des données sismiques, ces données sismiques étant traitées en vue
d'en tirer une information utile sur la géologie du sous-sol, les données
sismiques contenant un signal yo que l'on cherche à isoler, ce signal yo étant
présent dans les données sismiques noyé dans un signal additif. Le procédé
comprend des étapes consistant à: 1) calculer pour un entier p donné un filtre
d'erreur de prédiction a ayant p coefficients et annulant au mieux le signal
yo; 2) appliquer aux données sismiques le filtre d'erreur de prédiction auto-
déconvolué pour obtenir des données filtrées dans lesquelles le signal yo est
absent; 3) soustraire aux données initiales les données filtrées pour obtenir
des données traitées contenant le signal yo sans le signal additif.

Claims

40
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

41
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 :
<IMG>
et

42
<IMG>
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 :
<IMG>
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 <IMG>
2) Calcul par filtrage non-récursif pour n variant de 0 à N de e(n), tel que :
<IMG>
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 <IMG>
2) Calcul par filtrage récursif ascendant pour n variant de 0 à N de u(n) tel
que :
<IMG>
3) Calcul par filtrage récursif descendant pour n variant de 0 à N de e(n) tel
que :
<IMG>
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 :

43
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 à
<IMG>
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 à
<IMG>
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

44
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 <IMG>
3) Calcul de g(i) à partir de r(i), g(i) étant tel que <IMG>
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)
<IMG>
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 <IMG>
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
<IMG>
5) Calcul de la matrice de covariance COV u de u(n)
<IMG>

45
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
<IMG>
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
<IMG>

46
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)

47
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:
<IMG>
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:

48
.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).

49
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,

50
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
<IMG>
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
<IMG>
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.

Description

2 ~ 7 718 0 p~~~p1231
WO 96/09562
1
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
X177180
2
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
3
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
WO 96/09562 PCT/FR95/01231 -
2177180
4
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
PCT/FR95101231
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
6
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
2177180
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
M-1
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
PCT/FR95I01231
9
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
2177180
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
11
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
12
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

PCT/FR95I01231
~''~ wo ~'°~2 21 7 71 8 0
13
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
117180
14
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
PCT/FR95/01231
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
16
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
PCT/FR95/01231
17
. 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~
18
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
PCT/FR95J01231
19
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
2177180
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.

WO 96/09562
PCT/FR95/01231
2171180v
21
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

WO 96109562 PCT/FR95/01231
21~7~80
22
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
23
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.
0
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 ô ~
PCT/FR95/01231
24
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.

2177180
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

W0 96109562 PCT/FR95/01231
2177180
26
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
27
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
28
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
29
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

WO 96/09562 PCT/FR95/01231
2177180
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

PCT/FR95/01231
wo 9b/o9s62 217 718 0
31
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

WO 96/09562 PCT/FR95/01231
2177180
32
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é.

2177180
wo mo~2 rc~r~siom
33
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

WO 96109562 PCT/FR95/01231
2177180
34
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

217718Q
WO 96109562 PCT/FR95/01231
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
2177180
36
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
37
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
2171180
38
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
39
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

Representative Drawing