JMM @ EOST

Sismique : Ondes - Imagerie - Traitement - Propagation - Applications

Sismique réflexion : Introduction - Dispositifs d'acquisition, coordonnées et collections d'enregistrements
Migration à V constante et déport nul : Points diffractants - Réflecteurs plans - TF - Echantillonnage - Stolt - Kirchhoff
Migration à V constante et déport non nul : Point de tir - DMO
Migration à V(z) et déport nul : Relation temps-profondeur - Rayons et "phase-shift"
Migration à V(x,y,z) et déport nul : Faibles variations latérales - Equation iconale - Equations paraxiales - Différences finies 15° - Renversement temporel
Migration à déport non-nul : Equations de prolongement - Imagerie - Réfraction
Transformation : τ-p - Séparation montante-descendante
Documents : Figures et programmes - Sujets exams - Livres et liens - Documentation interne EOST

Les figures grand format s'ouvrent toutes dans la même fenêtre en cliquant sur les icônes. Les animations Java s'ouvrent en cliquant sur les liens en gras. Les programmes Matlab écrits pour réaliser les figures sont accessibles en cliquant sur le nom du programme (prog.m). Pour forcer l'affichage de la version la plus récente de la page, appuyez sur Shift en cliquant sur Actualiser. Pour imprimer correctement les formules et tableaux, utiliser le format paysage (feuille en travers).

Introduction

L'imagerie sismique du sous-sol par sismique réflexion à partir de la surface utilise les ondes réfléchies pour déterminer la géométrie des interfaces géologiques avec une bonne résolution latérale. On réalise une série d'expériences où les échos provoqués par une source ponctuelle sont enregistrés sur de nombreux récepteurs échantillonnant une distance (ou surface pour la 3D) limitée. Le traitement des données pour obtenir une image du sous-sol comporte plusieurs étapes . Un préalable indispensable est de pouvoir déceler et séparer les échos de faible amplitude en provenance de la profondeur parmi l'ensemble des ondes enregistrées. L'imagerie n'utilise que les réflexions primaires ayant effectué un trajet direct source-réflecteur-récepteur. Si l'on réussit à les isoler, la question suivante est : d'où proviennent les échos? Transformer les temps de parcours des ondes réfléchies en une image des réflecteurs en profondeur est le but premier de l'étape du traitement sismique appelée migration. On donne ici les bases essentiellement géométriques des différentes méthodes de migration pour un milieu de propagation acoustique. On ignore dans un premier temps les aspects concernant la préservation des amplitudes. On suppose aussi que la vitesse de propagation V est connue.

Les différences entre la structure du sous-sol en profondeur et l'image sismique obtenue en juxtaposant les enregistrements en temps faits par des couples source-récepteur confondus sont illustrées sur les figures suivantes (les temps correspondent aux temps aller-retour dits temps doubles).

reflecteur plan - reflecteur courbe - reflecteur concave - isis.m

Pour un segment de réflecteur plan horizontal, une conversion de l'axe des temps en profondeur suffit à ramener la réflexion à la position du réflecteur. Cependant, les extrémités du réflecteur se comportent comme des point diffractants qui renvoient l'onde incidente dans toutes les directions. Les hyperboles de diffraction doivent être focalisées sur les points diffractants par la migration.
Pour un segment de réflecteur plan penté, il y a un décalage latéral du au trajet oblique des rayons réfléchis alors que la réflexion est imagée en temps à la verticale des couples source-récepteur. L'extension latérale de la réflexion est plus grande que celle du réflecteur. De plus, la pente apparente en temps est différente du pendage du réflecteur.
Pour un réflecteur courbé (ici un arc de cercle), une forme concave focalise les rayons réfléchis. La réflexion résultante est donc d'extension latérale plus restreinte que le réflecteur. Au contraire, une forme convexe défocalise les rayons réfléchis. La réflexion résultante est donc d'extension latérale plus grande que le réflecteur.
Lorsque la courbure d'un réflecteur concave est telle que le foyer coïncide avec la surface, l'image en temps est focalisée en un point. Si le foyer est situé sous la surface, la courbure de la réflexion est inverse de celle du réflecteur. De plus, l'image est inversée, les points de réflexion situés sur la droite du réflecteur étant imagés sur la gauche de la réflexion.
Sur toutes ces figures, on observe que l'hyperbole de diffraction est tangente à la réflexion pour le couple source-récepteur pour lequel le rayon réfléchi coïncide avec un rayon diffracté issu du point de réflexion. Cette propriété est à la base des méthodes de migration qui reconstruisent la structure du sous-sol comme un ensemble de points diffractants.

Des expériences utilisant la propagation d'ondes à la surface de l'eau permettent de se familiariser avec les notions utilisées. Versez de l'eau dans un plat en verre de forme rectangulaire et munissez-vous d'une baguette et de deux peignes à dents très espacées.

Un point diffractant (point rouge) émet une onde diffractée sphérique lorsqu'il est atteint par le front d'onde sphérique émis par une source située en surface (point bleu). Sur les figures suivantes, les fronts d'onde et les enregistrements faits en surface sont représentés au temps où :
- l'onde directe atteint le point diffractant
- le front d'onde diffracté atteint la surface (point vert)
- le front d'onde diffracté atteint le récepteur à x = 900m (point vert).

diffraction - diffraction -

Inversement la rétropropagation du front d'onde diffracté fait converger le front d'onde sur la position du point diffractant au temps où il a été émis. On obtient ainsi l'image du point à partir de l'onde qu'il a émise. L'imagerie sismique est basée sur ce phénomène, appelé migration.

migration - migration -

Les interférences entre les ondes diffractées issues de points diffractants proches sont constructives sur la réflexion provenant du réflecteur formé par la juxtaposition des points diffractants.

diffraction-réflexion - sdr.m

Dispositifs d'acquisition, coordonnées et collections d'enregistrements

Le dispositif d'acquisition des données est constitué d'un ensemble de récepteurs (géophones enregistrant la vitesse de déplacement du sol à terre, hydrophones enregistrant la variation de pression en mer) disposés en ligne le long d'un profil (cas 2D) ou en lignes parallèles le long d'un andain (cas 3D). Un enregistrement élémentaire, fait par un groupe de récepteurs connectés pour fournir un seul signal électrique, est appelé une trace. Chaque tir est enregistré par plusieurs centaines de traces. La collection d'enregistrements résultant d'un même tir est dite collection de traces en tir commun.
En mer, un bateau tire un cable d'enregistrement de plusieurs kilomètres de long (plusieurs cables parallèles en sismique 3D) permettant l'enregistrement simultanés de centaines de traces séparées par un intervalle de l'ordre de la dizaine de mètres. La source est une impulsion de pression générée par un réseau de canons à air (voir par exemple cette
vidéo).
A terre, on procède de la même façon, mais le dispositif récepteur peut être disposé de façon symétrique par rapport à la source, souvent constituée de camions vibrateurs (voir par exemple cette vidéo).

Les données sont acquises selon le principe de la couverture multiple, afin que chaque point de réflexion soit illuminé selon un éventail d'angles d'incidence. La source est actionnée de façon à ce que l'intervalle entre tirs Δs soit un multiple de l'intervalle entre traces Δg. Ainsi, les points milieux géométriques des couples sources-récepteurs échantillonnent le profil avec un intervalle Δy = min(Δs , Δg)/2 régulier. Le tri des traces en point-milieu commun permet de constituer des collections de traces ayant le même point-milieu et des déports sources-récepteurs variant entre les distances source-récepteur minimum et maximum. Dans un milieu homogène ou stratifié à interfaces horizontales, les réflexions sur une collection de traces en point-milieu commun ont des points de réflexion situés à la verticale du point-milieu. Le nombre de traces ayant un point-milieu commun est l'ordre de couverture.

