FR2616920A1 - Inversion d'un profil sismique vertical en minimisant une fonction du type entropie - Google Patents

Abstract

L'invention concerne un procédé d'inversion de traces sismiques Sk en déterminant la série R{ri } représentative des coefficients de réflexion sismique des couches géologiques. Le procédé consiste : - à déterminer une solution approchée Rdegre(s) en minimisant une fonction H(R)= SIGMAf(ri ) du type entropie sous des contraintes liées aux champs d'ondes remontantes Uk et descendantes Dk, - à déterminer une valeur corrective DELTAR en minimisant une fonction linéaire J(DELTAR) sous des contraintes linéaires liées au bruit contenu dans les traces et à la courbe temps-profondeur extraite des traces, et - à calculer R=Rdegre(s)+DELTAR représentative des coefficients de réflexion. Selon une variante, il est possible d'orienter le choix de la solution approchée Rdegre(s) en englobant des informations contenues dans une ou plusieurs diagraphies, et de déterminer DELTAR pour que la série R respecte les vitesses d'addition qui résultent du traitement de la sismique de surface.

Description

4-

INVERSION D'UN PROFIL SISMIQUE VERTICAL EN MINIHISANT

UNB FONCTION DU TYPE ENTROPIE

L'invention concerne l'étude du Profil Sismique Vertical généralement connu sous le sigle PSV, et plus particulièrement

l'inversion d'un tel profil.

D'une manière générale, un signal sismique-recueilli par un détecteur est la réponse de la variation d'impédance acoustique des formations géologiques à l'émission d'une onde sismique. L'impédance acoustique d'une formation est définie comme étant le produit de la densité du milieu par la vitesse de propagation de l'onde sismique à

travers ce milieu.

L'acquisition d'un PSV consiste à enregistrer les signaux sismiques recueillis à différents niveaux d'un sondage en réponse à des excitations sismiques émises à la surface du globe terrestre. L'on peut également imaginer la configuration inverse, à savoir des détecteurs situés en surface du globe en combinaison avec des excitations

sismiques à différents niveaux d'un sondage.

L'étude du PSV consiste notamment à résoudre un problème d'inversion des signaux sismiques de manière à en extraire l'impédance acoustique de la formation et ainsi reconstruire le profil d'impédance

acoustique en fonction d'une échelle de temps et/ou de profondeur.

-2- Le profil d'impédance acoustique recherché correspond mathématiquement à une suite de dioptres dont il faut définir le nombre n, les positions géométriques et les valeurs d'impédance acoustique tAi, i=l,n. D'une manière équivalente, il est défini par une suite de coefficients de réflexion tri, izl,n/ placés aux interfaces des dioptres et dont la valeur est donnée par la formule: ri = ( Ai+l - Ai) / ( Ai+, + Ai) Une difficulté majeure rencontrée lors de l'inversion réside dans le fait que le nombre d'inconnues (par exemple les positions et valeurs des coefficients de réflexion) est très grand (typiquement de l'ordre du millier), car chaque échantillon représentant une profondeur est un candidat potentiel pour être l'emplacement d'un coefficient de réflexion, tandis que le nombre d'équations disponibles dépend de la largeur du spectre en fréquence des traces sismiques et est, en

général, très inférieur au nombre d'inconnues.

Parmi les techniques d'inversion connues, l'on soulignera celle décrite par A. Bamberger et al dans l'article intitulé "Inversion of normal incidence seismograms" et paru dans Geophysics (Vol. 47, mai 1982, pages 757-770). Cette technique est considérée comme l'une des plus satisfaisantes d'un point de vue strictement mathématique, mais sa mise en oeuvre présente néanmoins l'inconvénient d'être extrêmement

coûteuse en temps de calcul.

Une autre technique connue consiste, dans une première étape, à selectionner un nombre limité de positions possibles pour les interfaces, puis, dans une seconde étape, à calculer, pour chacune de ces positions déterminées, la valeur du coefficient de réflexion de façon à minimiser l'écart entre la trace enregistrée et la trace synthéthique obtenue à l'aide du profil trouvé. L'article intitulé "Inversion of Vertical Seismic Profile by Iterative Modelling" de P. Grivelet et paru dans Geophysics (Vol. 50, juin 1985, pages 924-930) -3- résume bien les différentes étapes d'un tel procédé. Cet article souligne d'autre part qu'il est souhaitable de réduire le nombre des solutions possibles en imposant des conditions-qui résultent de la connaissance a priori de la valeur de l'impédance acoustique à travers une ou plusieurs couches géologiques déterminées. I1 faut tout d'abord souligner que le procédé évoqué ci-dessus ne permet, dans sa présentation actuelle, d'inverser qu'une trace à la fois, limitation qui peut se révéler très gênante lorsque le signal faisant l'objet de l'inversion est affecté par du bruit. En effet, il est bien connu que le bruit risque d'introduire des interfaces

fictives, donc des couches géologiques inexistantes.

D'autre part, la seconde étape ne permet pas de redéfinir de nouvelles positions d'interface en cas de désaccord avec les équations sismiques, et reste donc prisonnière d'erreurs éventuelles qui peuvent

résulter de la première étape.

Enfin, ce type d'inversion accepte comme seule information additionnelle la connaissance a priori de l'impédance acoustique de

certaines couches géologiques.

En résumé, les techniques connues ne permettent pas d'optimiser le procédé d'inversion en y injectant des informations contenues dans tous les signaux sismiques recueillis, ainsi que celles connues préalablement et dont la nature est différente de celle de la

solution recherchée.

La présente invention vise un processus d'inversion qui permet d'englober simultanément des informations disponibles à propos des coefficients de réflexion. Une première source d'informations disponibles est constituée par des informations proprement acoustiques, notamment la mesure du temps de propagation des ondes entre les différents niveaux d'acquisition du PSV, les estimations du bruit -4- affectant les traces enregistrées, la connaissance éventuelle des vitesses d'addition fournies par le traitement de la sismique de surface. Une seconde source importante d'informations se trouve dans les mesures de diagraphie, telles que les mesures de résistivité, de porosité, de densité, etc..., qui ont l'avantage par rapport aux données sismiques d'être échantillonnées selon un pas beaucoup plus petit. L'étude de ces diagraphies permet de définir des positions probables d'interfaces majeures avec une très bonne précision selon la

direction longitudinale du puits.

Elle vise également à minimiser, autant que possible, le temps de calcul nécessaire à l'unité centrale d'un ordinateur pour effectuer l'inversion du PSV, et tout particulièrement de maintenir le temps sous

un seuil commercialement acceptable.

Selon un premier aspect de l'invention, le procédé consiste à identifier les couches géologiques souterrainnes à partir d'une pluralité de traces sismiques Sk(t) enregistrées dans un sondage à différents niveaux k de profondeur Zk, chaque trace Sk(t) étant composée de n valeurs discrètes échantillonnées en fonction du temps t selon un pas constant Et, en déterminant la série R = tri, i = 1,n1 représentative des coefficients de réflexion sismique des interfaces formées par les couches géologiques, i étant un indice relatif à la

profondeur échantillonnée en échelle temps selon ledit pas constant At.

