CA2790049A1 - Procede de pointage d'horizons sismiques discontinus dans des images sismiques - Google Patents

Abstract

L'invention vise en particulier un procédé de recherche d'un horizon sismique dans une image sismique du sous-sol. Le procédé comprend notamment la designation de deux points appartenant à l'horizon recherché, la recherche itérative de la meilleur solution pour l'équation de l'horizon sismique notamment grâce à une décomposition de cette solution en deux composantes, une composante pseudo-continue et une composante de saut.

Description

PROCÉDÉ DE POINTAGE D'HORIZONS SISMIQUES DISCONTINUS DANS
DES IMAGES SISMIQUES

La présente invention concerne le domaine de l'analyse d'images sismiques du sous-sol.

II est connu, notamment dans l'exploration pétrolière, de déterminer la position des réservoirs pétroliers à partir des résultats de mesures géophysiques effectuées depuis la surface ou dans des puits de forage. Dans la technique de la sismique en réflexion, ces mesures impliquent l'émission dans le sous-sol d'une onde et la mesure d'un signal comportant diverses réflexions de l'onde sur les 1o structures géologiques recherchées. Ces structures sont typiquement des surfaces séparant des matériaux distincts, des failles, etc. D'autres mesures sont effectuées depuis des puits. On envoie alors dans le sous-sol, des ondes acoustiques ou un rayonnement électromagnétique.

Les mesures sont traitées pour reconstituer un modèle du sous-sol, en général sous forme d'images sismiques, ou images échographiques. Ces images peuvent être bidimensionnelles (sections sismiques) ou tridimensionnelles (blocs sismiques). Une image sismique se compose de pixels dont l'intensité est représentative d'une amplitude sismique dépendant des variations locales d'impédance.

Les géophysiciens sont habitués à analyser de telles images sismiques porteuses d'information d'amplitude. Par l'observation visuelle, ils sont capables de séparer des zones du sous-sol ayant des caractéristiques différentes en vue de déterminer la structure du sous-sol.

Il existe des méthodes automatiques pour extraire de l'information structurelle d'images sismiques. En particulier, il est connu d'estimer des horizons sismiques en analysant par ordinateur les gradients d'amplitude dans l'image sismique.
Les horizons ainsi estimés sont appelés "horizons de synthèse" par opposition aux horizons déterminés par pointage manuel des images.

Une méthode possible pour estimer des horizons dans une image sismique

2 bidimensionnelle consiste, à partir d'un pixel de l'image, à rechercher la direction selon laquelle le gradient local d'amplitude est minimal. En se propageant le long de cette direction, on construit de proche en proche une ligne représentant un horizon de synthèse. Si l'image sismique est tridimensionnelle, on peut estimer des horizons sous forme de surfaces transversales à la direction verticale, par exemple au moyen de la méthode de propagation décrite dans le brevet français n 2 869 693.

Néanmoins, si le pointé automatique d'horizon sismique est basé sur le suivi de continuités dans une image ou un bloc sismique, la présence d'une 1o discontinuité, telle qu'une faille, peut empêcher le suivi de l'horizon ou peut introduire un biais dans le pointé de l'horizon.

Il y a ainsi un besoin pour faciliter le pointage d'horizon sismique en cas de présence de discontinuités.

A cet effet, la présente invention propose un procédé de recherche d'un horizon sismique dans une image sismique du sous-sol. Ce procédé comprend:
- designer deux points de coordonnées respectives x1, y1 et xN, YN dans l'image sismique comme étant deux points appartenant à l'horizon recherché, un nombre entier N d'abscisses discrètes x1, x2, ..., XN étant définies entre les abscisses x1, xN des deux points désignés;

- considérer une position de faille à une abscisse discrète xna où na est un entier compris entre 1 et N-1 et une amplitude de faille Cd dans un domaine discret;
- en se donnant deux fonctions f et g0 dérivables sur l'intervalle [x1, xN], initialiser une amplitude de faille dans un domaine continu à la valeur CO = Cd - [f(xna+1) -f(xna)], une composante pseudo-continue To(xn) et une composante de saut î0(xn) définies sur les N abscisses discrètes selon i0(x)=f(x) et i0 (xn) = g0 (xn) + C0.H(n-n(x) , où H(.) désigne la fonction de Heaviside;
- effectuer plusieurs itérations d'une séquence d'étapes de calcul pour un index d'itération k initialisé à k = 0 puis incrémenté par unités, la séquence

3 comprenant, pour un index k:

-calculer un résidu rk(xn)=oik(xn)+Vgk(xn)-p(xn,ik(xn)+îk(xn)), où V désigne l'opérateur gradient et p(xn, y) désigne la tangente d'un pendage estimé pour une position d'abscisse discrète xn et d'ordonnée y dans l'image sismique;
- résoudre une équation de Poisson A(âik) = -div(rk) pour déterminer un terme de mise à jour OTk (xn) , avec des conditions aux limites de Dirichlet: âto (x1) Y1 - f (x1) - go (x1) à la première itération SiO(xN) = YN -f(xN)-g0(xN)-CO