Le recouvrement des positions de tir et de récepteurs le long d'un profil rend nécessaire une représentation en plan de la géométrie du dispositif d'acquisition où les axes correspondent respectivement aux abscisses des récepteurs et des tirs le long du profil (coordonnées terrain) ou à l'abscisse des points-milieux et au déport source-récepteur (coordonnées utilisées pour le traitement). Voici des exemples de géométries pour des dispositifs d'acquisition mobiles symétriques, dissymétriques et fixes dans les coordonnées (récepteurs, tirs) et (points-milieux , tirs), le déport source-récepteur étant figuré en couleur.

dispositif d'acquisition symétrique - dispositif d'acquisition dissymétrique - dispositif d'acquisition fixe - sgyh.m

Sismique 2D : sources et récepteurs placés le long d'un profil, abscisse ou distance mesurée le long du profil sources
s = xs
récepteurs
g = xr
point milieu
y = (s+g)/2
demi-déport
h = (g-s)/2
collections ou sections avant sommation source commune
p(s0,g,t)
récepteur commun
p(s,g0,t)
point milieu commun
p(y0,h,t)
déport constant
p(y,h0,t)
Exemples de données acquises lors de la campagne sismique EW0108 Gulf of Corinth Profil GOC11, effectuée en 2001 sous la responsabilité scientifique de Brian Taylor (University of Hawai'i).
Les données non traitées (points de tir) sont archivées au Marine Geoscience Data System, Lamont-Doherty Earth Observatory of Columbia University. Les profils traités (stacks et migrations) sont archivés par Shipley, T., Gahagan, L., Johnson, K., and Davis, M., editors, 2005, Seismic Data Center, University of Texas Institute for Geophysics
EW0108-GOC11 pt8000   EW0108-GOC11 tr240

Les résultats scientifiques de la campagne sismique EW0108 Gulf of Corinth, destinée à l'imagerie des sédiments, du socle et des failles dans une zone d'extension active, sont décrits dans les publications suivantes :
Taylor, B.; Goodliffe, A.; Sachpazi, M.; Hirn, A.; EW0108 Science Party, 2001, First Results from EW0108 Marine Seismic Survey of the Gulf of Corinth, Eos Trans. AGU, 82(47), Fall Meet. Suppl.,
Zelt, B. C., Brian Taylor, B., Weiss, J. R., Goodliffe, A. M., Sachpazi, M., and Hirn, A., 2004, Streamer tomography velocity models for the Gulf of Corinth and Gulf of Itea, Greece, Geophysical Journal International, 159, p. 333–346
Taylor, B., Weiss, J.R., Goodliffe, A.M., Sachpazi, M., Laigle M. and Hirn A., 2011, The structures, stratigraphy and evolution of the Gulf of Corinth rift, Greece, Geophysical Journal International, 185, p. 1189–1219

Points diffractants explosants

Une inhomogénéité localisée en M dans un milieu de vitesse constante V agit comme une source ponctuelle secondaire. Elle émet une onde sphérique au temps où elle est atteinte par l'onde incidente. Lorque le déport est nul, les trajets aller-retour SMR peuvent être remplacés par des trajets simples MR dans un milieu de vitesse moitié. Tout se passe comme si une source placée en M explosait à t=0. La migration consiste à concenter l'onde diffractée au point M au temps t=0 où elle a été émise. En répétant cette procédure pour chaque point du sous-sol, on obtient une image du plan (x,z).

cône de diffraction - cone.m hyperbole de diffraction en z=0 position du point diffractant à t=0
V2t2 = (x-xm)2+(z-zm)2 td = (1/V)(zm2+(x-xm)2)½ x = xm , z = zm
migration de p(x,t) par sommation des diffractions émises à t=0 par chacun des points du sous-sol sommation sur la courbe de diffraction image de M = onde(xm,zm,t=0)
algorithme de sommation sur trajectoire de diffraction dans le plan xt : ∀(xm,zm) S = ∑xp(x,td(x,xm,zm)) p(xm,zm) = S

La sommation est constructive si :

De façon équivalente à la sommation sur les trajectoire de diffraction, chaque échantillon du plan (x,t,z=0) peut être placé sur le réflecteur semi-circulaire dont il peut provenir dans le plan (x,z,t=0). La superposition des demi-cercles produit l'image migrée.

cone de migration - migration par superposition 1/2cercles - migsdc.m réflecteur semi-circulaire à t=0 correspondant à une impulsion en x = x0 , t = t0 à z = 0
V2t2 = (x-xm)2+(z-zm)2 zm = (V2t02-(x0-xm)2)½
algorithme de superposition de 1/2 cercles dans le plan xz : ∀(x0,t0) ∀xm : p(xm,zm) = p(xm,zm)+p(x0,t0)

Voici des exemples de migration de réflexions obtenues par la sommation de diffractions pour des points diffractants alignés sur des segments de droite et sur un synclinal.
migration par superposition 1/2cercles - migration par superposition 1/2cercles - migration par superposition 1/2cercles

Voici la correspondance entre un réflecteur en forme de synclinal et la réflexion correspondante pour un milieu de vitesse constante V = 1000 m/s et des temps simples de propagation. Les segments courbés du synclinal sont des arcs de cercle.
synclinal

Le segment de courbe de diffraction qui contribue constructivement à la migration d'une réflexion correspond à la zone de Fresnel dans laquelle l'interférence entre une réflexion plane et une diffraction harmonique est constructive dans le plan (x,t,z=0). Sa largeur s'obtient en écrivant que le déphasage correspondant à la différence de longueur δl entre le trajet réfléchi et le trajet diffracté est inférieur à une demi-période (kδl < π). Soit xF l'écart entre les points d'émergence en z=0 des rayons réfléchi et diffracté issus d'un même point diffractant à la profondeur zm et α l'angle sous lequel xF est vu depuis la profondeur zm :

zone de Fresnel - fresnel.m tgα = xF/zm ≈ α δl ≈ xFsinα ≈ xF2/zm xF2 = πzm/k = λzm/2
zone de Fresnel tgα = xFcos2θ/zm ≈ α δl ≈ xFcosθsinα ≈ xF2cos3θ/zm xF2 = πzm/(kcos3θ)

Babcock, J. M. ; Harding, A. J. ; Kent, G. M. ; Orcutt, J. A. , 1998 An examination of along-axis variation of magma chamber width and crustal structure on the East Pacific Rise between 13°30'N and 12°20'N montrent l'effet des diffractions sur les rugosités de la croûte océanique sur les images sismiques à l'axe des dorsales.
Zhu, J., L. Lines and S. Gray, 1998, Smiles and frowns in migration/velocity analysis illustrent l'effet d'une mauvaise vitesse de migration. Voir aussi ce sujet d'examen et ces figures.

point diffractant migré avec VM ≠ V0 - segment migré avec VM ≠ V0 - migration par superposition 1/2cercles, VM ≠ V0 - migration par superposition 1/2cercles, VM ≠ V0 - migration par superposition 1/2cercles, VM ≠ V0

Le traitement des données du profil EW0108-GOC11 fournit une coupe sismique par transformation des temps de trajet des réflexions primaires à déport nul et sommation (stack) en point milieu commun le long du profil. On y observe simultanément les réflexions visibles sur les points de tir et sur les coupes à déport constant dans l'intervalle entre la réflexion primaire au fond de l'eau et les réflexions multiples dans l'eau : les réflexions dans la partie supérieure des sédiments syn-rift sont vues à faible déport, celles dans la partie profonde (dont l'interface sédiments syn-rift-socle) à grand déport. Quelques pentes correspondant à différents pendages de réflecteurs et diffractions associées aux failles sont indiquées sur la figure suivante. La conversion pendages-dt/dx est faite en utilisant la vitesse constante dans l'eau V = 1500 m/s. Comme la couche d'eau est épaisse au centre du graben, cette approximation convient pour les sédiments peu profonds. La migration à la vitesse de l'eau fournit une image correcte à faible profondeur, où l'on distingue clairement les blocs découpés par les failles normales. Par contre la partie profonde est sous-migrée.