Le procédé comprend les étapes suivantes: - établir la courbe tempsprofondeur reliant, pour chaque niveau k, la profondeur Zk au temps de transit Tk relevé sur la trace Sk(t); - extraire des traces Sk(t), le champ d'ondes montantes Uk(t) et le champ d'ondes descendantes Dk(t); calculer UmoY(t), umin(t) et Umax(t) représentatives respectivement des valeurs moyenne, minimale et maximale du champ d'ondes montantes Uk(t-Tk), et Dm Y(t) représentatif de la valeur

moyenne du champ d'ondes descendantes Dk(t+Tk).

26 16920

-5- Le procédé est caractérisé par le fait qu'il comprend en outre: a) dans une première étape, déterminer une série R ={r i= représentative d'une valeur approchée de la série R en minimisant une fonction H(R)= E f(ri) dans laquelle f(ri) est une fonction croissante de Iril, symétrique en abscisse et concave, sous la contrainte Dmoy * R = Umoy dans laquelle l'opérateur * représente un opérateur de convolution; b) dans une seconde étape, déterminer une série àR = {6ri/ représentative d'une valeur corrective de la série R en minimisant une fonction J(AR) linéaire en I5ri] sous les contraintes: U mn(t) < [(R + àR) * Dm Y(t)] < umax(t) et IZk - Zks(R +àR)l < ck dans laquelle Zks(R +âR) est une fonction linéarisée en 6R, représentative de la profondeur Zks calculée d'une manière synthétique au niveau k à partir de la solution initiale R , et la grandeur úk étant représentative de l'incertitude sur les grandeurs de Zk et Tk; c) calculer la série R = R +&R tri = r i+6ri} représentative

des coefficients de réflexion sismique.

De préférence, les valeurs Um Y(t) et Dm Y(t) sont tranformées dans le domaine des fréquences de manière à obtenir les valeurs Um Y(w) et Dm Y(w) et les contraintes de la première étape sont exprimées dans le domaine des fréquences, à savoir: F. R = UmoY(W)/ Dm Y() dans laquelle F est la matrice associée à la transformée de Fourier,

restreinte à la bande de fréquence utile des traces sismiques.

De préférence, la fonction J(&R) est de la forme

J(aR) = E aienri.

dans laquelle a. est-un coefficient de pondération dépendant de r10.

-6- Une variante du procédé selon l'invention comprend en outre les étapes suivantes: - - à l'aide de ladite courbe temps-profondeur Zk, Tkt, transposer au moins une courbe de diagraphie Y(d) enregistrée préalablement en fonction de la profondeur d dans ledit sondage sur un intervalle de profondeur qui recouvre au moins les profondeurs Zk, en une diagraphie discrète Y(t). iyi(t) échantillonnée en échelle temps selon le même pas àt; - calculer la différentielle discrète Y'(t)= {y'i(t)! de Y(t) , ainsi que sa transformée dans le domaine des fréquences Y'(X); et la première étape du procédé selon l'invention consiste à déterminer deux séries R =/r iQ et Y'=/y'ii respectivement représentatives d'une solution initiale approchée et de la différentielle discrète de la A diagraphie zonée Y(t) et pour lesquelles les valeurs nulles de r i et Y'i portent le même indice i, en minimisant simultanément les fonctions H(R) E f(r.) A H(Y') > e(y') dans laquelle e(y'i) est une fonction croissante de y'il, symétrique en abscisse et concave, sous les contraintes F.R = Umo Y() /DmoY()

E.Y' Y'(w).

dans laquelle E est la matrice associée à la transformée de Fourier, restreinte à un nombre identique de fréquences à celui de la bande de

fréquence utile des traces sismiques.

Une autre variante du procédé selon l'invention consiste à minimiser la fonction J(MR) sous une contrainte additionnelle Va - Vas(R +àR) < cv dans laquelle Va sont les vitesses d'addition correspondant aux temps sa obtenues suite au traitement de la sismique de surface, et Vas(RO+àR) est une fonction linéarisée en àR, représentative des vitesses d'addition Vas calculées d'une manière synthétique à partir de la solution initiale R , et la grandeur Ev étant représentative de

l'incertitude sur les vitesses d'addition Va.

-7- Selon un second aspect de l'invention, le procédé consiste à identifier les couches géologiques successives traversées par un A sondage en effectuant le zonage d'une diagraphie discrète Y(t) yi(t), i=l, n échantillonnée selon un pas constant at, de manière à en extraire la série Y' y'i, i. 1, n} indicative de la position des interfaces formées par les couches géologiques successives. Le procédé comprend les étapes suivantes: - calculer la différentielle discrète Y',. { i de Y, ainsi que sa transformée dans le domaine des fréquences Y'(X); A - déterminer la série Y' indicative de la position des interfaces en minimisant une fonction AA H(Y')- I e(y'.) dans laquelle e(y'i) est une fonction croissante de ly'il, symétrique en abscisse et concave, sous la contrainte A E.Y' = Y'(w) dans laquelle E est la matrice associée à la transformée de Fourier,

restreinte à un certain nombre de fréquences sélectionnées.

De préférence, les fonctions e et f sont identiques: Iril f(r.) (1 - e) i Ayf 1 -Iy'il e(Y'i) - (1 - e) D'autres caractéristiques et avantages de l'invention

apparaîtront à l'examen de la description détaillée ci-après, ainsi que

des dessins annexés, sur lesquels: - la figure 1 représente un ordinogramme général du procédé selon l'invention qui consiste à identifier les couches géologiques successives traversées par un sondage, uniquement à partir des informations contenues dans les traces Sk(t) de la sismique de fond; -8- - la figure 2 représente un ordinogramme détaillé illustrant un moyen particulier pour accomplir la minimisation de la première étape qui a été décrite en référence à la Fig. 1; - la figure 3 représente un ordinogramme général d'une variante du procédé selon l'invention qui consiste à identifier les couches géologiques successives traversées par un sondage, non seulement à partir des informations contenues dans les traces Sk(t) de la sismique de fond, mais également à partir d'informations provenant de la sismique de surface et/ou d'une ou plusieurs diagraphies; - la figure 4 représente un ordinogramme détaillé illustrant un moyen particulier pour accomplir la minimisation de la première étape qui a été décrite en référence à la Fig. 3; - la figure 5 représente, à l'aide d'un exemple illustratif, les traces Sk(t), les champs d'ondes montantes Uk(t-Tk) et descendantes Dk(t+ Tk) déconvoluées, et ies traces moyennes montante Um Y(t) et descendante Dm Y(t); - les figures 6A, 6B, 6C et 6D représentent respectivement, à l'aide du même exemple illustratif, - la série de coefficients de réflexion ri obtenue i selOn le procédé décrit en référence à la Fig. 1; - la série de coefficients de réflexion r. obtenue selon la variante décrit en référence à la Fig. 3; - le profil d'impédance acoustique obtenu selon la variante décrit en référence à la Fig. 3; - la diagraphie échantillonnée Y(t) et la diagraphie

Y(t) obtenue après l'opération de zonage.

Pour résoudre, d'une manière simple et rapide, le problème complexe que pose l'inversion, l'on fait habituellement les hypothèses suivantes: d'une part, uniquement les phénomènes acoustiques sont effectivement pris en compte, - d'autre part, toutes les couches géologiques forment des plans horizontaux, ce qui conduit à résoudre un problème à une dimension, - enfin que la récupération de gain effectuée compense les

effets de divergence sphérique.