btk (x1) = gk-1(x1) - gk (x1) et a chaque itération Sik(xN) = gk-1(xN)-gk(xN)+Ck-1 -Ck d'index k > 1, où gk est une fonction dérivable sur l'intervalle [x1, xN];

- mettre à jour la composante pseudo-continue selon ik+1(xn) = ik(xn)+Sik(xn);

- mettre à jour l'amplitude de faille dans le domaine continu selon Ck+1 =Cd - hk+1(xna +1) - Tk+1(xna ), - mettre à jour la composante de saut selon tk+1(xn) = gk+1(xn) + Ck+1 =H(n-na) ;

- si une valeur finale de l'index d'itération k est atteinte, calculer une fonction ina Cd(xn) représentant l'ordonnée dans l'image sismique d'un horizon estimé en fonction de l'abscisse discrète xn, à partir d'une somme de la composante pseudo-continue et de la composante de saut; et - si la valeur finale de l'index d'itération k n'est pas atteinte, exécuter l'itération suivante de la séquence.

Pour toute valeur de l'index d'itération k, les fonctions gk sont prises comme étant des fonctions dérivables sur l'intervalle [x1, xN]. Le choix de ces fonctions peut être libre et peut dépendre de l'implémentation dans un mode de réalisation particulier. En particulier, il n'existe pas nécessairement de formule générale

4 d'évolution de la suite des fonctions gk.

Le pointé automatisé d'un horizon sismique discontinu repose sur une segmentation en deux parties de l'horizon par une faille supposée quasi-verticale.
Il faut connaître deux points délimitant l'horizon à reconstruire dans l'image sismique, et connaître ou faire une hypothèse sur le lieu et l'amplitude de la discontinuité. Le pointé s'obtient par la résolution d'une équation aux dérivées partielles non-linéaire liant la solution recherchée à l'orientation locale de l'horizon. Le procédé met en jeu un algorithme précis, robuste au bruit et présentant des propriétés de convergence rapide. Il permet, si nécessaire, de 1o mettre en place un schéma ad hoc efficace sur des données sismiques réelles sans connaissance a priori du lieu et de l'amplitude de la discontinuité.

Dans un mode de réalisation du procédé, la fonction f est choisie comme une fonction constante. Elle peut notamment être prise égale à Y1 + R-(YN - Y1), R
étant un nombre réel choisi entre 0 et 1. Dans ce cas, chaque fonction gk peut aussi être prise constante, en particulier de valeur -t3.Ck.

Pour résoudre l'équation de Poisson, on peut utiliser une transformée de Fourier discrète DFT et sa transformée de Fourier discrète inverse DFT_1 selon âik = DFT-1 _ DFT[div(rk)] En particulier, on peut utiliser une transformée en DFT[A]

sinus discrète de type I.

L'invention a également pour objet un programme d'ordinateur pour un système de traitement d'images sismiques du sous-sol, le programme comprenant des instructions pour mettre en oeuvre les étapes d'un procédé tel que défini ci-dessus lors d'une exécution du programme par un calculateur du système de traitement d'images sismiques.

D'autres caractéristiques et avantages de l'invention apparaîtront encore à la lecture de la description qui va suivre. Celle-ci est purement illustrative et doit être lue en regard des dessins annexés sur lesquels :
- la figure 1 est une portion d'image sismique synthétisée, en noir et blanc;
- la figure 2 est un graphique illustrant un horizon sismique ayant une discontinuité et introduisant des notations utilisées dans un mode de réalisation selon l'invention;
- la figure 3 est un organigramme d'un exemple d'algorithme de résolution itératif dans le domaine continu pour la détermination d'un horizon

5 sismique présentant une discontinuité;
- la figure 4 est un graphique illustrant un horizon sismique dans le domaine discret, présentant une discontinuité;
- la figure 5 est un organigramme d'un exemple d'algorithme de résolution itératif dans le domaine discret pour la détermination d'un horizon sismique présentant une discontinuité;
- la figure 6 est un organigramme d'algorithme généralisé de résolution itératif dans le domaine discret pour la détermination d'un horizon sismique présentant une discontinuité;
- la figure 7 est un organigramme d'un exemple d'algorithme de résolution équivalent à celui de la figure 6, dans lequel les index d'itérations k ne figurent pas explicitement; et - la figure 8 montre un organigramme d'un exemple de procédé de recherche d'horizon sismique mettant en oeuvre l'un des algorithmes précédemment illustrés.