EW0108-GOC11 stack - EW0108-GOC11 migration V0

Réflecteurs plans vibrants

La décomposition en points diffractants traite tous les pendages pour chacun des points du sous-sol. Une décomposition en ondes planes traite globalement chaque pendage. Lorque le déport est nul, les trajets aller-retour des ondes réfléchies peuvent être remplacés par des trajets simples issus de sources placées sur les réflecteurs et explosant à t=0. L'image en profondeur est donnée par le champ d'onde prolongé en profondeur au temps t=0. L'utilisation d'ondes planes harmoniques permet de prolonger le champ d'onde en profondeur par un simple déphasage.

. onde(x,z,t) enregistrement en z = 0 image à t = 0
relation de dispersion kx2+kz2 = ω2/V2 kx = ωsinθ/V kz = -ω(1/V2-kx22)½
(- car onde montante)
onde harmonique = 1 pendage p(x,z,t) = Pei(kxx+kzz-ωt) p(x,t) = Pei(kxx-ωt) p(x,z) = Pei(kxx+kzz)
superposition d'ondes harmoniques = tous pendages p(x,z,t) = ∫∫P(kx,ω)ei(kxx+kzz-ωt)dkx p(x,t) = ∫∫P(kx,ω)ei(kxx-ωt)dkx p(x,z) = ∫∫P(kx,ω)ei(kxx+kzz)dkx

Transformée de Fourier

Les facteurs constants, absorbés dans les gains appliqués pour la représentation des traces, sont omis.

TF inverse et directe p réel p réel symétrie translation dilatation convolution en x convolution en k dérivation
p(x) = ∫P(kx)eikxxdkx p(x) p(-x) P(x) p(x-x0) p(ax) p∗q(x) = ∫p(ξ)q(x-ξ)dξ p(x)q(x) p'(x)
P(kx) = ∫p(x)e-ikxxdx = |P(kx)|eiχ(kx) P*(kx) = P(-kx) P*(kx) p(-kx) e-ikxx0P(kx) (1/|a|)P(kx/a) P(kx)Q(kx) P∗Q(kx) ikxP(kx)

La TF détecte les périodicités : la TF de fonctions harmoniques est constituée d'impulsions aux nombres d'ondes correspondants. Une somme infinie d'harmoniques successives produit par interférence un peigne d'impulsions à la longueur d'onde fondamentale.

p(x) eik0x 1 cosk0x sink0x ∑eink0x eink0x = ↑λ0
P(kx) δ(kx-k0) δ(kx) δ(kx-k0)+δ(kx+k0) -i(δ(kx-k0)-δ(kx+k0)) ∑δ(kx-nk0) δ(kx-nk0) = ↑k0

Voici quelques transformées intervenant en sismique

mesure sur [-a,a] autocorrélation signe(kx) échelon
p(x)Π[-a,a] ∫p(x')p(x+x')dx' = p(x)∗p(-x) i/x H(x) = (1+signe(x))/2
P(kx)∗(2/kx)sin(kxa) P(kx)P*(kx) = |P(kx)|2 signe(kx) δ(kx)+1/(ikx)

La décomposition d'une onde sphérique en ondes planes harmoniques s'obtient explicitement par le calcul des composantes P(kx,ky) dans le plan z=0 puis par prolongement de chaque composante en z>0.

r2 = x2+y2+z2 k2 = ω2/V2 = kx2+ky2+kz2 eikr/r = ∫∫P(kx,ky)eikzzei(kxx+kyy)dkxdky = ∫∫(i/kz)eikzzei(kxx+kyy)dkxdky
r2 = x2+y2 , x = rcosθ κ2 = kx2+ky2 , kx = κcosη P(kx,ky) = ∫∫eikr/re-i(kxx+kyy)dxdy = ∫∫eikre-iκrcos(θ-η)drdθ
0≤r<∞ , -π≤θ≤π kz réel donc κ/k<1 P(kx,ky) = ∫dθ∫ei(k-κcos(θ-η))rdr = (i/k)∫dθ/(1-(κ/k)cos(θ-η)) = 2iπ/(k22)½ = 2iπ/kz
a = κ/k , t = tg(α/2) ∫dα/(1-acosα) = 2∫d(tgα/2)/((1-a)+(1+a)tg2α/2) = 2/(1-a)∫dt/(1+((1+a)/(1-a))t2) = 2/(1-a2)½Arctg(((1+a)/(1-a))½t) = 2π/(1-a2)½

Echantillonnage : J F o u r i e r - JFourierJFourier - J u e F r r o i

La TF échantillonnée utilise un nombre fini d'harmoniques pour représenter le spectre. La longueur d'onde fondamentale correspond à la longueur d'enregistrement X, celle de l'harmonique la plus élevée à 2 fois le pas d'échantillonnage Δx (nombre d'onde de Nyquist kxN = π/Δx). Si des longueurs d'onde plus petites sont présentes, la périodisation associée à l'échantillonnage superpose kx = kxN+2nπ/X avec kx = -kxN+2nπ/X (repliement spectral ou aliasing).

onde sur intervalle X périodisation en x → échantillonnage en k échantillonnage en x → périodisation en k longueurs et nombres onde
p(x) p∗↑X↑ p∗↑X↑.↑Δx↑ λx = ∞,X,X/2,..,2Δx
P(kx) P.↑2π/X↑ P.↑2π/X↑∗↑2π/Δx↑ ±kx = 0,2π/X,4π/X,..,π/Δx

Voici l'illustration de la TF échantillonnée pour des harmoniques complexes et réelles et des diracs réels :

TF harmoniques complexes - TF harmoniques reelles - TF diracs reels
proptf.m

et de l'effet de l'aliasing spatial dans les plans (x,t) et (kx,ω) pour une réflexion harmonique :

aliasing spatial
alias.m

Biondi B., 2001, Kirchhoff imaging beyond aliasing illustre l'effet de l'aliasing sur les images et les opérateurs de migration.

Le module de la TF2D d'un point de tir du profil EW0108-GOC11 montre que les arrivées à la vitesse de l'eau (1520 m/s) sont aliasées aux fréquences supérieures à 30 Hz pour un pas d'échantillonnage Δg = 25 m. Il est possible de déaliaser partiellement l'enregistrement en reconstituant le plan f-k correspondant à un échantillonnage spatial deux fois plus dense Δg = 12.5 m, compte-tenu de la dissymétrie de l'enregistrement.

EW0108-GOC11 PT8000 TF2D DG = 25 m - EW0108-GOC11 PT8000 TF2D DG = 12.5 m - EW0108-GOC11 PT8000 DG = 25 ET 12.5 m

Migration de Stolt