Ces hypothèses conduisent à formuler l'équation suivante: Uk = R *Dk + Bk (I) dans laquelle - Uk représente l'onde montante enregistrée au niveau k; - Dk représente l'onde descendante enregistrée au niveau k; - R représente les coefficients de réflexion exprimés en échelle temps; - Bk représente le bruit défini au sens mathématique comme étant tout ce qui ne correspond pas à une onde réfléchie primaire, c'est à dire qui ne contient pas de réflexions multiples; et

- le signe * est un opérateur de convolution.

Il s'agit donc de résoudre l'inconnue R à partir des ondes

remontante Uk et descendante Dk.

Dans le domaine des fréquences, l'équation I devient: Uk(w) = R(w). Dk(w) + Nk(w) (II) dans laquelle - le terme Nk(w) est le bruit défini au sens mathématique comme étant tout ce qui ne correspond pas à une réflexion primaire

telle que définie par l'équation des ondes à une dimension.

- 10- L'ensemble des solutions pour l'équation (II) est infini. Par conséquent, il est souhaitable de limiter le domaine des solutions possibles à l'aide de contraintes qui résultent de la connaissance d'informations provenant d'origines diverses et dont la nature n'est pas forcément la même que celle de la solution recherchée. D'après la théorie de l'entropie, aussi appelée théorie de l'information, une solution particulière et préférable dans le domaine

des solutions possibles est celle qui présente un minimum d'entropie.

Cependant, étant donné que les variables manipulées ne sont pas statistiques, on considérera une notion élargie du concept d'entropie qui en retient les propriétés essentielles, notamment la croissance de

l'entropie en fonction de la dispersion des variables.

Les travaux entrepris par la demanderesse l'ont conduit à définir d'une manière générale des fonctions du type entropie qui sont croissantes en valeur absolue, symétriques en abscisse et concaves au sens large; les fonctions concaves au sens large incluent notamment les fonctions linéaires en valeur absolue, telle que celle proposée dans l'article intitulé "Reconstruction of a sparse spike train from a portion of its spectrum and application to high resolution deconvolution" de S. Levy et P.K. Fullagar et paru dans Geophysics (Volume 46, pages 1235-1243) n B (X) E x I i i=l Parmi les fonctions du type entropie qui sont croissantes en valeur absolue, symétriques en abscisse et concaves au sens large, la demanderesse préfère tout particulièrement la fonction strictement concave suivante: -11- n B(R) = ( 1 - e i) = f(ri) i=l qui présente également l'avantage d'être nulle à l'origine. Procédé d'inversion ( Fig. 1) La Figure 1 représente schématiquement l'ordinogramme général du procédé selon l'invention. Les traces sismiques Sk(t) représentées sur la Fig. 5 ont été enregistrées dans un puits à différents niveaux k de profondeur Zk et sont échantillonnées selon un pas constant At. Les traces Sk(t) (1) servent à établir, d'une manière connue, la courbe temps- profondeur (2) sous la forme d'une série de couples de valeurs Zk,Tk} dans laquelle Tk représente le temps de transit ou d'arrivée de

la première onde relevé sur la trace Sk(t) correspondante.

Les traces Sk(t) subissent ensuite un traitement préalable (3) qui consiste à extraire du champ de traces Sk(t) les champs d'ondes montantes Uk(t) et descendantes Dk(t) à l'aide d'un filtre de vitesse connu, tel que celui décrit dans le brevet EP-0,053,525 de Seeman & Horowicz par exemple. Il est souhaitable d'appliquer aux champs un opérateur de déconvolution classique de manière à obtenir des champs d'ondes montantes Uk(t) et descendantes Dk(t) qui ne contiennent que

des ondes primaires.

La Fig. 5 représente également les champs d'ondes montantes Uk(t-Tk) et Dk(t+Tk) respectivement décalés vers la droite et vers la gauche de la valeur de Tk, de façon à respectivement aligner les

événements remontants et descendants.

A partir des champs d'ondes montantes Uk(t-Tk) et descendantes Dk(t+Tk) ainsi obtenus, on sélectionne une fenêtre correspondant par exemple à un intervalle de 100 ms, pour déterminer (4), d'une manière connue, les valeurs suivantes: -12- - Um Y(t) représentative d'une valeur moyenne de l'onde montante dans le domaine des temps, et illustré en Fig. 5, Umoy(w) représentative d'une valeur moyenne de l'onde montante dans le domaine des fréquences; - Dm Y(t) représentative d'une valeur moyenne de l'onde descendante dans le domaine des temps et illustré en Fig. 5; Dm Y(w) représentative d'une valeur moyenne de l'onde descendante dans le domaine des fréquences; - Umax(t) représentative d'une valeur maximale de l'onde montante dans le domaine des temps; min - Umin(t) représentative d'une valeur minimale de l'onde

montante dans le domaine des temps.

Il est à noter que les quantités "valeur moyenne, valeur minimale et valeur maximale" peuvent aussi bien être déterminées arithmétiquement qu'être déduites d'une analyse statistique; par exemple on peut prendre pour valeurs extrêmes des quantités déduites de

la valeur moyenne et de l'écart-type calculé.

Les valeurs numériques utilisées dans l'exemple ci-après ne sont données qu'à titre purement illustratif et, par conséquent, ne peuvent donc pas être interprétées comme une limitation quant à la portée de l'invention. Pour l'exemple traité, le pas d'échantillonage

At des traces Sk(t) est de lms sur un intervalle de 256 ms.

Dans une première étape (5) dans laquelle le bruit N n'est pas pris en compte, c'est à dire qu'il est considéré comme- étant nul, il s'agit de déterminer une estimation initiale R =/r ii (6) de la solution à l'intérieur du domaine des solutions en minimisant une fonction du type entropie, croissante en valeur absolue, symétrique en abscisse et concave au sens large, et préférablement la fonction H(R) suivante: -13- n _Iri I H(R) = ( 1 - e > = t f(ri) i=l avec n = 256 pour l'exemple considéré, tout en vérifiant les contraintes linéaires qui résultent de l'équation (II): Um Y(t) = R * Dm Y(t) ou encore dans le domaine des fréquences: F. R = Um Y(w)/DmoY(X) dans laquelle - F est une matrice (256x256) associée à la transformée de Fourier rapide "FFT" de Cooley-Tuckey et restreinte à la bande de fréquence utile; - R une matrice (lx256) composée de la série <ri, i étant un indice variant de 1 à- 256 et relatif à la profondeur échantillonée selon une échelle temps à intervalle régulier At; - Um Y(o) /Dm Y(w) une matrice (lx256) composée des 256 échantillons u m Y/dimoy. Les contraintes représentent donc un ensemble i i de 256 équations linéaires pour chaque pulsation X ou fréquence f

entière correspondante.

Il faut à ce propos souligner que cette contrainte n'est significative que pour la bande de fréquence [fmin, fmax] effectivement utilisable, cette dernière étant liée à la nature et à la qualité de la source sismique même. Généralement, l'on considère qu'aux sources sismiques couramment utilisées correspond une bande de fréquence utile

située entre 10 et 80 Hz. Pour l'exemple utilisé dans la description,

nous prendrons à titre illustratif une bande de fréquence utile située entre 10 et 50 Hz; pour les fréquences inférieures à 10 Hz et supérieures à 50 Hz, l'on considérera que D(X) = 0, c'est-à-dire du

bruit.