En référence à la figure 1, une image sismique se présente couramment sous la forme d'une matrice de pixels dont l'intensité correspond à une amplitude sismique. Les variations d'amplitude sismique traduisent des variations d'impédance vis-à-vis de la propagation des ondes acoustiques dans le sous-sol.
La représentation est souvent en deux couleurs (non rendues sur la figure), une pour les amplitudes négatives, l'autre pour les amplitudes positives. L'image, lorsqu'elle est bidimensionnelle, peut correspondre à une intégration à deux dimensions de traces sismiques verticales. Elle peut aussi correspondre à une coupe verticale dans un bloc sismique obtenu par intégration à trois dimensions de traces sismiques verticales.

Dans la portion d'image sismique représentée à titre d'exemple sur la figure 1, les abscisses x correspondent à une direction horizontale, tandis que les

6 ordonnées y peuvent correspondre soit à la direction verticale si l'image est migrée en profondeur, soit au temps si l'image est migrée en temps. Il est possible d'y observer des horizons sismiques qui, dans cet exemple, s'étendent de part et d'autre d'une faille verticale située à une abscisse x = . Un pointé
manuel d'horizons peut être effectué par un opérateur par une observation des textures visibles sur l'image, encore qu'un tel pointé puisse parfois être ambigu en fonction de la complexité de l'image, spécialement de ses discontinuités, et du bruit qu'elle contient.

Il est souhaitable d'effectuer un pointé d'horizons automatisé à l'aide d'une 1o technique de pointé informatisée afin d'améliorer l'efficacité et la productivité
dans le traitement des images sismiques. Toutefois les techniques de propagation couramment mises en oeuvre pour réaliser un tel pointé automatisé
peuvent être mises en défaut du fait des discontinuités engendrées notamment par des failles.

L'opérateur (interprétateur d'images sismiques) est typiquement sollicité pour indiquer deux points 11, 12 de coordonnées respectives (x1, y,) et (XN, YN) dans l'image, supposés appartenir au même horizon de part et d'autre d'une faille 10.
La difficulté est alors de déterminer l'horizon le plus plausible (représenté
par une ligne interrompue 15 sur la figure 1) passant par ces deux points 11, 12 compte tenu des données de l'image.