En effectuant un changement de variable dans le domaine de Fourier, on remplace le calcul d'une intégrale double par celui d'une TF2D inverse. Le changement de variable implique une interpolation des valeurs complexes de P (source d'artefacts) ainsi qu'une pondération des amplitudes en cosθ (optionnelle). A cause de la périodisation, la section migrée inclut les symétriques des portions de réflecteurs migrés au-dela des limites latérales du domaine de calcul. De plus, le prolongement vers le bas de l'onde montante depuis la surface z=0 jusqu'à la profondeur z est équivalent au prolongement vers le haut de l'onde descendante depuis le bas de la section z=Z jusqu'à la profondeur z. La réponse impulsionnelle comporte donc deux demi-cercles tangents centrés respectivement en z=0 (réflecteur) et z=Z (artefact).

TF2D chgt. variable → interpolation (+ pondération) image à t = 0 par TF2D inverse
P(kx,ω) = ∫∫p(x,t)e-i(kxx-ωt)dxdt ω = -Vkz(1+kx2/kz2)½
P(kx,kz) = P(kx,ω(kx,kz))
dω/dkz = A(kx,kz) = Vcosθ
p(x,z) = ∫∫(AP)(kx,kz)ei(kxx+kzz)dkxdkz

e-ikz(Z-z) = ei(-2π+kzz) = eikzz

Voici l'illustration des étapes de la migration de Stolt p(x,t) → P(kx,ω) → P(kx,kz) → p(x,z) pour :

onde plane monochromatique - somme d'harmoniques - aliasing spatial - stolt.m
une onde plane monochromatique , une somme d'harmoniques pour un pendage fixé (kx/ω = cte) sans et avec aliasing spatial

impulsion impulsion impulsion
une impulsion en (x=0 , t=0.25s) , en (x=0 , t=1s) , au centre de la section.

Voici ce qui se passe si on migre une même section avec différentes vitesses de migration (voir ce sujet d'examen)

impulsion migré par Stolt avec vitesse différente de V - réflexion  migrée par Stolt avec vitesse différente de V

Stolt, R.H., 1978, Migration by Fourier transform.

Migration de Kirchhoff

La migration dans le domaine de Fourier manipule les ondes planes conformément à l'équation des ondes. Une intégrale de Fourier représente une sommation de sinusoïdes qui n'est constructive qu'au voisinage des points stationnaires de l'argument de l'exponentielle complexe. La méthode de la phase stationnaire est basée sur cette propriété (pour une onde, la phase est stationnaire dans la direction de propagation). On passe ainsi de l'expression de la migration dans le domaine de Fourier à celle de la migration par sommation selon les hyperboles de diffraction. Des termes correctifs affectant la phase et l'amplitude apparaissent lorsqu'on fait un développement à l'ordre 2 : ils sont pris en compte dans la migration de Kirchhoff.

calcul de I = ∫eiφ(k)dk par développement de φ(k) autour de k0 tel que φ'(k0) = 0 ordre 0 : φ(k) ≈ φ(k0)
I ≈ eiφ(k0)
ordre 2 : φ(k) ≈ φ(k0)+((k-k0)2/2)φ''(k0)
I ≈ eiφ(k0)∫ei((k-k0)2/2)φ''(k0)dk = eiφ(k0)(2π/|φ''(k0)|)½esignφ''(k0)iπ/4
p(xm,zm) = ∫∫P(kx,ω)ei(kxxm+kz(kx,ω)zm)dkxdω = ∫dx∫P(x,ω)dω∫eiφ(kx)dkx
φ(kx) = -kx(x-xm)+kz(kx)zm

rd(x) = ((x-xm)2+zm2)½ , td(x) = rd(x)/V , tgθd = (x-xm)/zm
kx0 = ωsinθd/V , kz(kx0) = -ωcosθd/V , φ(kx0) = -ωtd(x)

eiφ(kx0) = e-iωtd(x)

p(xm,zm) = ∫dx∫e-iωtdP(x,ω)dω = ∫p(x,td(x))dx = sommation de p sur hyperbole de diffraction

φ''(kx) = kz''(kx)zm = zmV/|ω|(1-kx2V22)-3/2
φ''(kx0) = V(zm/cos3θd)/|ω| = V(rd/cos2θd)/|ω|
∫eiφ(kx)dkx = Vcosθdrdeiπ/4|ω|½e-iωtd(x)

p(xm,zm) = V∫dxcosθdrd∫dωeiπ/4|ω|½e-iωtd(x)P(x,ω)

La formule de Green fournit directement l'expression de la migration de Kirchhoff en 3D. La solution élémentaire G(x,xm) de l'équation des ondes pour une source impulsionnelle δ(x-xm) placée en xm est l'onde sphérique issue de xm et enregistrée à la surface en x. Le champ d'onde réfléchi P(x) se propage selon l'équation des ondes homogène. La formule de Green permet d'exprimer P(xm) en fonction de l'intégrale de surface de P(x,y,z=0) et (∂/∂z)P(x,y,z=0). A cause de la surface libre en z=0, la solution élémentaire est la différence entre les ondes sphériques provenant de la source impulsionnelle en z>0 et virtuelle en z<0. Elle s'annule en z=0, ce qui permet d'exprimer P(xm) uniquement en fonction de P(x,y,z=0), seul enregistré, et gradG (l'angle entre gradG et l'axe vertical est l'angle de diffraction θd). Pour la migration, il faut rétropropager l'onde de x à xm le long du cône de diffraction en remontant dans le temps de td à t=0 (par un déphasage de -ωtd). On trouve que la quantité à sommer le long de l'hyperboloïde de diffraction est ∂P/∂t avec un terme de pondération en amplitude.

ΔP+k2P(x,ω) = 0
ΔG+k2G(x,xm,ω) = δ(x-xm)
rd = |x-xm| = Vtd , erd = (x-xm)/rd
G(x,xm,ω) = -eikrd/rd
gradG(x,xm,ω) = erd-(ik-1/rd)eikrd/rd
gradG(x,xm,ω) ≈ -ikerdeikrd/rd pour krd >> 1
∫∫∫P(x,ω)δ(x-xm)dx = ∫∫∫(PΔG-GΔP)dx = ∫∫∫div(PgradG-GgradP)dx
P(xm,ω) = ∫∫(PgradG-GgradP)(x,y,z=0,ω).dS
G(x,xm,ω) = -eikrd(zm)/rd(zm)+eikrd(-zm)/rd(-zm)
P(xm,ω) = ∫∫(PgradG)(x,y,z=0,ω).dS = -2∫∫ikP(x,y,z=0,ω)eikrdcosθd/rddxdy
p(xm, t=0) = (2/V)∫∫dxdycosθd/rd∫-iωP(x,y,z=0,ω)e-iωtd
p(xm, t=0) = (2/V)∫∫(∂/∂t)P(x,y,z=0,td)cosθd/rddxdy

Schneider, W.A., 1978, Integral formulation for migration in two and three dimensions.

Migration par point de tir

A déport non-nul, les points diffractants dans le plan (x,z) sont atteints par l'onde issue de la source au temps direct source-point diffractant tsd. Le temps d'arrivée en z = 0 est la somme de ce temps et du temps diffracté depuis le point diffractant. La migration peut être faite par sommation des amplitudes diffractées le long de cette trajectoire de diffraction.
Pour un couple source-récepteur, la surface dans le plan (x,z) correspondant à un temps t est une demi-ellipse (une ellipse est engendrée avec une corde de longueur constante fixée en deux points). Un échantillon au temps t peut donc être placé le long de cette demi-ellipse. En superposant les demi-ellipses correspondant aux différents échantillons, on réalise la migration.

cône de diffraction - cone.m hyperbole de diffraction en z=0 position du point diffractant à t = tsd
temps direct source-point diffractant : tsd = (1/V)((xm-xs)2+zm2)½ td = tsd + (1/V)((x-xm)2+zm2)½ x = xm , z = zm
migration de p(x,t) par sommation des diffractions émises par chacun des points du sous-sol au temps tsd sommation sur la courbe de diffraction image de M = onde(xm,zm,tsd)
algorithme de sommation sur trajectoire de diffraction dans le plan xt : ∀(xm,zm) S = ∑xp(x,td(x,xs,xm,zm)) p(xm,zm) = S

Voici un exemple de migration d'un point de tir pour un point diffractant, des points diffractants alignés sur un segment de droite et un synclinal et pour un réflecteur penté :

migration point de tir - point diffractant - migration point de tir - segment horizontal - migration point de tir - segment penté - migration point de tir - synclinal - migration point de tir - réflecteur penté - migsdc.m

Dip Move-Out

Pour un réflecteur plan horizontal, les temps de réflexion t enregistrés avec un déport non nul sont ramenés à déport nul par une correction dynamique (NMO) t→tn. Pour un réflecteur plan de pendage α, l'hyperbole de réflexion observée sur un point de tir est décalée à la verticale de la source virtuelle. Sa courbure et ses asymptotes correspondent à la vitesse de propagation V du milieu. Sur un point milieu commun, l'hyperbole de réflexion est centrée sur la valeur de déport nulle mais sa courbure et ses asymptotes sont fonction de V/cosα. Dans le triangle source virtuelle Sv - point d'intersection O du réflecteur avec la surface - récepteur G, on a :

(SvG)2 = (OSv)2+(OG)2-2(OSv)(OG)cos(2α) (Vt)2 = s2+g2-2sgcos(2α) point milieu fixé y0 , t0 = 2y0sinα/V
s = y-h ; g = y+h (Vt)2 = (y-h)2+(y+h)2-2(y2-h2)(cos2α-sin2α) = 4y2sin2α+4h2cos2α t2 = (4y02sin2α+4h2cos2α)/V2 = t02+4h2cos2α/V2

Voici des exemples de trajectoires de réflexion en point de tir, points milieux et déport nul pour un milieu comportant un réflecteur penté et un réflecteur horizontal se croisant en un point :
Tir, point-milieu, déport nul - Tir reflecteur horizontal-penté - Point-milieu reflecteur horizontal-penté

Les points de réflexion sur un réflecteur penté pour un point milieu commun ne sont pas confondus comme dans le cas où le pendage est nul. Pour s'affranchir de ce phénomène il faut effectuer une migration partielle appelée Dip Move Out (DMO). Si on ne fait pas d'hypothèse sur le pendage, une impulsion sur une section à déport constant h correspond à un réflecteur elliptique dont les foyers sont la source et le récepteur. Les temps à déport nul t0 correspondant aux différents pendages sont obtenus en construisant les rayons normaux au réflecteur elliptique. Ceux-ci intersectent l'axe z=0 en des points-milieux y0 différents de celui de l'impulsion initiale. L'ellipse DMO est la courbe t0(y0) : en reportant l'amplitude de chaque échantillon de la section à déport constant le long de l' ellipse DMO associée et en sommant, on obtient la section temps à déport nul quelques soient les pendages.

correction NMO ellipse de migration E pour une impulsion en y=0 à t intersection rayon normal à E avec z=0 longueur rayon normal entre point courant M de E et (y0,0) correction DMO
dmo.m
tn2 = t2-4h2/V2 ym2/t2+zm2/tn2 = V2/4 y0 = ym(1-tn2/t2) V2t02/4 = (ym-y0)2+zm2 t02/tn2+y02/h2=1

Voici des figures illustrant l'opérateur DMO t(t0, x0) pour un demi-déport h = 1000 m dans le cas où V = V0 = cte et dans le cas où V(z) = V0+az, en distinguant la partie de l'opérateur correspondant aux rayons réfléchis vers le haut et celle aux rayons réfléchis vers le bas (voir paragraphe suivant) :

correction DMOV0 - correction DMOVZ - correction DMOVZ
dmovz.m

Liner, C.L., 1999, Concepts of normal and dip moveout
Artley C. and D. Hale, 1994, Dip-moveout processing for depth-variable velocity
Black J.L., K.L. Schleicher and L. Zhang, 1993, True-amplitude imaging and dip moveout
Deregowski, S.M. and Rocca F., 1981, Geometrical optics and wave theory of constant offset sections in layered media: Geophys. Prosp., 29, 374-406.

Relation temps-profondeur pour V(z)

Quand la vitesse dépend de la profondeur, la relation temps de propagation verticale - profondeur n'est pas linéaire. De même, la relation fréquence - longueur d'onde verticale dépend de la profondeur. L'analyse des temps de propagation des hyperboles de réflexion sur les points milieux fournit VRMS(t). On en déduit V(t) par dérivation (formule de Dix). t(z) peut être mesuré directement si l'on dispose d'un puits, ce qui permet de "caler" en profondeur les réflexions au niveau du puits.

Par exemple, pour une augmentation linéaire de la vitesse avec la profondeur, on a :

V(z) = V0+az dt = dz/V(z)
ω = V(z)kz(z)
t(z) = ∫(dz/V(z)) = (1/a)Log(1+az/V0) VRMS2(z) = (1/t(z))∫(V(z)dz) = (z/t(z))(V0+az/2) V(z) = d(t(z)VRMS2(z))/dz
V(t) = V0eat z(t) = (V0/a)(eat-1) VRMS2(t) = (1/t)∫(V2(t)dt) = (V02/(2at))((e2at-1) V2(t) = d(tVRMS2(t))/dt

temps-profondeur

Migration pour V(z) par tracé de rayons ou "phase-shift"

Quand la vitesse ne dépend que de la profondeur, la vitesse apparente horizontale le long d'un rayon est constante. Pour calculer les courbes de diffraction exactement, il faut tracer des rayons à partir de chaque profondeur en utilisant la loi de Snell. L'augmentation de la vitesse avec la profondeur est favorable pour l'imagerie des forts pendages. En effet, l'angle d'incidence diminuant vers la surface, le rayon concave parcourt une distance horizontale inférieure à celle dans un milieu de vitesse constante. Il est donc moins sensible aux variations latérales de vitesse. De plus, il est possible d'utiliser les rayons diffractés vers la bas pour imager les réflecteurs pentés par en dessous (par exemple les flancs des dômes de sel). Pour une faible ouverture angulaire, l' approximation VRMS est utilisable.

migration VZ - migvz.m
migration VZ - migsvz.m

La loi de Snell s'exprime exactement dans le domaine de Fourier en exprimant kz(kx,ω,z) en fonction de V(z) : la migration par "phase-shift" est une méthode itérative en z basée sur cette propriété.

migration "phase-shift" kz(jΔz) = -ω(1/V(jΔz)2-kx22)½ = -ωcosθ(jΔz)/V(jΔz)) P(kx,ω,nΔz) = P(kx,ω,z=0)∏j(eiΔzkz(jΔz)) = P(kx,ω,z=0)eiΔz∑j(kz(jΔz))

Voici l'illustration des étapes initiale, intermédiaire et finale de la migration phase-shift pour une vitesse constante et une vitesse V(z) :

migration phase-shift V0 - migration phase-shift VZ

Pour modéliser une section sismique à déport nul lorsque la vitesse dépend de la profondeur, on peut utiliser le même principe. La TF du modèle p(x,z,t=0) donne P(kx,z,t=0). On remonte vers la surface depuis la profondeur maximum par des déphasages successifs kz(kx,ω,z)Δz et en ajoutant à chaque pas en profondeur la contribution de P(kx,z,t=0) à cette profondeur.

modélisation "phase-shift" kz(jΔz) = ω(1/V(jΔz)2-kx22)½ = ωcosθ(jΔz)/V(jΔz)) P(kx,ω,z-Δz) = P(kx,ω,z)eiΔzkz(z)+P(kx,t=0,z-Δz)