-14- Par conséquent, seulement un nombre m d'équations de contraintes sont à considérer, ce qui signifie que seuls m coefficients peuvent être calculés de façon unique, m étant défini par la formule n (fmax - fmin + 1) 256 (50-10+1) m =-=-- 21 < n fnyq 500 dans laquelle - n représente le nombre d'échantillons pour chacune des traces sismiques, soit 256 dans l'exemple décrit; - fmax et fmin définissent la bande de fréquence utile, soit à 50 Hz dans l'exemple décrit; - fnyq représente la fréquence de Nyquist dérivée du pas d'échantillonnage des signaux sismiques, soit 500 Hz dans l'exemple décrit, soit m = 21 interfaces possibles pour un intervalle analysé de 256 ms, ou encore sous une autre forme, un maximum de 21 coefficients de

réflexion r. non nuls.

Le processus de minimisation (5) de la fonction H(R) peut s'effectuer par toute méthode primale travaillant dans le domaine faisable ainsi que démontré dans le livre intitulé "Linear and Nonlinear Programming" de D.G. Luenberger, 1984 paru chez "Addison-Vesley publishing Company". En ce qui concerne les méthodes primales, l'on décrira, d'une manière détaillée en référence à la Fig. 2, un processus de minimisation à l'aide d'un algorithme de gradient réduit modifié, étant entendu que ce processus pourrait également s'effectuer, d'une manière similaire, à l'aide d'un algorithme de

gradient projeté modifié.

-15- Le résultat de la minimisation conduit à une solution initiale (6) approchée R0 = <rOi avec i= 1 à 256, dont 21 coefficients de

réflexion non nuls, les 235 autres coefficients étant nuls.

Nous abordons maintenant la seconde étape du procédé selon l'invention dans laquelle l'on affine la solution initiale R en posant l'équation suivante: R = R + AR ou ri = r i + Sri avec i = 1, 256 Il s'agit donc de calculer l'incrément AR = /6ri en prenant en compte le bruit Bk et en respectant les informations relatives à la courbe temps-profondeur >Zk, Tk. A cet effet, l'on minimise (7) la fonction linéaire pondérée J(àR) suivante: n J(AR)= Lai I{ri} avec AR = 6ri} i=! dans laquelle ai représente un facteur de pondération qui pénalise les échantillons pour lesquels r i est nul en leur affectant un poids plus

fort que celui affecté aux échantillons pour lesquels r i est non nul.

Plus particulièrement, l'on préférera pour ai un nombre réel positif correspondant à la plus grande des deux valeurs suivantes l/r min,

1/rI..

i min a max / 1/r mIn /r

avec r min = valeur minimale des coefficients r i non nuls.

Cette minimisation est soumise aux contraintes définies, d'une

part par le bruit et, d'autre part par la courbe temps-profondeur.

26 16920

-16- D'une part, les contraintes liées au bruit découlent des valeurs maximales et minimales Umax(t) et Umin(t) déterminées précédemment (4) qui définissent une plage à l'intérieur de laquelle la solution finale R doit se trouver. Autrement dit, ces contraintes peuvent se mettre sous la forme de l'inégalité suivante linéaire en AIR umin(t) < ( R + R) * Dm Y(t) < U a (t) D'autre part, les contraintes liées à la courbe temps- profondeur découlent de la série de valeurs >Tk, Zk} déterminées précédemment (2). En effet, il est possible de poser que la différence entre la profondeur réelle Zk et la profondeur synthétique Zks, calculée d'une manière indépendante à partir des vitesses de propagation, doit être inférieure à une grandeur ek qui se déduit de l'incertitude sur les erreurs de mesure sur Zk et Tk. Cette contrainte

peut s'exprimer mathématiquement sous la forme de l'inégalité sui-

vante: IZk - Zks (Ro+AR) < ck Sous une forme développée, cette inégalité devient: Nk Zk - et Vi < ck i = 1 dans laquelle At représente le pas d'échantillonnage, soit lms dans l'exemple traité, Nk représente le temps de la première arrivée à la profondeur Zk exprimée en nombre d'échantillons, et Vi la vitesse de propagation de l'onde sismique dans

l'intervalle [ i, i+l de 1 ms.

l'intervalle [ i, i+l I de 1 ns.

-17- Si l'on considère que la densité des formations est constante, l'on peut écrire: Nk 1+r. |Zk - Vl.àt < ck

1 - r.

i=1 j=1 3 dans laquelle V1 est la vitesse de la couche supérieure déterminée à

partir de la trace enregistrée au niveau le moins profond.

En différenciant l'équation ci-dessus autour de la solution initiale r et en négligeant les termes du second ordre, la contrainte liée à la courbe temps-profondeur se résume à l'inégalité suivante linéaire en r.: l Nk-1 Nk Zk- 2.t Sri ( V j) = IZk - Zks(R +àR) < ck i=l j=i+l dans laquelle V . est la vitesse de propagation calculée à partir des coefficients r 1 obtenus lors de la première étape d'après la formule: j- 1 v j = V1 1 + r 1=1 r En résumé, la seconde étape consiste donc à minimiser la fonction n n J(6R) = Z ai I6ri avec AR = {rit i=1 -18- sous les contraintes linéaires en àR suivantes: Umin(t) < ( R + R) * Dm Y(t) <Umax(t) et jIZk - Zks(R + àR) < ck Toutefois, la norme J(AR) étant exprimée en fonction de valeurs absolues et ne constituant donc pas une fonction dérivable de ri, il faut décomposer l'incrément Sri en deux valeurs positives 6ri+ i i i et &ri., et ce quelle que soit la valeur de l'indice i, de sorte que: Sri = Sri 6ri avec Sri+ >0 et 6ri >0 pour i= 1, n Il est aisé de-démontrer que le minimum de J(AR) est atteint lorsqu'au moins une des valeurs soit Sri, soit +ri, est égale à zéro pour chaque valeur de l'indice i. Par conséquent l'on peut écrire:

1r1 = r i I 6 ri- i ri ri-

Le problème se ramène finalement à trouver deux séries R+= 6ri+ et ÈR- = /ri i qui minimisent la fonction n J(AR) = Z ai (6ri+ ri-) - i=l tout en vérifiant les équations de contraintes Umin(t) < ( R + R+ - R). Dm Y(t) < Umax(t) et IZk - Zks(R + R+ - R-)| < ck -19- Le processus de minimisation de la fonction J(AR) peut être effectué à l'aide d'un programme d'ordinateur disponible dans le commerce qui codifie l'un des nombreux algorithmes de programmation linéaire approprié. A titre illustratif, le demandeur fait usage du logiciel IMSL (International Mathematical & Scientific Libraries) qui

codifie l'algorithme bien connu qu'est celui du simplexe.

La solution finale R est obtenue en additionnant (9) à la solution initiale R = t it obtenue lors de la première étape (6), l'incrément àR =.6rii obtenu lors de la seconde étape (8), soit

