PROCEDE 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.
Il 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 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 98085 4 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 5 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 10 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 15 horizon sismique dans une image sismique du sous-sol. Ce procédé comprend: - designer deux points de coordonnées respectives x1, yi 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 xi, xN des deux points désignés; 20 - considérer une position de faille à une abscisse discrète xtIOE 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 go dérivables sur l'intervalle [x-t, xN], initialiser une amplitude de faille dans un domaine continu à la valeur 25 Co = Cd - [f(Xna+i) - f(xn«)], une composante pseudo-continue to (xn ) et une composante de saut î0(xn ) définies sur les N abscisses discrètes selon to(xn) = f(xn) et î 0(xn ) = go(xn)+Co.H(n-na), 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 rk (xn ) = Vtk (xn )+ Vgk (xn ) - p(xn , tk (xn )+-tk (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 etk )= -div(rk ) pour déterminer un terme de mise à jour ôtk(xn), avec des conditions aux limites de Dirichlet: { 80(xi ) = yi -f(x1)-go(xi) à la première itération 5t0(xN)= YN -f(xN)-g0(xN)-00 { ÔTk(x1)- 9k-1(x1)-gk (xi) et à chaque itération 5tk(xN)= 9k_1(xN)-gk(xN)+ Ck_i -Ck d'index k _. 1, où gk est une fonction dérivable sur l'intervalle [x1, xN]; - mettre à jour la composante pseudo-continue selon k-1-1(Xri) = 1.((Xn)-F> ek(Xn); - mettre à jour l'amplitude de faille dans le domaine continu selon Ck+1 =- Cd -K+1(Xnc,+1)- tk+1(Xn« )1 - 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 tn«,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 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 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 à Yi +13.(YN - Yi), R étant un nombre réel choisi entre 0 et 1. Dans ce cas, chaque fonction gk peut 15 aussi être prise constante, en particulier de valeur -13.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 &Ci( ---- DFT-1 - DFT[0] - D FT[div(rk)] . En particulier, on peut utiliser une transformée en sinus discrète de type I. 20 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. 25 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 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 -Io 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 15 é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. 20 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 25 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. 30 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 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 = a. 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 1 o 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. 15 L'opérateur (interprétateur d'images sismiques) est typiquement sollicité pour indiquer deux points 11, 12 de coordonnées respectives (x1, y1) 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 20 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 25 T(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(xi) = yi et t(xN) = yN. La fonction t présente un saut d'amplitude Cc = t+ - t_ à l'abscisse a E ]x1, xN[, où t+ = lim t(x) est la valeur limite de la fonction t lorsque x tend vers a par x -m+ valeurs supérieures et T__ = lim t(x) est la valeur limite de la fonction t lorsque x x--->cc- 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 t compte tenu des données de l'image sismique. Ces données peuvent être prétraitées pour estimer un pendage 0(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 Dxs = (u, v) suivant laquelle les valeurs de l'image sismique présentent un gradient d'amplitude maximum et à représenter le pendage 0(x, y) par sa tangente p(x, y), c'est-à-dire l'opposé de la cotangente de l'angle formé par la direction Dxs 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 15 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 '1 et la tangente p du pendage local sont liés par l'équation suivante Vt(x) = p(x, t(x)) (1) 20 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 25 fait de la discontinuité, il est proposé de décomposer la fonction t qui est de classe C2 sur [x1, xN] - {a} en la somme d'une composante continue inconnue tr de classe C2 sur [x1, xN] vérifiant l'équation (1) et d'une composante de saut discontinue î constante sur [x1, a[ et constante également sur ]a, xN]: Î(X) = -p.Cc.H(a-x)+ (1-(3).Cc.H(x-a) (2) H étant la fonction de Heaviside (H(x) = 0 si x 0 et H(x) = 1 sinon), et [3 é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 5 3. L'algorithme a pour objet de calculer une suite de composantes continues tk 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 to de la suite, par exemple à une valeur constante: 10 tb(x) = Y1 ± (3.(YN-Y1) (3) prise entre les ordonnées yi et yN des points 11, 12 si le paramètre (3 est choisi entre 0 et 1. La fonction 't(x) est alors fixée à la valeur t(x) = -13.Cc.H(a-x)+ (1-p).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 15 consiste en une évaluation du résidu: rk (x) = Vtk (x) - p(x, Ti( (x) + î(X)) (4) Ensuite une fonction de mise à jour Stk(x) est calculée à l'étape 22 en résolvant l'équation de Poisson: Llek - div(rk ) (5) 20 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 ôtk(xi) et ôtk(xN), qui prennent en compte l'amplitude de la discontinuité à la première itération et assurent l'existence d'une solution unique. Ainsi : 3t0(x1) = -[3.(YN-Y1-Cc pour k = 0 8To(xN)= (1 -(3).(YN-Yi-Cc) (6) Mk(xi)= 0 mk(xN) = 0 pour k > 0 (7) La composante tk+1 (élément k+1 de la suite) est calculée à l'étape 23: Tkil(x) = tk(x)+ ek(x) (8) Si (test 24) un critère de fin des itérations est rempli (par exemple: nombre 5 d'itérations atteignant une valeur maximale prédéfinie, distance entre les composantes t'k+1 et tk inférieure à un seuil fixe ou adaptatif, etc.), la fonction finale est calculée à l'étape 26: t(x) = (x)+ t(x) (9) Dans le cas contraire, l'index d'itération k est incrémenté à l'étape 25 et l'itération 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 15 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, xiy] est échantillonné sur N points numérotés de 1 à N sur la figure 4. La notation t(xn) désigne la valeur de la fonction 'I' à la n-ième 20 abscisse discrète xn (1 .5 n 5_ N), avec yi = t(x1) et yN = t(xN). La discontinuité réelle Cc, en x = a, n'est pas directement observable sur les données discrétisées. À la place, on peut observer la discontinuité "discrète", ou amplitude de faille dans le domaine discret, notée Cd t(Xdo), la position a de la faille étant comprise entre les abscisses discrètes xr4+1 et xna. 25 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, yi) 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 St(xn" ) et 5T(xnu+i) ont a priori des valeurs différentes, 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 Ôt(xn« ) et ôt(xn,+1).
A l'initialisation 30, où l'index k est pris égal à zéro, l'algorithme fixe la composante pseudo-continue to , une valeur d'initialisation Co pour l'amplitude de faille dans le domaine continu et la composante de saut 'Co , par exemple selon: { to(xn) = y, + r3.(yN-Y1) to(xn) ---- ---P.co.H(na+i-n)+ (1-(3).co.F1(n-na) Co = Cd (10) avec un paramètre E3 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 ) = Vtk (xn ) - p(xn, tk (xn )+tk (xn )) (11) Ensuite, la fonction de mise à jour 5% est calculée à l'étape 32 en résolvant l'équation de Poisson dans le domaine discret: AfStk -- . -div(rk ) (12) avec des conditions aux limites de Dirichlet: {{8Yo(xi) - 13.(YN - Yi -00) Sto(xN) = (1-13).(YN - Yi -Co) ôtk (x1) = -E3.(Ck-1 - Ck) Ck ) pour k > 0 3t1( (XN ) = (1- (3).(Ck--1 - 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 5 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): Stk =DFT-1 DFT[cliv(rk )1 DFT[3,1 _ (15) 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: 10 Ôtk = DST-1 DST[div(rk )1 _ DST[3.1 _ (16) À l'étape 33 intervient la mise à jour de la composante pseudo-continue selon: Tk-o(xn)=--- tk(xn)+Stk(xn) (17) de l'amplitude de faille dans le domaine continu selon: Ck+1 = Cd - K+1(Xn,+1)- 'Ic-F1(Xna )1 (18) 15 et de la composante de saut selon: ek+1 ( xn ) = -idek-o .1"i(na +1-n)+ (1-eCk+1-11(n-na.) (19) À 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 20 composantes pseudo-continues et tk-1 inférieure à un seuil fixe ou adaptatif, etc.), la fonction t finale est calculée à l'étape 36: pour k = 0 (13) et (14) T(Xn)= tk(Xn)+Îk(Xn) (20) Dans le cas contraire, l'itération suivante est effectuée à partir de l'étape 31 précitée. On observe que l'expression (2) de la composante de saut î, dans le cas 5 continu, peut être généralisée sous la forme: tk (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) = Vt'k(x)+Vgk(x)-p(x,k(x)+ 't(x)) (4') (le terme Vgk est nul lorsque gk est une fonction constante). Si l'initialisation de la composante continue (étape 20 de la figure 3) est 15 également effectuée suivant une fonction f quelconque choisie arbitrairement de classe C2 sur l'intervalle [x1, xN]: to(x) = f(x) , l'expression (6)-(7) des conditions aux limites de Dirichlet est généralisable: .f âto(xi)-- Y1 -f(x1)-go(x1) isto(xN)- (YN -cc)-f(xN)-go(xN) pour k = 0 (6') Ôtk(x1)=gk-1(x1)--9k(x1) (7') Stk(xN)= gk--i(xN)-gk(xN) Pour k > 0 20 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 la figure 6 fixe la composante pseudo-continue to , la valeur d'initialisation Co pour l'amplitude de faille dans le domaine continu et la composante de saut to selon: to(xn)= f(xn) (10') Co = Cd - (f(Xn +1)- f(xna to(xn)= 9o(xnc;+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 ) Vic(xn ) Vgk (Xn )-P(Xn, tk (Xn )-1-tk(xn Ensuite, la fonction de mise à jour 5% est calculée à l'étape 42 de résolution de l'équation de Poisson dans le domaine discret (12) avec des conditions aux limites de Dirichlet dont l'expression est ici: 3to(x1)=y1 -f(x1)-go(x1) 5T0(xN)= yN -f(xN)-go(xN)-co à la première itération (k = 0), et: -3k (X1) = 9k--1(X1)- gk (X1 ) Ôtk(X1\1) = gk-1(XN)- gk(XN)+ Ck-1 -Ck aux itérations suivantes (k > 0). À 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: tk+1(xn ) = gk+1(xn ) Ck+1.11(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 1r 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: (13') (14') t(x,)= (x,)+ tk(x,,)- 9k(x' )-F 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 tk+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) = y1 + vR(yN - et chaque fonction gk est également constante: gk(x) = 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 7 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 = Dans la suite, on note tnacd la fonction finalement estimée à l'étape 36 ou 45 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 (ne, Cd) peuvent être fournies comme étant supposées connues ou, plus couramment, elles peuvent être inconnues a 20 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 tria,cd obtenue au bout de K itérations. On peut considérer successivement plusieurs candidats pour le 25 couple (ne, 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 'r, par exemple selon le procédé présenté par la figure 5, 6 ou 7, sur soumission (étape 301) des valeurs initiales (na, Cd). Une fois la ligne d'horizon sismique définie par cette fonction 'r, un vecteur u est défini (étape 303) comme étant le vecteur représentant la discontinuité, c'est- à-dire u = [(na, t(na )),(na+1,t(na +1))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 gradient (xi,yi) de l'image sismique complète, avec i E [1. Li 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)cov(X,Y) var(Y) où cov est l'opérateur covariance et var l'opérateur variance. Les valeurs propres Xi et X de cette matrice y sont alors calculées (étape 'Xi -x2i x1+ x2 (é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 20 homogène dans le segment [0; 21t(x1)--t(xN)II et avec na compris dans [1..M. Au final, le couple de valeurs (na, Cd) retenu est celui dont la valeur de cohérence c(na, 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é. 25 Il est également possible de réaliser le procédé en considérant que Cd est 306). Enfin, la valeur de cohérence c(na,Cd) est définie comme étant constante et prend par exemple la valeur ft(x1)---c(xN)1. 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 ICk+i - Cki 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 o 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 15 mis en oeuvre pour un faible nombre d'itérations K (2 ou 3). L'amplitude Cd pour laquelle la différence ICk+i - 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 20 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 25 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.