modélisation phase-shift VZ

Si la méthode phase-shift traite correctement la phase quand la vitesse dépend de la profondeur, elle ignore la variation d'amplitude liée à la réfraction des rayons. La méthode WKBJ de résolution de l'équation d'onde pour V(z) donne l'effet sur l'amplitude.

2P/∂z2+kz(kx,ω,z)2P(kx,ω,z) = 0 P(kx,ω,z) = A(z)ei∫kz(z)dz 2kzdA/dz+Adkz/dz = 0 A(z)/A(z0) = (kz(z0)/kz(z) )½

Voici la correspondance entre un réflecteur en forme de synclinal et la réflexion correspondante pour un milieu comportant une couche de vitesse V0 = 500 m/s sur un milieu de vitesse V1 = 1000 m/s et des temps simples de propagation. Les segments courbés du synclinal sont des arcs de cercle.
synclinalV01

Gazdag, J., 1978, Wave equation migration with the phase-shift method
Hale D., N. R. Hill, and J. Stefani, 1992, Imaging salt with turning seismic waves montrent comment utiliser les ondes réfractées par un gradient de vitesse pour imager les réflecteurs par dessous.
Zhang, J. and McMechan, G. A., 1997, Turning wave migration by horizontal extrapolation montrent comment combiner prolongement vertical et horizontal du champ d'onde pour imager tous les pendages.
Margrave G.F., 2001, Direct Fourier migration for vertical velocity variations généralise la migration de Stolt au cas où V=V(z)