R = ri: = /ri + &ri/ = {r i + ri+ - Sri- La Fig. 6A représente en échelle temps, la série de coefficients r.

obtenue à partir des traces Sk(t) illustrées dans la

Fig. 5.

Minimisation de la fonction H(R) (Fig. 2)

La description détaillée de la minimisation de la fonction

H(R) va maintenant être faite en référence à la figure 2 et à l'aide

d'un exemple pratique pour illustrer sa mise en oeuvre.

L'équation matricielle des contraintes F. R = Um Y()/Dm Y(0) peut se représenter dans R256 comme un polyèdre dont chaque sommet est défini par 256 coordonnées, parmi lesquelles m=21 sont non nulles et appelées variables de base, tandis que les 235 coordonnées restantes sont nulles et appelées variables hors base. Un théorème classique de la programmation linéaire, qui en l'occurrence est applicable à une fonction concave telle que H(R), montre que la fonction H(R) admet un minimum absolu qui est situé sur l'un des sommets du polyèdre, c'est à dire que les coordonnées de la série R recherchée correspondent à un sommet du polyèdre des solutions dont m= 21 coordonnées sont non nulles

et n-m=235 coordonnées sont nulles.

-20- Ainsi, la minimisation de la fonction H(R) va consister à considérer successivement des séries R sommets du polyèdre des solutions pour aboutir à la solution RI tel que la fonction H(R )

atteigne son minimum absolu.

En s'appuyant sur la méthode du gradient réduit qui est décrite dans les pages 345-358 du livre intitulée "Linear and Nonlinear Programming" de D. G. Luenberger, 1984 paru chez "Addison-Wesley publishing Company", il est donc possible de segmenter l'ensemble R en deux sous-ensembles RO et R1 comprenant respectivement 235 coefficients r. nuls et 21 coefficients r. non nuls. L'on peut donc reposer le

1 1

problème de minimisation de la fonction H(R) sous la forme suivante: minimiser H(R) = H(RO,R1) = H(RO) + H(R1) sous la contrainte F. R = (F1. R1, FO. RO) = Um Y () /Dm Y () dans laquelle les sous-matrices F0 et F1 correspondent respectivement

aux sous-ensembles ROet R1 pour les fréquences envisagées.

A partir des ondes montante Umoy(w) et descendante Dmoy(w) (4), on établit (10) une segmentation d'indice O en déterminant 21 positions d'interfaces, puis en calculant la série R10 correspondante en résolvant l'équation R10 [F10-l1. Um Y(w)/Dm Y("), dans laquelle Fl0 est la sous- matrice de F0 relative à cette segmentation. On déterminera cette segmentation d'indice 0 de façon à ce que Fl0 soit une matrice numériquement inversible, par exemple en utilisant la méthode préconisée dans l'article précité "Inversion of Vertical

Seismic Profile by Iterative Modelling" de P. Grivelet (page 926).

Puis, selon un procédé itératif, il s'agit de trouver à l'itération p une nouvelle estimation Rlp+1, ROp+1 à partir de la segmentation Flp, FOp (11) : -21- - pour laquelle la valeur de la fonction H(Rp+i) est inférieure à celle de H(Rp); - qui comprenne la même proportion de coefficients non nuls et nuls (21 et 235 respectivement), c'est à dire qui corresponde à un nouveau sommet du polyèdre; - et qui diffère de la solution précédente en ce qu'un coefficient de Rl devient nul, tandis qu'un coefficient de RO devient

P P

non nul.

A cet effet, -l'on calcule (12) la valeur de H(Rlp) et son

gradient VH(Rlp) pour chacune des valeurs de base de Rlp.

Simultanément, l'on calcule (13) la matrice G qui est le produit de la i

matrice inversée [Flp] 1 par la matrice FOp.

p P' En effectuant (14) le produit G. VH(Rl), l'on obtient un vecteur g composé de 235 coordonnées, indicatif de la variation de la

H(Rp) selon ces 235 coordonnées.

- Le vecteur g est alors soumis à un test (15) pour savoir s'il est nul ou non. Si toutes les 235 coordonnées de g sont nulles, l'on est en présence du minimum absolu de la fonction H(R), ce qui correspond à la solution R0 recherchée (16). Si toutes les 235 valeurs - de g ne sont pas nulles, l'on sélectionne (17) la coordonnée maximale gj parmi les 235 coordonnées de g qui représente, très probablement, la direction de plus grande pente de la fonction H(R) et, par conséquent, celle qui, très vraisemblablement, est à suivre pour minimiser le plus

rapidement la fonction H(R).

La direction dj de la pente gj est un vecteur de dimension 21, calculée (18) de la manière suivante: dj = G. gj = {dj(i)., i = 1 à 21} -22- En suivant cette direction, il s'agit alors de déterminer (19) la distance X à parcourir sur cette pente pour atteindre un autre sommet du polyèdre. Autrement dit, il s'agit d'annuler une variable de base pour la rendre hors base, et inversement de transformer une variable hors base en une variable de base non nulle. Pour cela, il existe 21 valeurs possibles de X. xi, i=l à 21} (24) qui permettent d'annuler le coefficient de reflexion correspondant rll(i) non nul, Xi étant donné par la formule: i = rll(i) / dj(i) Pour chacune des 21 valeurs de X. ainsi obtenues, on tire (20) i une nouvelle série RI 1 selon la formule: Pour chacune des 21 séries Rp+l ainsi obtenues, on calcule i Ip+1

(21) la valeur de H(Rlp+l) correspondante.

Parmi les 21 valeurs de H(R1p+l1), on retiendra (22) uniquement celle qui est la plus faible d'entre elles, c'est à dire celle qui permet de réduire la valeur de la fonction H(R) le plus rapidement, et on la rebaptise alors H(Rlp+1), avec p *25 H(Rl +1) = minimum des /H(R pi), i= 1,21/ p îp+i Un test (23) permet alors de vérifier si H(Rlp+1) est plus petit que H(Rlp) obtenue au cycle précédent. Si le test est positif, la nouvelle segmentation Rlp+1 est adoptée pour un nouveau cycle avec p=p+l (26), et l'on vient à nouveau se placer au bloc (11). Par contre si H(Rlp+ 1) est plus grand que H(Rlp), on annule la coordonnée gj de g (25), et l'on vient à nouveau se placer au dessus du test (15) pour

- verifier si oui ou non toutes les coordonnées du vecteur g sont nulles.

-23- Lton s'apercevra aisément que les solutions successives Rp, y compris la solution initiale R , comprennent toujours, dans l'exemple

considéré; 21 coefficients non nuls et 235 coefficients nuls.

Variantes du procédé d'inversion (Fig. 3A) Nous abordons maintenant des variantes du procédé d'inversion qui permettent en outre d'injecter des informations provenant d'une ou plusieurs diagraphies dans la première étape et/ou celles provenant de la sismique de surface dans la seconde étape, Selon une première variante, il est possible d'orienter le choix de la solution initiale Ra=-r 0i en utilisant les informations contenues dans une ou plusieurs diagraphies, telle que la diagraphie de résistivité, de densité, etc.. En effet, la probabilité est grande pour qu'à une large variation d'amplitude observée sur l'une ou plusieurs diagraphies corresponde une interface, c'est-à-dire une variation d'impédance acoustique. Il est donc tout à fait souhaitable que le profil d'impédance acoustique dérivé de l'inversion du PSV soit cohérent avec les informations contenues dans les diagraphies et

relatives à la localisation des interfaces.

Selon une seconde variante, il est souhaitable que la solution finale R = ril soit cohérente avec les informations disponibles à partir de l'étude sisimique de surface, notamment les vitesses d'addition Va obtenues lors du traitement des lignes sismiques de surface, ou encore avec des informations précises concernant une ou plusieurs couches géologiques, notamment la vitesse de propagation à

travers une ou plusieurs couches déterminées.

Une troisième variante consiste à combiner simultanément les

première et seconde variantes évoquées ci-dessus.

-24- Référence est maintenant faite à la Fig. 3 qui représente l'ordinogramme de la troisième variante, étant entendu que les première et seconde variantes s'en déduisent aisément en supprimant respectivement soit les informations provenant de la sismique de surface, soit les informations provenant des diagraphies. Pour des raisons de simplicité et de concision de l'exposé, on utilisera une seule diagraphie Y(d); toutefois, le procédé d'inversion peut également accepter plusieurs diagraphies selon une chaine de traitement strictement identique et parallèle à celle qui va

être exposée à propos de la diagraphie Y(d).

Les données dont on dispose comprennent: - les vitesses d'addition Va correspondant aux temps ta (301) de la sismique de surface; - les traces sismiques Sk(t) (101) du PSV; - une diagraphie Y(d) (201) exprimée en fonction de la

profondeur du puits.

A partir des traces Sk(t) (101) du PSV, l'on détérmine, d'une manière connue et identique à celle décrite en référence à la figure 1, les valeurs suivantes: - la courbe temps-profondeur >Zk, Tk} (102); - les champs d'ondes montantes Uk(t) et descendantes Dk(t) qui ne contiennent que des ondes primaires (103); - les valeurs um Y(t), Umoy(), Dm Y(t), Dm Y(), umax( et