Un horizon sismique ayant une discontinuité entre les deux abscisses x1, xN
des points 11, 12 pointés par l'opérateur de part et d'autre de la faille 10 peut être modélisé par une courbe illustrée par la Figure 2, représentant une fonction T
définie pour des abscisses x comprises entre les abscisses x1 et xN. La valeur i(x) désigne la profondeur de l'horizon sismique à l'abscisse x dans le cas où
l'image est migrée en profondeur, ou sa position temporelle à cette abscisse x dans le cas où l'image est migrée en temps, avec T(x1) = Y1 et T(xN) = yN. La fonction i présente un saut d'amplitude Cc = i+ - i_ à l'abscisse a E ]x1, xN[, où
T+ = lim T(x) est la valeur limite de la fonction i lorsque x tend vers a par x-+a+

7 valeurs supérieures et i_ = lim i(x) est la valeur limite de la fonction r lorsque x X->-tend vers a par valeurs inférieures. L'amplitude Cc est l'amplitude de la faille 10 dans un domaine continu.

Le pointé d'horizon automatisé consiste donc à déterminer la meilleure fonction i compte tenu des données de l'image sismique. Ces données peuvent être prétraitées pour estimer un pendage 6(x, y) à chaque position x, y de l'image sismique. De façon classique, ce prétraitement peut par exemple consister à
rechercher, en chaque position x, y, la direction DX y = (u, v) suivant laquelle les valeurs de l'image sismique présentent un gradient d'amplitude maximum et à
1o représenter le pendage 9(x, y) par sa tangente p(x, y), c'est-à-dire l'opposé de la cotangente de l'angle formé par la direction DX y de gradient maximum par rapport à l'axe des abscisses de l'image (p(x, y) = -u/v), comme illustré
schématiquement pour le point 16 sur la figure 1. Naturellement, les valeurs de pendage obtenues de cette manière au niveau de la discontinuité 10 ne sont pas fiables.

Pour la recherche de l'horizon sismique, il est proposé de résoudre une équation aux dérivées partielles de manière itérative. Sauf en x = a, la fonction T
et la tangente p du pendage local sont liés par l'équation suivante VT(x) = p(x, i(x)) (1) où V est l'opérateur gradient. Aux fins de l'explication, les fonctions T et p sont supposées être respectivement de classe C2 et Cl sauf en x = a. Une fonction est de classe Cl si elle est continue. Une fonction est de classe C2 si elle est dérivable et sa dérivée de classe Cl. L'équation (1) étant non-linéaire, le recours à un algorithme itératif est nécessaire pour déterminer la solution recherchée. Du fait de la discontinuité, il est proposé de décomposer la fonction i qui est de classe C2 sur [x1, xN] - {(x} en la somme d'une composante continue inconnue i de classe C2 sur [xl, xN] vérifiant l'équation (1) et d'une composante de saut discontinue î constante sur [xl, a[ et constante également sur ]a, xN]:

8 î(x) = -(3.Cc.H(a-x) + (1-I3).Cc.H(x-a) (2) H étant la fonction de Heaviside (H(x) = 0 si x< 0 et H(x) = 1 sinon), et p étant un paramètre réel, par exemple choisi entre 0 et 1.

Ainsi, un exemple d'algorithme de résolution itératif est représenté sur la figure 3. L'algorithme a pour objet de calculer une suite de composantes continues îk convergeant efficacement vers la solution de l'équation (1) pour un index d'itération k = 0, 1, 2, ...

A l'initialisation 20, où l'index k est pris égal à zéro, l'algorithme fixe la fonction îp de la suite, par exemple à une valeur constante:

ti0 (x) = Y1 + P (YN-Y1) (3) prise entre les ordonnées y, et YN des points 11, 12 si le paramètre f3 est choisi entre 0 et 1. La fonction î(x) est alors fixée à la valeur î(x)=-(3.Cc.H(a-x)+(1-f3).Cc.H(x-a) comme décrit précédemment. Dans chaque itération d'index k de l'algorithme, à partir de k = 0, la première étape 21 consiste en une évaluation du résidu:

rk(x) = oik(x)-p(x,îk(x)+î(x)) (4) Ensuite une fonction de mise à jour 8îk(x) est calculée à l'étape 22 en résolvant l'équation de Poisson:

ABîk - div(rk) (5) où A est l'opérateur Laplacien et div l'opérateur divergence. L'équation (5) est résolue à l'étape 22 avec des conditions aux limites de Dirichlet sur 8îk(x1) et 8îk(xN), qui prennent en compte l'amplitude de la discontinuité à la première itération et assurent l'existence d'une solution unique. Ainsi :

8î0(x1) _ -R=(YN-Y1-CJ pour k = 0 (6) { bi0(XN) _ (1-R).(YN-Y1-Cc)

9 8ik(X1) = 0 pour k > 0 (7) 8ik(XN) = 0 La composante ik+1 (élément k+1 de la suite) est calculée à l'étape 23:
ik+1(x) = ik(x) + 8rk(X) (8) Si (test 24) un critère de fin des itérations est rempli (par exemple: nombre d'itérations atteignant une valeur maximale prédéfinie, distance entre les composantes ik+1 et Tk inférieure à un seuil fixe ou adaptatif, etc.), la fonction T
finale est calculée à l'étape 26:

T(x)= Tk(x)+T(x) (9) Dans le cas contraire, l'index d'itération k est incrémenté à l'étape 25 et l'itération 1o suivante est effectuée à partir de l'étape 21 précitée.

Pour des raisons d'implémentation informatique, la résolution des équations différentielles est effectuée dans un domaine discret, du moins pour les abscisses x.

En effet, comme présenté sur la figure 4, les horizons sismiques sont le plus souvent échantillonnés afin de pouvoir être facilement analysés et stockés par des moyens informatiques. Ainsi, la fonction T ne peut pas être considérée comme continue par morceaux. C'est une fonction d'une variable discrète.

Dans ce cas, l'intervalle [x1, xN] est échantillonné sur N points numérotés de à N sur la figure 4. La notation i(xn) désigne la valeur de la fonction T à la n-ième abscisse discrète xn (1 <_ n< N), avec Y1 = T(x1) et YN = r(xN). La discontinuité

réelle Cc, en x = a, n'est pas directement observable sur les données discrétisées. A la place, on peut observer la discontinuité "discrète", ou amplitude de faille dans le domaine discret, notée Cd = i(xna+1) - ti(xna), la position a de la faille étant comprise entre les abscisses discrètes xna+1 et xna.

Un algorithme de pointé adapté à une résolution dans un domaine discret est présenté sur la figure 5. Pour expliquer cet algorithme, on suppose dans un premier temps que les valeurs de na et de Cd sont connues après désignation des points 11, 12 de coordonnées respectives (x1, y,) et (XN, YN) par l'opérateur.
Par rapport à l'algorithme du cas continu présenté en référence à la figure 3, les termes de mises à jour 8ti(xna) et 8i(xna+1) ont a priori des valeurs différentes, 5 ce qui nécessite d'adapter l'algorithme pour le cas discret. Il est notamment nécessaire de mettre à jour la composante de saut Tk de la fonction T au fur et à
mesure des itérations, ainsi que l'estimation Ck de l'amplitude de faille dans le domaine continu, pour tenir compte de cette différence entre 6(x) et bi(xna +1)

10 A l'initialisation 30, où l'index k est pris égal à zéro, l'algorithme fixe la composante pseudo-continue ip , une valeur d'initialisation Co pour l'amplitude de faille dans le domaine continu et la composante de saut ip , par exemple selon:

io(xn) = Yi +P.(YN-Y1) CO = Cd (10) tiO (xn) = -(3.C0.H(n(x+1-n) + (1-P).C0.H(n-na ) avec un paramètre p pouvant être choisi entre 0 et 1. La notation H(n) désigne ici la valeur de la fonction de Heaviside pour un argument entier n (H(n) = 0 si n < 0 et H(n) = 1 sinon).

Dans chaque itération d'index k de l'algorithme, à partir de k = 0, la première étape 31 consiste à nouveau en une évaluation du résidu:

rk(xn)=Vik(xn)-p(xn,ik(xn)+Ik(xn)) (11) Ensuite, la fonction de mise à jour 8ik est calculée à l'étape 32 en résolvant l'équation de Poisson dans le domaine discret:

08ik = -div(rk) (12) avec des conditions aux limites de Dirichlet:

11 8i0(x1) = -R=(YN - Y1 - CO) pour k = 0 (13) 8i0(xN) = (1- R).(YN - Y1 - CO) et 8zk (x1) _ -p.(Ck-1 - Ck) l pour k > 0 (14) 8ik (XN) (1- 0)=(Ck-1 - Ck) le paramètre Ck pour k >- 1 ayant été déterminé à l'itération précédente.

A titre d'exemple, l'équation de Poisson (12) peut être résolue dans un domaine fréquentiel à l'étape 32, en utilisant une méthode par transformée de Fourier discrète (DFT) et transformée de Fourier discrète inverse (DFT-1):

8tik = DFT-1 - DFT[div(rk )] (15) DFT[A]

La transformée de Fourier discrète peut notamment être une transformée en sinus discrète (DST) de type I, auquel cas l'expression (14) devient:

8Tk = DST-1 _ DFT[div(rk)] (16) DST[A]
A l'étape 33 intervient la mise à jour de la composante pseudo-continue selon:

Tk+l(xn)= Tk(xn)+8rk(xn) (17) de l'amplitude de faille dans le domaine continu selon:

Ck+1 = Cd - [Tk+1(xna+1) - ik+1(xna )j (18) et de la composante de saut selon:

'k+l (xn) = -P.Ck+1 =H(na +1-n)+ (1-R).Ck+1.H(n-na) (19) A l'étape 34, l'index d'itération k est incrémenté d'une unité.

Si (test 35) le critère de fin des itérations est rempli (par exemple: nombre d'itérations atteignant une valeur maximale prédéfinie, distance entre les composantes pseudo-continues ik et ik_1 inférieure à un seuil fixe ou adaptatif, etc.), la fonction ti finale est calculée à l'étape 36:

12 T(xn) = ik(xn)+Tk(xn) (20) Dans le cas contraire, l'itération suivante est effectuée à partir de l'étape précitée.

On observe que l'expression (2) de la composante de saut î, dans le cas continu, peut être généralisée sous la forme:

ik (x) = gk (x) + Cc .H(x-a) (2') où gk est une fonction quelconque, de classe C2 sur l'intervalle [x1, xN], choisie arbitrairement pour chaque itération d'index k. L'équation (2) correspond au cas particulier où gk(x) = -(3.Cc pour tout index entier k.

Dans ce cas plus général, l'expression (4) du résidu (étape 21 de la figure 3) devient:

rk (x) = vik (x) + ogk (x) - p(x, tk (x) + I(x)) (4') (le terme ogk est nul lorsque gk est une fonction constante).

Si l'initialisation de la composante continue 'r0 (étape 20 de la figure 3) est également effectuée suivant une fonction f quelconque choisie arbitrairement de classe C2 sur l'intervalle [x1, xN]: i0(x) = f(x), l'expression (6)-(7) des conditions aux limites de Dirichlet est généralisable:

J 810(x1) = Y1 -f(x1)-g0(x1) pour k = 0 (6') 8i0(xN) = (YN -Cc)-f(xN)-g0(xN) 8ik (x1) = 9 k-1(x1) - 9k (x1) 8ik (xN) = gk-1(xN) - gk (xN) pour k > 0 (7') Cette forme plus générale de l'algorithme (avec des fonctions f et gk quelconques) peut également se transposer au cas d'une résolution dans le domaine discret. C'est ce qui est représenté sur la figure 6.

A l'initialisation 40, où l'index k est pris égal à zéro, l'algorithme généralisé de

13 la figure 6 fixe la composante pseudo-continue i0 , la valeur d'initialisation Co pour l'amplitude de faille dans le domaine continu et la composante de saut T0 selon:

ti0(xn)=f(xn) CO = Cd -(f(xna,+1)-f(xna )) (10') Io(x) = g0 (xn) + CO.H(n-na ) Dans chaque itération d'index k de l'algorithme, à partir de k = 0, la première étape 41 consiste en l'évaluation du résidu selon:

rk(xn) =7ik(xn) -vgk(xn)-p(xn,tk(xn)+ik(xn)) (11') Ensuite, la fonction de mise à jour Sik est calculée à l'étape 42 de résolution de l'équation de Poisson dans le domaine discret (12) avec des conditions aux 1o limites de Dirichlet dont l'expression est ici:

Si0(x1)=Y1-f(x1)-90(x1) (13') SiO(xN) = YN -f(xN)-g0(xN)-CO

à la première itération (k = 0), et:

bik(x1) = gk-1(x1)-gk(x1) (14') Sik(xN) = gk-1(xN)-gk(xN)+Ck-1 -Ck aux itérations suivantes (k > 0).

A l'étape 43 intervient la mise à jour de la composante pseudo-continue selon (17) de l'amplitude de faille dans le domaine continu selon (18).

En outre, la composante de saut est mise à jour à l'étape 46 selon:

tik+1(xn) = gk+1(xn) + Ck+1.H(n-na) (19') puis l'index d'itérations k est incrémenté d'une unité à l'étape 47 Si (test 44) l'index d'itération k a atteint sa valeur maximale notée K-1, la fonction i finale est calculée à l'étape 45 à partir de la somme de la composante pseudo-continue et de la composante de saut, par exemple selon:

14 Z(XI) Zk (xn) + Tk (xn) - gk (xn) + gk-1 (xn (20') Le test 44 peut être similaire au test 35 de la figure 5. Une possibilité est de faire porter ce test sur une comparaison de la précédente estimation Ck de l'amplitude de faille avec la nouvelle estimation Ck+1 tenant compte de la pente locale de la fonction ik+1 . Il est également possible de prédéfinir le nombre total d'itérations K.

La figure 5 peut être vue comme un cas particulier la figure 6 dans lequel la fonction f est constante: f(x) = yi + P(yN - Y1), et chaque fonction gk est également constante: gk(x) = -(3.Ck.

Dans la pratique, il n'est pas indispensable de conserver en mémoire toutes les valeurs calculées au cours des itérations successives d'index k. La figure représente ainsi un organigramme d'un algorithme de résolution équivalent à
celui de la figure 6, dans lequel les index d'itérations k n'interviennent plus. La fonction g peut être différente d'une itération à une autre, par exemple g = -P.C.

Dans la suite, on note Tna,Cd la fonction T finalement estimée à l'étape 36 ou lorsqu'on s'est donné les valeurs de la position discrète na de la faille et son amplitude Cd dans un domaine discret.

Les valeurs à retenir pour le couple (na, Cd) peuvent être fournies comme étant supposées connues ou, plus couramment, elles peuvent être inconnues a priori. Dans ce cas, une procédure de recherche du meilleur couple peut être entreprise, par exemple l'une de celles décrites ci après.

Une procédure possible comporte l'évaluation d'un critère de cohérence c(na, Cd) à partir de l'image sismique et de la fonction Tna,Cd obtenue au bout de K itérations. On peut considérer successivement plusieurs candidats pour le couple (n(,, Cd) et retenir celui pour lequel le critère de cohérence s'avère optimal.

En référence à la figure 8, le procédé calcule (étape 302) la fonction i, par exemple selon le procédé présenté par la figure 5, 6 ou 7, sur soumission (étape 301) des valeurs initiales (n(,, Cd).

Une fois la ligne d'horizon sismique définie par cette fonction i, un vecteur u est défini (étape 303) comme étant le vecteur représentant la discontinuité, c'est-5 à-dire u = [(na,t(na)),(na+1,T(na+1))]. En fonction de la discrétisation de l'espace de travail, plusieurs points peuvent être définis le long de ce vecteur. A
titre d'illustration, il est supposé pouvoir être définis L points avec L un entier positif strictement.

En chaque point (parmi ces L points) est calculé (étape 304) un vecteur 1o gradient de l'image sismique complète, avec i e [1..L] entier.

Une matrice de covariance u de ce champ de gradient estimé sur L points est alors calculés (étape 305) comme étant : u = var(X) cov(X,Y) où cov est cov(X,Y) var(Y) l'opérateur covariance et var l'opérateur variance.

Les valeurs propres ? et X2 de cette matrice u sont alors calculées (étape ~ +~

15 306). Enfin, la valeur de cohérence c(na,Cd) est définie comme étant 1 2 (étape 307).

Il est possible dans un mode de réalisation de tester toutes les couples (na, Cd) possibles dans un espace donné, éventuellement discrétisé de valeurs.
Par exemple, il est possible de tester 20 valeurs de Cd reparties de manière homogène dans le segment [0; 21T(xl)-T(xN)I] et avec na compris dans [1..N].

Au final, le couple de valeurs (na, Cd) retenu est celui dont la valeur de cohérence c(n(X, Cd) est minimale parmi tous les couples testés. En effet, si les valeurs propres sont proches, cela signifie en pratique que le vecteur u est positionné sur une discontinuité.

Il est également possible de réaliser le procédé en considérant que Cd est

16 constante et prend par exemple la valeur IT(x1)-z(XN)I. Ainsi, une minimisation de c(na, Cd) est effectuée en faisant varier na tandis que Cd ne varie pas. Puis, une fois la valeur de na permettant de minimiser c(na, Cd) sélectionnée, le procédé
est effectué à nouveau en considérant que na est constant est que Cd varie.

Cette adaptation permet d'optimiser le temps de calcul et le temps de détermination du couple (na, Cd) optimum.

De plus, il est constaté que la convergence vers zéro des conditions incrémentales fCk+1 - Ckl est d'autant plus rapide que l'hypothèse sur Cd est proche de la valeur réelle. Ainsi, une méthode ad hoc de détection des 1o paramètres de discontinuité peut consister, dans un premier temps, à
évaluer le paramètre Cd optimal pour un certain nombre de valeurs possibles de na.

Il est donc possible d'envisager, par exemple, trois hypothèses pour le paramètre Cd centrées sur une valeur a priori Co et régulièrement espacées d'un pas hc. L'algorithme de suivi d'horizon présenté sur la figure 5, 6 ou 7 est alors mis en oeuvre pour un faible nombre d'itérations K (2 ou 3). L'amplitude Cd pour laquelle la différence ICk+1 - Ckl est minimale est alors substituée à Co. Le procédé peut être répété M fois en divisant le pas hc par 2 de sorte à
atteindre une précision satisfaisante.

Bien entendu, la présente invention ne se limite pas aux formes de réalisation décrites ci-avant à titre d'exemples ; elle s'étend à d'autres variantes.

Par exemple, la valeur de cohérence c(na, Cd) peut être calculée de manière différente que celles mentionnées ci-dessus. L'objet de ce calcul est de trouver une décorrélation entre le gradient vertical et le gradient horizontal.

Le procédé qui vient d'être décrit est typiquement mis en oeuvre dans un ordinateur ou une station de travail dont un processeur exécute les étapes ci-dessus sous le contrôle d'un programme dont les instructions sont exécutées par le processeur sur des images sismiques chargées depuis une mémoire de stockage telle qu'un disque dur.

Claims (13)

REVENDICATIONS

1. Procédé de recherche d'un horizon sismique dans une image sismique du sous-sol, le procédé comprenant:
- designer deux points (11, 12) de coordonnées respectives x1, y1 et x N, y N
dans l'image sismique comme étant deux points appartenant à l'horizon recherché, un nombre entier N d'abscisses discrètes x1, x2, ..., x N étant définies entre les abscisses x1, x N des deux points désignés;

- considérer une position de faille à une abscisse discrète x n.alpha. où
n.alpha. est un entier compris entre 1 et N-1 et une amplitude de faille C d dans un domaine discret;
- en se donnant deux fonctions f et g0 dérivables sur l'intervalle [x1, x N], initialiser une amplitude de faille dans un domaine continu à la valeur C0 = C d -[f(x n .alpha.+1) -f(x n .alpha.)], une composante pseudo-continue ~0(x n) et une composante de saut ~0(x n) définies sur les N abscisses discrètes selon ~0(x n)=f(x n) et ~0(x n) = g0 (x n) + C0.H(n-n .alpha.), où H(.) désigne la fonction de Heaviside;
- effectuer plusieurs itérations d'une séquence d'étapes de calcul pour un index d'itération k initialisé à k = 0 puis incrémenté par unités, la séquence comprenant, pour un index k:

-calculer un résidu r k(x n)= .gradient. ~ k(x n)+ .gradient. g k(x n)-p(x n,~k(x n)+~ k(x n)), où .gradient. désigne l'opérateur gradient et p(x n, y) désigne la tangente d'un pendage estimé pour une position d'abscisse discrète x n et d'ordonnée y dans l'image sismique;
- résoudre une équation de Poisson .DELTA.(.delta.~ k )=-div(r k) pour déterminer un terme de mise à jour .delta. ~ k(x n), avec des conditions aux limites de Dirichlet: à la première itération et à chaque itération d'index k >= 1, où g k est une fonction dérivable sur l'intervalle [x1, x N];

- mettre à jour la composante pseudo-continue selon ~k+1(x n) = ~k(x n)+.delta.~k(x n);

- mettre à jour l'amplitude de faille dans le domaine continu selon ;

- mettre à jour la composante de saut selon ~k+1(x n)=g k+1(x n)+C k+1.cndot.H(n-n.alpha.) ;

- si une valeur finale de l'index d'itération k est atteinte, calculer une fonction .tau.n.alpha., C d(x n) représentant l'ordonnée dans l'image sismique d'un horizon estimé en fonction de l'abscisse discrète x n, à partir d'une somme de la composante pseudo-continue et de la composante de saut; et - si la valeur finale de l'index d'itération k n'est pas atteinte, exécuter l'itération suivante de la séquence.

2. Procédé selon la revendication 1, dans lequel la fonction f est choisie comme une fonction constante.

3. Procédé selon la revendication 2, dans lequel la fonction f est prise égale à
Y1 + .beta..(.gamma.N - .gamma.1), .beta. étant un nombre réel choisi entre 0 et 1.

4. Procédé selon la revendication 3, dans lequel chaque fonction g k est choisie comme une fonction constante de valeur -.beta..C k.

5. Procédé selon l'une quelconque des revendications 1 à 3, dans lequel chaque fonction g k est choisie comme une fonction constante.

6. Procédé selon l'une quelconque des revendications précédentes, dans lequel la résolution de l'équation de Poisson utilise une transformée de Fourier discrète DFT et sa transformée de Fourier discrète inverse DFT-1 selon

7. Procédé selon la revendication 6, dans lequel la transformée de Fourier discrète est une transformée en sinus discrète de type I.

8. Procédé selon l'une quelconque des revendications précédentes, comprenant en outre:
- évaluer un critère de cohérence pour les valeurs considérées de l'abscisse discrète n.alpha. et de l'amplitude de faille C d dans le domaine discret.

9. Procédé selon la revendication 8, dans lequel l'évaluation du critère de cohérence comprend:
- calculer les valeurs propres .lambda.1 et .lambda.2 d'une matrice de covariance d'un champ de gradient estimé à partir de l'image sismique en un nombre donné de points du segment [(n.alpha., .tau.n.alpha.C d)), (n.alpha.+1,.tau.n.alpha.,C d(n.alpha.+1))]; et - évaluer le critère de cohérence c(n.alpha., C d) pour les valeurs considérées de l'abscisse discrète n.alpha. et de l'amplitude de faille C d dans le domaine discret selon c(n.alpha., C d) =

10. Procédé selon la revendication 8 ou 9, comprenant en outre:
- identifier un couple de valeurs de l'abscisse discrète n.alpha. et de l'amplitude de faille C d dans le domaine discret pour lequel le critère de cohérence est optimal; et - pointer un horizon défini par la fonction .tau. n .alpha. C d(n) calculée pour le couple donnant le critère de cohérence optimal.

11. Procédé selon l'une des revendications 1 à 7, comprenant en outre:
- évaluer un premier critère de cohérence pour les valeurs considérées de l'abscisse discrète n .alpha., l'amplitude de faille C d dans le domaine discret étant fixe ;
- identifier la valeur de l'abscisse discrète n .alpha. dans le domaine discret pour lequel le premier critère de cohérence est optimal;
- évaluer un deuxième critère de cohérence pour les valeurs considérées de l'amplitude de faille C d dans le domaine discret, la valeur de l'abscisse discrète n .alpha. étant la valeur de n .alpha. identifiée;

- identifier la valeur de l'amplitude de faille C d dans le domaine discret pour lequel le deuxième critère de cohérence est optimal.

12. Procédé selon l'une des revendications 1 à 7, comprenant en outre:
a) évaluer un critère de cohérence pour un ensemble de valeurs parmi les valeurs considérées de l'abscisse discrète n .alpha. et des valeurs de l'amplitude de faille C d dans le domaine discret espacées d'un pas h c et centrées sur une valeur d'amplitude centrale C d0 ;

b) calculer ¦C K - C K-1 ¦ pour un nombre K d'itérations;
c) réitérer M fois les étapes a) à b), - en remplaçant la valeur d'amplitude centrale C d0 par une valeur d'amplitude de faille pour laquelle la valeur de ¦C K+1 - C K¦ est minimale parmi les valeurs de l'amplitude de faille C d évaluée, et - en remplaçant la valeur de distance de h c par une valeur inférieure à la valeur de distance à l'itération précédente.

d) identifier la valeur de l'abscisse discrète n .alpha. parmi les valeurs considérées de l'abscisse discrète pour laquelle le critère de cohérence est optimal, la valeur d'amplitude de faille C d étant une valeur pour laquelle la valeur de ¦C K+1 - C K¦ est minimale parmi les valeurs de l'amplitude de faille C d évaluée lors de l'itération M.

13. Programme d'ordinateur pour un système de traitement d'images sismiques du sous-sol, le programme comprenant des instructions pour mettre en uvre les étapes d'un procédé selon l'une quelconque des revendications précédentes lors d'une exécution du programme par un calculateur du système de traitement d'images sismiques.