Migration pour V(x,y,z) par tracé du rayon image ou "phase-shift" avec correction

Si les variations latérales de vitesse sont petites par rapport aux variations avec la profondeur, on peut se contenter d'apporter des corrections aux méthodes précédentes.

Le temps minimum entre un point diffractant et la surface correspond au rayon image, orthogonal à la surface (puisque le front d'onde atteint en premier la surface selon une incidence nulle). Si on effectue une sommation sur trajectoire de diffraction hyperbolique (sans tenir compte des variations latérales de vitesse), la somme doit être placée en profondeur non pas à la verticale du point d'émergence du rayon image, mais le long du rayon image au temps correspondant au temps de l'apex de la trajectoire de diffraction. Ceci nécessite de tracer uniquement les rayons image.

Black J.L. and M. A. Brzostowski, 1994, Systematics of time-migration errors

On peut aussi utiliser la méthode de phase-shift avec une vitesse moyenne Vm(z) et effectuer une correction dans le domaine spatial tenant compte approximativement de l'écart V(x,z)-Vm(z). On fait un développement au premier ordre de kz et on applique un déphasage vertical variant latéralement en plus du terme tenant compte de la propagation sous tous les angles à la vitesse moyenne.

migration "phase-shift" avec correction kz ≈ -ω/V(x,z)+ω/Vm(z)-ω(1/Vm(z)2-kx22)½ p(x,z+Δz) = ∫dωe-iω(1/V(x,z)-1/Vm(z))Δz{∫P(kx,ω,z)e-iω(1/Vm(z)2-kx22)½Δzeikxxdkx}

Voici la correspondance entre un réflecteur en forme de synclinal et la réflexion correspondante pour un milieu comportant une couche de vitesse V0 = 500 m/s sur un milieu de vitesse V1 = 1000 m/s et des temps simples de propagation. L'interface a un pendage de 10°. Les segments courbés du synclinal sont des arcs de cercle.
synclinalVxz

Opérateur de différences finies pour l'équation iconale

Si les trajectoires de diffraction sont trop déformées par les variations latérales de vitesse pour qu'une sommation sur une trajectoire hyperbolique soit constructive, il est nécessaire de les déterminer avant d'effectuer la sommation. Ceci nécessite de déterminer les temps de trajet à partir de chaque point diffractant. La résolution numérique de l'équation iconale par différences finies depuis chaque position source permet de déterminer les temps directs entre tout point de la surface et tous les points en profondeur.

équation iconale différences finies centrées opérateur pour Δx = Δz = a onde conique
(∂T/∂x)2+(∂T/∂z)2 = 1/V(x,z)2 (∂T/∂x) = (T(x+Δx,z)-T(x,z)+T(x+Δx,z+Δz)-T(x,z+Δz))/2Δx
(∂T/∂z) = (T(x,z+Δz)-T(x,z)+T(x+Δx,z+Δz)-T(x+Δx,z))/2Δz
T(x+Δx,z+Δz) = T(x,z)+(2a2/V(x,z)2-(T(x+Δx,z)-T(x,z+Δz))2)½ T(x+Δx,z+Δz) = min(T(x+Δx,z),T(x,z+Δz))+a/V(x,z)

Voici des figures illustrant les temps calculés par différences finies pour une source ponctuelle dans un milieu de vitesse constante ou à gradient vertical et latéral sans interface.

iconale V0 - iconale V0 - iconale Vz - iconale Vz - iconaleVxz -
iconal.m

Noble, M., A. Gesret and N. Belayouni, 2014, Accurate 3-D finite difference computation of traveltimes in strongly heterogeneous media
Mo L.W. and J. M. Harris, 2002, Finite-difference calculation of direct-arrival traveltimes using the eikonal equation

Equations d'onde paraxiales

Une autre façon de tenir compte des déformations des fronts d'onde par les variations de vitesse est d'utiliser une résolution numérique d'une équation d'ondes. On obtient des équations d'onde adaptées à la migration (propagation d'ondes montantes dans un domaine d'ouverture angulaire limité autour de l'incidence verticale : kx ≈ 0 , kz ≈ -ω/V) en faisant des développements limités de la relation de dispersion des ondes montantes valables pour des domaines restreints de pendage. Ces équations, obtenues dans le domaine de Fourier, peuvent être écrites dans le domaine spatial et temporel en identifiant les nombre d'ondes et pulsations avec des dérivées partielles. On obtient ainsi des équations d'onde paraxiales que l'on utilise pour prolonger le champ d'onde en profondeur par des méthodes de différences finies.

relation dispersion kz = -ω(1/V2-kx22)½ terme de diffraction éq. paraxiale (x,z,ω) éq. paraxiale (x,z,t)
développement 15° kz ≈ -ω/V(1-V2kx2/2ω2) (ikz)(-iω) = Vkx2/2 ∂P/∂z = (V(x,z)/2iω)∂2P/∂x2 2p/∂z∂t = -(V(x,z)/2)∂2p/∂x2
développement 45° kz ≈ -ω/V(1-(V2kx22)/(2-V2kx2/2ω2)) 2-V2kx2/4)(ikz) = (iω)Vkx2/2 2+(V2/4)∂2/∂x2)∂P/∂z = -(iωV/2)∂2P/∂x2 (∂3/∂t2∂z-(V2/4)∂3/∂x2∂z)p = -(V/2)∂3p/∂t∂x2

Les relations de dispersion des équations paraxiales ne sont pas isotropes : la vitesse dépend de la direction de propagation. La réponse impulsionnelle de la migration n'est plus un demi-cercle. Pour le développement à 15°, c'est une ellipse. On le montre à partir de l'expression des composantes de la vitesse de groupe pour la relation de dispersion.

relation dispersion 15° phase stationnaire vitesse de groupe = grad(ω(kx,kz)) équation du front d'onde elliptique
kz ≈ -ω/V(1-V2kx2/2ω2) x = (∂ω/∂kx)t
z = (∂ω/∂kz)t
∂ω/∂kx = (V2kx/ω)/(1+V2kx2/2ω2)
∂ω/∂kz = -V/(1+V2kx2/2ω2)
x2/2+(z+Vt/2)2 = V2t2/4

Ristow D. and T. Ruhl, 1994, Fourier finite-difference migration
Gray, S.H., J. Etgen, J. Dellinger and D. Whitmore, 2001, Seismic migration problems and solutions font une présentation générale des méthodes de migration et de leur utilisation pour l'estimation des vitesses.

Opérateurs de différences finies pour l'équation paraxiale à 15°

Les opérateurs de différences finies de l'équation paraxiale à 15° s'obtiennent sous forme explicite ou implicite dans les différents domaines. La formulation explicite est instable. La formulation implicite, toujours stable, nécessite la résolution d'un système d'équations linéaires à matrice tridiagonale à chaque pas d'extrapolation (voir algorithme de Thomas). Dans le domaine (x,z,ω), on a :