Umin (t) (104).

A l'aide de la courbe temps-profondeur {Zk, Tk/, l'on effectue un changement d'échelle pour transposer, d'une manière connue, la diagraphie Y(d) en une diagraphie Y(t) - {yi(t)I (202) exprimée en fonction d'une échelle temps et échantillonnée selon un pas identique au pas d'échantillonnage de Sk(t), soit un pas de lms pour l'exemple traité. La Fig. 6D illustre la diagraphie Y(t) ainsi obtenue. Puis l'on calcule (203) la différentielle discrète Y'(t)=ty'i(t)t de la diagraphie Y(t) de la manière suivante: -25- Y'(t) = t Y'i(t) = Yi,+(t) - Yi(t)l

et sa transformée Y'(w) dans le domaine des fréquences.

Il existe une infinité de possibilités en ce qui concerne le nombre des interfaces et leur localisation. Comme il a été exposé précédemment, une solution particulière et préférable est celle qui présente un minimum d'entropie. De plus, il est souhaitable que la localisation des interfaces soit cohérente avec les transitions que l'on peut observer sur une diagraphie zonée en couches géologiques

relativement homogénes.

Nous désignerons par Y=-yi la diagraphie zonée et par

Y'=y'i. sa différentielle discrète.

En reprenant l'approche suivie par la demanderesse à propos de la théorie de l'entropie, le zonage d'une diagraphie peut s'exprimer sous la forme mathématique qui consiste à minimiser une fonction du type entropie, croissante en valeur absolue, symétrique en abscisse et concave au sens large:

A A

H(Y') = E f(yi) et tout particulièrement la fonction suivante

n -''.

^ \. -lY il H(Y') = ( i - e) i=l sous les contraintes linéaires

F. Y' = Y'(@)

-26- Il est à noter, que cette nouvelle approche du problème de zonage est applicable, d'une manière générale, à toute opération de

zonage d'une courbe, indépendamment de l'inversion du PSV.

Dans le cas de l'inversion du PSV, cette approche conduit à un résultat intermédiaire supplémentaire que constitue la diagraphie zonée. Dans le cas de l'inversion du PSV et en vertu de la limitation imposée par la bande de fréquence utilisable [fmin, fmaxl, le zonage devra se réduire à un problème de m =21 interfaces possibles pour l'exemple traité. L'on comprendra aisément que la série Y' devra par conséquent être composée de 21 valeurs non nulles correspondant à 21

interfaces et de 235 valeurs nulles correspondant aux zones continues.

Le problème consiste donc à trouver deux séries R =/r i/ et A A Y'=ty'ii qui présentent la caractéristique selon laquelle les mêmes indices correspondent respectivement aux valeurs nulles et aux valeurs non nulles, et qui minimisent simultanément les deux fonctions du type entropie suivantes: H(R) = E f(ri) A H(Y') = ú e(yi) dans lesquelles f et e sont deux fonctions croissantes en valeur absolue, symétriques en abscisse et concaves au sens large sous les contraintes F. R = Um Y( >/DmoY(W) A E À Y' = Y'(o) -27- dans lesquelles F et E sont les matrices associées à la transformée de Fourrier rapide "FFT" de Cooley-Tuckey, respectivement restreintes à un nombre identique de fréquences sélectionnées. En effet, il est essentiel que chacune des deux matrices soient restreintes à un nombre identique de fréquences sélectionnées de manière à ce que les deux équations de contraintes définissent respectivement deux polyèdres de la même famille, c'est à dire dont les sommets sont définis par n coordonnées, parmi lesquels au plus m sont non nulles, tandis que les

n-m restantes sont nulles.

Pour simplifier la description, on envisagera par la suite

uniquement le cas dans lequel les fonctions f et e sont identiques et égales à la fonction f, et pour lequel les matrices F et E couvrent la même bande de fréquence utile, à savoir celle de la matrice F. Dans un tel cas, il s'agit donc de minimiser simultanément (105) les deux fonctions: H(R) = E f(ri) H(Y') = E f(yi) ou, d'une manière préférentielle, les deux fonctions: n H(R)= ( 1 - e rI) n B(Y') = ( 1 - e) i=l sous les contraintes linéaires suivantes F. R = Um Y(w)/DmoY()

F. Y' =Y'(@)

-28- Le processus de minimisation (105) simultanée de H(R) et de A H(Y') sera décrit en détail en référence aux Figs. 4A et 4B à l'aide d'un algorithme de gradient réduit modifié, étant entendu que ce processus pourrait également s'effectuer, d'une manière similaire, à l'aide d'un algorithme de gradient projeté modifié. Le résultat de cette minimisation (105) conduit à une solution A initiale RD et à la différentielle discrète Y' (106). L'intégration de la différentielle discrète Y' ainsi obtenue, suivie d'un recalage à une constante près, permet de reconstruire la diagraphie zonée Y dont une

représentation graphique se trouve dans la Fig. 6D. Cette dernière fait -

clairement apparaître les différentes couches géologiques.

Dans la deuxième étape (107), il s'agit de minimiser, comme décrit en référence à la Fig. 1, la même fonction linéaire J(AR) suivante n J(R) Zai I ril avec àR = { rit i soumise aux mêmes contraintes que celles exposées précédemment, à savoir: Umin(t) < ( RO + R) * Dmo Y(t) < Umax(t) et IZkZks(R + àR)1 < ck En outre, la minimisation de J(àR) peut également être soumise à une ou plusieurs contraintes additionnelles qui sont linéaires et qui proviennent, soit des vitesses d'addition obtenues suite au traitement de la sismique de surface, soit des contraintes liées à la connaissance préalable de la vitesse de propagation à travers une ou plusieurs

couches géologiques déterminées, soit des deux.

-29- En ce qui concerne les vitesses d'addition, on désignera par Va les vitesses d'addition correspondant aux temps ta. Il est connu que les vitesses d'addition sont très proches des racines carrées du moment d'ordre deux de la série des vitesse Vi. Par conséquent, de manière similaire à la contrainte issue de la courbe temps-profondeur, il est possible de poser que la différence entre la vitesse d'addition réelle Va et la vitesse d'addition synthétique Vas calculée d'une manière indépendante à partir des vitesses de propagation Vi, doit être

inférieure à une grandeur cv qui se déduit de l'incertitude sur Wa.

Cette contrainte peut donc prendre la forme de l'inégalité suivante: Na INa. Va2- E Vi2| < cv i=l dans laquelle Na représente le temps Ta exprimé en nombre d'échantillons selon le pas d'échantillonage Et, soit 1 ms dans l'exemple traité, et Vi la vitesse de propagation de l'onde dans

l'intervalle [i, i+l].

En reprenant le raisonnement mathématique utilisé précédemment dans le cadre de la linéarisation de la contrainte liée à la courbe tempsprofondeur >Zk, Tk}, la contrainte liée aux vitesses d'additions se résume à l'inégalité suivante linéaire en Sri: Na Na-1 Na 6t V Vi2+ 4.t L ri ( Vj) = Vas2(R%+AR) i=l i=l j=i+l ou encore Va - Vas(R%+R) < cv -30- dans laquelle V j est la vitesse de propagation calculée à partir des

coefficients r 1 obtenues lors de la première étape d'après la for-

mule: j-1 V ' = V1 n 1 + r 1 J 1 - r 1 1=1 En résumé, la variante de la seconde étape consiste donc à minimiser la fonction: n J(:R) Z ai I6ril avec AR = Srit i=l sous les contraintes linéaires en &R suivantes: Umin(t) < ( R + AR DmoY(t) < umax(t) jZk - Zks(R + AR)j < ck |wa - was(R + aR) | < cv Dans le cadre de cette variante, le processus de minimisation de la fonction J(AR) s'effectue également à l'aide d'un algorithme de programmation linéaire, tel que l'algorithme du simplexe codifé dans.le

logiciel IMSL évoqué précédemment.

La Fig. 6B représente la série de coefficients de réflexion ri obtenue selon la variante qui a été décrite ci-dessus. On observera tout particulièrement la cohérence entre la position des coefficents de réflexion de la Fig. 6B et la position des interfaces de la diagraphie

zonée Y(t) sur la Fig. 6D.

-31- ^

Ninimisation simultanée de H(R) et H(Y') (Figs. 4A, 4B).

L'on décrira maintenant, d'une manière détaillée en référence aux Figs. 4A et 4B, le processus de minimisation simultanée des A fonctions H(R) et H(Y') qui correspond à la variante de la première étape.

Par analogie à la description en référence à la Fig. 3, les

équations de contraintes peuvent se représenter comme deux polyèdres, les sommets de chacun de ces polyèdres étant définis par 256 coordonnées. Ainsi le processus de minimisation simultanée va consister

256 A 256

à parcourir successivement les sommets des deux polyèdres R256 et Y256

pour aboutir aux séries R et Y' de faible entropie.

Le traitement préalable décrit en référence à la Fig. 3 conduit aux ondes moyennes montante Umoy(w) et descendante Dmoy(o) (104), ainsi qu'à la différentielle discrète Y'(X) (203). L'on détermine alors les maxima de la différentielle Y'(X) en vue de diriger la localisation des interfaces et d'obtenir une segmentation initiale

d'indice p=O de la matrice F en deux sous-matrices F00, Flo (112).

A partir de cette segmentation initiale, on calcule simultanément les séries Ro10 =.[Flo]-l.Um Y(w)/Dm Y(w) (115) et Y'10=[F10]-.Y'(w) (116) dont chacune correspond à un sommet de son polyèdre respectif et qui vont servir à démarrer un processus itératif

similaire à celui décrit en référence à la Fig. 2.

Dans le processus itératif, il s'agit de trouver à l'itération p deux nouvelles séries <Rlp, R0pt et 'lpy, Y' 0p correspondant à la segmentation Flp, FOp, à partir de la segmentation précédente Flp_ FOp_1 - qui réduisent simultanément les valeurs des fonctions H(R) et H(Y'); -32- qui comprennent chacune un nombre de coefficients nuls et non nuls identique à celui de la série d'indice p-l; - pour lesquelles il existe une correspondance stricte à propos des indices concernant les coefficients nuls et non nuls; - qui annulent un coefficient de chacune de deux séries portant un même indice, tout en permettant à un autre coefficient nul de chacune des deux séries portant également un même indice de devenir non-nul. / A cet effet, l'on calcule la valeur de la fonction H(R1) et de son gradient VH(Rlp) (115), ainsi que celle de la fonction H(Y'lp) et de son gradient VH(Y'lp) (215). Simultanément on calcule (116) la matrice G qui est le produit de la matrice inversée [Fl] 11 par la p

matrice FO.

p En effectuant les produits G.VH(Rlp) (117) et G.VH(Y'lp), l'on obtient respectivement deux vecteurs gr et gy, chacun des vecteurs étant composé de 235 valeurs, respectivement indicatifs de la

variation, selon ces 235 coordonnées, de la fonction H correspondante.

L'on calcule un vecteur g (118) dont les 235 coordonnées sont égales au produit des coordonnées respectives de g et gy. Puis le r Y vecteur g est soumis au test (119) pour savoir si il est nul ou non. Si toutes les 235 coordonnées de g sont nulles, l'on est en présence des ^ solutions recherchées R , Y' (120). Au contraire, si toutes les coordonnées de g ne sont pas nulles, l'on sélectionne (121) la

coordonnée maximale en valeur absolue gj du vecteur g.

Puis par analogie au processus décrit en référence à la Fig. 2, on évalue successivement: - les directions drj=G.gr,j (122) et dyj=G.g yj (222), chacune des directions étant représentée par un vecteur de 21 coordonnées. -33- - pour chacune des valeurs de l'indice i (130), les distances Xr,i = rlp(i)/dri(i) (123) et lp(i)/dr,j(i) Yi = Y'îp(i)/dyj(i) (223) à parcourir selon les directions respectives drj(i) et-dyj(i) pour atteindre le sommet suivant du polyédre correspondant; - pour chacune des valeurs de l'indice i (130), les séries i = R1 -X. 14,e Rlp+1 i Rlp Xr,i dr,j(i) (124), et Y'p+l Y'lp - xy,.dyj(i) (224) - pour chacune des valeurs de l'indice i (130), la valeur des fonctions H(Rlp i) (125) et p+ 1.

H(Y'p+l) (225).

À p i - pour chacune des valeurs de l'indice i (130), la valeur p+1 (126) produit des fonctions H(R1) et H(Y') soit i. p+1 . p+l ' p+l1 = H(Rlp+il). H(Y'p+1) Parmi les 21 valeurs de p+l, on retiendra uniquement celle qui a la plus faible valeur d'entre elles (126) que l'on rebaptise alors Pp+1, soit Pp+1 = minimum des {Pp+l pour i=1 à 21 Un test (128) permet alors de vérifier si l'inégalité P+l <

P est vérifiée ou non.

p Si le test (128) est positif, la nouvelle segmentation Rlp1 est adoptée pour un nouveau cycle avec p=p+l (129) et l'on revient se placer au bloc (114). Si le test (128) est négatif, l'on annule la coordonnnée gj (130), et l'on vient se replacer au-dessus du test (119) pour vérifier si oui ou non toutes les coordonnées du vecteur g sont nulles.

26 1 6920

-34- Une alternative à ce dernier test (128) consiste, non pas à vérifier que la valeur du produit des deux fonctions H(R) et H(Y') diminue par rapport à sa valeur à l'itération précédente, mais à vérifier, d'une manière plus contraignante, que la valeur de chacune des deux fonctions, pour un même indice i, diminue simultanément par

rapport à sa valeur à l'itération précédente.

- Comme évoqué précédemment, la description des variantes-a

porté uniquement sur le cas dans lequel les fonctions f et e sont identiques et égales à la fonction'f, et pour lequel les matrices F et E couvrent la même bande de fréquence utile, à savoir celle de la matrice F. Néanmoins, il est tout à fait possible de minimiser les entropies définies par deux fonctions f et e distinctes qui soient toutes les deux croissantes, symétriques en abscisse et concaves; il est également possible de restreindre les matrices F et E à des fréquences distinctes, aussi longtemps que le nombre de fréquences

sélectionnées est identique pour les deux matrices.

-35-

Claims (10)

REVENDICATIONS

1. Procédé pour identifier les couches géologiques souterraines à partir d'une pluralité-de traces sismiques Sk(t) enregistrées dans un sondage à différents niveaux k de profondeur Zk, chaque trace Sk(t) étant composée de- n valeurs discrètes échantillonées en fonction du temps t selon un pas constant At, le procédé consistant à déterminer la série R = {ri, i = 1,n' représentative des coefficients de réflexion sismique des interfaces formées par les couches géologiques successives, i étant un indice relatif à la profondeur échantillonnée en échelle temps selon ledit pas constant At, ledit procédé comprenant les étapes suivantes: - établir la courbe temps-profondeur reliant, pour chaque niveau k, la profondeur Zk au temps de transit Tk relevé sur la trace Sk(t); - extraire des traces Sk(t), le champ d'ondes montantes Uk(t) et le champ d'ondes descendantes Dk(t);

- calculer les traces umoY(t), umin(t) et umax(t) représen-

tatives respectivement d'une valeur moyenne, minimale et maximale du champ d'ondes montantes Uk(t-Tk), et la trace Dm Y(t) représentatif d'une valeur moyenne du champ d'ondes descendants Dk(t+Tk); caractérisé par les étapes suivantes qui consistent à: a) dans une première étape, déterminer une série R ={r i représentative d'une valeur approchée de la série R en minimisant une fonction B(R)= ú f(ri) dans laquelle f(ri) est une fonction croissante de Iril, symétrique en abscisse et concave, sous la contrainte R * DmoY(t) = UmoY(t) dans laquelle l'opérateur * représente un produit de convolution; b) dans une seconde étape, déterminer une série AR=/Sri/ représentative d'une valeur corrective de la série R en minimisant une fonction J(àR) linéaire en l6ril sous les contraintes: -- 36- min Umax()e U (t) < I(R + R) * Dm Y(t)] < umax(t) et IZk - Zks(R +.) | < ck dans laquelle Zks(R +àR) est une fonction linéarisée en àR, représentative de la profondeur Zks calculée d'une manière synthétique au niveau k à partir de la solution initiale R , et la grandeur ek étant représentative de l'incertitude sur les grandeurs de Zk et Tk; c) calculer la série R = R +AR <ri = r i+&ri/ représentative

des coefficients de réflexion sismique.

2. Procédé selon la revendication 1, caractérisé par le fait qu'il comprend en outre l'étape consistant à transposer les valeurs Um Y(t) et Dm Y(t) dans le domaine des fréquences de manière à obtenir les valeurs UmoY(w) et Dm Y(w), et que les contraintes de la première étape sont exprimées dans le domaine des fréquences, à savoir: F.R = UmoY ()/ Dm Y () dans laquelle F est la matrice associée à la transformée de Fourier,

restreinte à la bande de fréquence utile des traces sismiques.

3. Procédé selon l'une quelconque des revendications 1 ou 2,

caractérisé par le fait que la fonction J(MR) est de la forme - J(AR) = E ai6ril dans laquelle ai est un coefficient de pondération dépendant de ri. ,

4. Procédé selon l'une quelconque des revendications 2 ou 3,

caractérisé par le fait qu'il comprend en outre les étapes suivantes: - à l'aide de ladite courbe temps-profondeur >Zk, Tk}, transposer au moins une courbe de diagraphie Y(d) enregistrée préalablement en fonction de la profondeur d dans ledit sondage sur un intervalle de profondeur qui recouvre notamment les niveaux de profondeur Zk, en une diagraphie discrète Y(t) = <yi(t)" échantillonnée en échelle temps selon le même pas At; -37- - calculer la différentielle discrète Y'(t)= /y'i(t)/ de Y(t), ainsi que sa transformée dans le domaine des fréquences Y'(X); et par le fait que la première étape consiste à déterminer deux séries R {r ii et Y'-y'i} respectivement représentative d'une solution initiale approchée et de la différentielle discrète de la diagraphie zonée Y(t) et pour lesquelles les valeurs nulles de r i et y'i portent le même indice i, en minimisant simultanément les fonctions H(R) =.Zf(r) (Y') Ze(y'i) dans laquelle e(y'i) est une fonction croissante de IY'i[, symétrique en abscisse et concave, sous les contraintes F.R = Um Y()/DmoY(w) A

E.Y' Y'(X)

dans laquelle E est la matrice associée à la transformée de Fourier, restreinte à un nombre identique de fréquences à celui de ladite bande utile.

5. Procédé selon l'une des revendications 3 ou 4, caractérisé par

le fait que la seconde étape consiste à minimiser la fonction J(MR) sous une contrainte additionnelle IVa - Vas(R +AR)I < cv dans laquelle Va sont les vitesses d'addition correspondant aux temps a, obtenues suite au traitement de la sismique de surface, et Vas(R +6R) est une fonction linéarisée en AR,' représentative de la vitesse d'addition Vas calculée d'une manière synthétique à partir de la solution initiale R , et la grandeur cv étant représentative de

l'incertitude-sur la vitesse d'addition Va.

6. Procédé selon l'une des revendications précédentes,

caractérisé par le fait que la fonction f est la fonction: -Irfi f(ri) = (1 - e) -38-

7. Procédé selon l'une quelconque des revendications 3 à 6, -

caractérisé par le fait que les matrices E et F sont restreintes à la

même bande de fréquence et, par conséquent, identiques.

8. Procédé selon l'une quelconque des revendications 3 à 7,

caractérisé par le fait que le coefficient de pondération ai est égal à la plus grande des. deux valeurs 1/r min et 1/r i, la valeur r min

étant la valeur minimale des coefficients r i non nuls.

9. Procédé pour identifier les couches géologiques souterraines en effectuant une opération de zonage d'une diagraphie discrète Y(t)=/yi(t), i = 1, ni échantillonée selon un pas constant At, de i AA manière à en extraire la, série Y'=ly'i, i = 1, ni indicative de la position des interfaces formées par les couches géologiques successives, caractérisé par les étapes suivantes: - calculer ia différentielle discrète Y'= ly'Ii de Y, ainsi que sa transformée dans le domaine des fréquences Y'(X); déterminer la série Y' indicative de la position des interfaces en minimisant une fonction A H(Y')= Z e(y'i) dans laquelle e(y'i) est une fonction croissante de i[, symétrique en abscisse et concave, sous la contrainte

E.Y' = Y'()

dans laquelle E est la matrice associée à la transformée de Fourier,

restreinte à un nombre de fréquences sélectionnées.

10. Procédé selon la revendication 9, caractérisé par le fait que la fonction e est la fonction: ^-y' i e(y'.3 = (1 - e