coefficients des opérateurs de différences finies de l'équation paraxiale à 15° ∑aiPi = 0 analyse stabilité pour P = eikxx
explicite (x,z,ω) P(x,z+Δz) = P(x,z)+(V/2iω)(Δz/Δx2)(P(x-Δx,z)-2P(x,z)+P(x+Δx,z)) P(z+Δz) = P(z)[1-4αsin2(kxΔx/2)]

|P(z+Δz)| > |P(z)|

α = (V/2iω)(Δz/Δx2) P(x-Δx, ., ω) P(x, ., ω) P(x+Δx, ., ω)
P(., z, ω)
P(., z+Δz, ω)

0
2α-1
1

0
implicite (x,z,ω) P(x,z+Δz) = P(x,z)+(V/2iω)(Δz/2Δx2)(P(x-Δx,z)-2P(x,z)+P(x+Δx,z)+P(x-Δx,z+Δz)-2P(x,z+Δz)+P(x+Δx,z+Δz)) P(z+Δz)[1+4αsin2(kxΔx/2)] = P(z)[1-4αsin2(kxΔx/2)]

|P(z+Δz)| = |P(z)|

α = (V/2iω)(Δz/2Δx2) P(x-Δx, ., ω) P(x, ., ω) P(x+Δx, ., ω)
P(., z, ω)
P(., z+Δz, ω)

2α-1
2α+1

Voici des figures illustrant le résultat de la migration par différences finies 15° (x, z, ω) implicite pour des réflecteurs plans de pendage croissants (la solution géométrique est donnée par la ligne rouge) et les réponses impulsionnelles pour une vitesse constante V = 2000 m/s :

migration dif. finies - migration dif. finies -
dif15.m

L'opérateur implicite dans le domaine (x,z,t) comporte 2x2x3 = 12 termes. On l'obtient à partir de celui à 2x2 termes dans le domaine (kx,z,t) en développant la dérivée seconde en x sur 3 termes.

implicite (kx,z,t) P(z+Δz,t+Δt)-P(z+Δz,t)-P(z,t+Δt)+P(z,t) = (ΔtΔz/4)(Vkx2/2)(P(z+Δz,t+Δt)+P(z+Δz,t)+P(z,t+Δt)+P(z,t))
α = (Vkx2/2)(ΔtΔz/4) P(kx, . , .)
P(kx, z , t)
P(kx, z , t+Δt)
P(kx, z+Δz , t)
P(kx, z+Δz , t+Δt)
α-1
α+1
α+1
α-1
implicite (x,z,t)
α = -(V/2)(ΔtΔz/4Δx2) p(x-Δx, ., .) p(x, ., .) p(x+Δx, ., .)
p(., z , t)
p(., z , t+Δt)
p(., z+Δz , t)
p(., z+Δz , t+Δt)
α
α
α
α
-2α-1
-2α+1
-2α+1
-2α-1
α
α
α
α

Si on représente des enregistrements de sismique de puits p(z,t) en fonction de z et t' = t+z/V (temps dit réduit), on fait apparaitre les ondes montant verticalement à la vitesse V (dt/dz = -1/V) à temps constant (celui en z = 0). La longueur d'onde verticale des ondes montantes devient infinie (k'z = 0) sans modifier la période. Ce changement de variable est aussi utile dans la résolution numérique de l'équation paraxiale : en temps retardé, le champ d'onde change peu avec Δz, et l'extrapolation vers le bas est facilitée. Pour une vitesse variable, on utilise la vitesse moyenne à la profondeur z , VM(z).

temps réduit effet sur les pentes et sur les dérivées partielles kz ≈ -ω/V , k'z ≈ 0 éq. paraxiale 15°
t' = t+z/V dt'/dz = dt/dz + 1/V
k'z/ω = kz/ω+1/V
k'z = kz+ω/V k'z = Vkx2/2ω
t' = t+∫dz/VM(z) P(z,t) = P'(z,t'(z,t))
∂P/∂z = ∂P'/∂z+(1/VM(z))∂P'/∂t' , ∂P/∂t = ∂P'/∂t'
kz = k'z-ω/VM(z)
k'z = kz+ω/VM(z) k'z = -ω(1/V(x,z)-1/VM(z))+Vkx2/2ω

Pour éviter les réflexions parasites sur les limites latérales de la grille, il faut utiliser des conditions aux limites "anti-reflets". Dans le cas de la résolution de l'équation d'ondes 1D par différences finies, les réflexions peuvent être supprimées en utilisant l'équation des ondes se propageant dans un seul sens sur les bords de la grille. Pour l'équation des ondes 2D, on peut utiliser des équations paraxiales pour les ondes allant vers les x croissants ou décroissants. Pour l'équation des ondes montantes, Clayton et Engquist (1980) ont proposé une relation hyperbolique entre kx/ω et kz/ω approximant le quart de cercle de la relation de dispersion correspondant aux ondes allant vers les x croissants ou décroissants. Les 3 paramètres a,b,c de cette relation sont obtenus en écrivant l'égalité des relations de dispersion pour 3 angles de propagation, par exemple θ = 30°, 60° et 90°.

eq. ondes vers gauche eq. ondes eq. ondes vers droite
1D (1/V)∂P/∂t = ∂P/∂x
-ω/V = kx
(1/V2)∂2P/∂t2 = ∂2P/∂x2
(ω/V)2 = kx2
(-1/V)∂P/∂t = ∂P/∂x
ω/V = kx
2D kx ≈ -ω/V(1-V2kz2/2ω2) (ω/V)2 = kx2+kz2 kx ≈ ω/V(1-V2kz2/2ω2)
paraxiale Vkz/ω = -(a-bVkx/ω)/(1-cVkx/ω)
cosθi = (a+bsinθi)/(1+csinθi)
kz = -ω/V(1-V2kz22)½ Vkz/ω = -(a-bVkx/ω)/(1-cVkx/ω)
cosθi = (a-bsinθi)/(1-csinθi)

Clayton R.W. and B. Engquist , 1980, Absorbing boundary conditions for wave-equation migration

Une autre méthodes consiste à entourer la grille avec des couches absorbantes sans contraste d'impédance avec le milieu de propagation. La méthode PML (Perfectly Matched Layer) introduit l'atténuation uniquement dans la direction perpendiculaire à l'interface. Le nombre d'onde dans la direction parallèle à l'interface n'est pas modifié. Il y a transmission totale dans la couche absorbante quelle que soit la période et l'incidence. L'atténuation est introduite par un changement de variable dans la couche. Pour une interface verticale en x = 0 séparant le mileu de propagation en x<0 d'un demi-espace absorbant en x>0, le calcul des coefficients de réflexion-transmission pour une onde plane SH montre qu'il y a bien transmission totale.

interface en x = 0 eq. paraxiale déplacement et contrainte continuité en x = 0
x < 0 x' = x kz = Vkx2/2ω uy(x) = (eikxx+Re-ikxx)eikzz
σxy(x) = μikx(eikxx-Re-ikxx)eikzz
1+R = T
(1-R) = T
R = 0 , T = 1
x > 0 x' = x(1+iα/ω) , α>0 k'x = kx/(1+iα/ω)
kz = Vk'x2/2ω
uy(x') = Teikxx'eikzz
σxy(x') = (μikx/(1+iα/ω))Teikxx'eikzz

Collino, F, 1997, Perfectly matched absorbing layers for the paraxial equations
Pan, G., A. Abubakar, T. M. Habashy, (2012), An effective perfectly matched layer design for acoustic fourth-order frequency-domain finite-difference scheme

Renversement temporel

Le prolongement vers le bas du champ d'onde construit le cône de diffraction P(x,z,t) en reconstituant les hyperboles dans les plans P(x,z = iΔz,t) par le calcul de kz ou ∂P/∂z. Alternativement, il peut être construit en parcourant l'axe des temps en sens inverse pour reconstituer les fronts d'onde dans les plans P(x,z,t = T-iΔt), T étant le temps maximum d'enregistrement. Dans les deux cas, la pointe du cône (focalisation de l'onde diffractée) est atteinte en P(x,z,t=0).

ω/V = -kz(1+kx2/kz2 1/V∂P/∂t = ikz(1+kx2/kz2)½P(kx,kz,t)

Baysal, E., D.D. Kosloff, and J. W. C. Sherwood , 1983, Reverse time migration

Géométrie de réflexion et équations de prolongement à déport non-nul

coordonnées sources - récepteurs coordonnées points-milieux - déports
θs = θ-α θg = θ+α y = (s+g)/2 h = (g-s)/2
∂t/∂s = -sinθs/V ∂t/∂g = sinθg/V ∂t/∂y = ∂t/∂s+∂t/∂g = 2sinαcosθ/V ∂t/∂h = -∂t/∂s+∂t/∂g = 2cosαsinθ/V
∂t/∂zs = -cosθs/V = -(1-(V∂t/∂s)2)½/V ∂t/∂zg = -cosθg/V = -(1-(V∂t/∂g)2)½/V ∂y/∂z = -1/tgα ∂h/∂z = -1/tgθ
∂t/∂z = ∂t/∂zs+∂t/∂zg = -2cosαcosθ/V = -((1-(V∂t/∂s)2)½+(1-(V∂t/∂g)2)½)/V ∂t/∂z = -((1-(V(∂t/∂y-∂t/∂h)/2)2)½+(1-(V(∂t/∂y+∂t/∂h)/2)2)½)/V
kz/ω = -((1-(Vks/ω)2)½+(1-(Vkg/ω)2)½)/V kz/ω = -((1-(V(ky-kh)/2ω)2)½+(1-(V(ky+kh)/2ω)2)½)/V
∂P/∂z = (-iω/V)((1-(Vks/ω)2)½+(1-(Vkg/ω)2)½)P(ks,kg,z,ω) ∂P/∂z = (-iω/V)((1-(V(ky-kh)/2ω)2)½+(1-(V(ky+kh)/2ω)2)½)P(ky,kh,z,ω)
∂P/∂z = (( (-iω/V(s))2-∂2/∂s2)½+( (-iω/V(g))2-∂2/∂g2)½)P(s,g,z,ω)

Imagerie à déport non-nul

L'image d'un point du sous-sol est obtenue en faisant le rapport dans le domaine fréquentiel entre les ondes montantes M (le champ d'onde réfléchi enregistré au récepteur) et les ondes descendantes D (le champ d'onde émis par les sources) après prolongation vers le bas de ces champs, ce qui fournit le coefficient de réflexion R(x,z,ω). Dans le domaine temporel, cela correspond à une déconvolution du champ montant par le champ descendant (comme en sismique de puits). La déconvolution étant potentiellement instable, on la remplace en pratique par la corrélation entre les champs montant et descendant, l'effet sur la phase étant identique. L'image I(x,z) d'un point peut être obtenue avec différentes formules en sommant sur les fréquences et les positions de sources et en stabilisant les divisions par l'illumination et un ajout de bruit blanc.

  déconvolution déconvolution stabilisée corrélation corrélation normalisée par l'illumination
R(x,z,ω) M(x,z,ω)/D(x,z,ω) M(x,z,ω)D*(x,z,ω)/(D(x,z,ω)D*(x,z,ω)) M(x,z,ω)D*(x,z,ω)
I(x,z) sω(M(x,z,s,ω)/D(x,z,s,ω)) sω(M(x,z,s,ω)D*(x,z,s,ω))/(D(x,z,s,ω)D*(x,z,s,ω)+ε2)) sωM(x,z,s,ω)D*(x,z,s,ω) s(∑ω(M(x,z,s,ω)D*(x,z,s,ω))/∑ω(D(x,z,s,ω)D*(x,z,s,ω)))

Voici l'exemple de la corrélation entre une onde harmonique sphérique émise par une source ponctuelle en surface (onde descendante) avec l'onde sphérique émise par la source virtuelle correspondant à un réflecteur plan (onde réfléchie montante rétropropagée). Le réflecteur (ligne blanche) coïncide avec un maximum de la corrélation : temps source-réflecteur = temps source virtuelle -réflecteur. Les autres maximums de corrélation correspondent à la périodicité de l'onde harmonique : phases égales modulo 2π.
corrélation montante-descendante

Schleicher, J., J. C. Costa, and A. Novais, 2008, A comparison of imaging conditions for wave-equation shot-profile migration

Imagerie d'un réfracteur par les ondes coniques

La migration peut être adaptée au cas de l'onde conique le long d'un réfracteur de pendage α entre deux milieux de vitesse V1 et V2 (sinθc = V1/V2). Le lieu des points tels que la somme du temps d'arrivée de l'onde transmise en xm le long du réfracteur et du temps diffracté depuis les points d'abscisse xm vers la surface soit constante pour une position de récepteur xg correspond à la position migrée pour ce temps et cette trace. Le front d'onde transmis étant orthogonal au réfracteur, le temps d'arrivée sur un plan perpendiculaire au réfracteur est constant. Si le pendage est faible, ce plan peut être assimilé au plan vertical. d étant la distance source-réfracteur, on a :

temps d'arrivée en xm le long du réfracteur temps diffracté depuis le réfracteur migration par superposition de fronts d'onde
t0 = d/cosθc/V1+(xm-dsin(θc-α)/cosθc)/cosα/V1 td = ((xg-xm)2+zm2)½/V1 p(xm, zm) = p(xm, zm)+p(x, (t0+td)(x, xm, zm))

migration onde conique -
migrefr.m

Transformée τ-p

directe inverse
v(τ,p) = ∫u(t=τ+px, x)dx = ∫∫U(ω,x)e-iω(τ+px)dωdx = ∫U(ω,kx=ωp)e-iωτ w(t,x) = ∫v(τ=t-px,p)dp = ∫∫U(ω,kx=ωp)e-iω(t-px)dωdp = ∫∫(U(ω,kx)/|ω|)e-iωteikxxdωdkx = u(t,x)∗ = ∫(U(ω,x)/|ω|)e-iωt
U(ω,x) = |ω|W(ω,x)

Séparation des ondes montantes et descendantes dans l'eau

L'enregistrement simultané à une profondeur h dans l'eau de la pression P et de la composante verticale Vz de la vitesse de déplacement des particules v fournit des informations spectrales complémentaires. En effet, les signes du coefficient de réflexion sur la surface libre étant opposés pour la pression et la vitesse de déplacement, un trou dans le spectre de pression dus aux interférences entre les ondes montantes et les réflexions sur la surface libre (2kzh+π = (2n+1)π ou f = nV/(2hcosθ)) correspond à un maximum pour la vitesse de déplacement et inversement. De plus cela permet de séparer les ondes montantes et descendantes au niveau des récepteurs. En supposant que le coefficient de réflexion pour la pression sur la surface libre est R ( = -1 s'il fait beau, mais ≠ -1 s'il y a des vagues), on a :

eq. ondes et potentiel de vitesse vitesse de déplacement pression
ρ∂v/∂t = -gradp v = gradφ p = -ρ∂φ/∂t
Φ = (Φmd)e-iωt = (e-ikzz+Reikzz)e-iωt Vz = ikz (-e-ikzz+Reikzz)e-iωt P = iωρ(e-ikzz+Reikzz)e-iωt
potentiel montant et descendant Φm = (1/2)(P/(iωρ)-Vz/(ikz)) Φd = (1/2)(P/(iωρ)+Vz/(ikz))

Bibliographie et liens