JMM @ EOST

Sismique : Ondes - Imagerie - Traitement - Propagation - Applications

Ondes élastiques : Caractéristiques - Elasticité - Equations d'ondes - Ondes planes - Ondes sphériques
Sismique géométrique : Fronts d'ondes et rayons - Source ponctuelle - V(z) - V(x,y,z)
Réflexion-transmission : Continuité - SH - P/SV surface libre - P/SV interface - Approximations
Ondes guidées : Rayleigh par surface libre et Stoneley par interface - Love et Rayleigh dans couche - Vitesse de groupe - Love dans stratification
Ondes élastiques + : Ondes cylindriques - Ondes P et S dues à une force ponctuelle - Anisotropie - Atténuation - Saturation en fluide
Réflexion-transmission + : Gradient vitesse - Couche et stratification - Monitoring sismique - Stratification cyclique - Relations réflexion-transmission - Surface ondulée
Documents : Figures, programmes, animations - Sujets exams - Sujets TD - Livres et liens
Programmes pour contrôle continu : ondePS.m - hodostrat.m - refltransPSV.m

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 (ajouter au préalable l'URL de la page dans l'onglet sécurité du panneau de configuration java pour permettre leur ouverture). Les programmes Matlab écrits pour réaliser les figures sont accessibles en cliquant sur le nom du programme (prog.m) sous l'icône de la figure. Pour imprimer entièrement les formules et tableaux, utiliser le format paysage (feuille en travers).

Introduction

Les ondes sismiques sont des ondes élastiques qui se propagent dans la Terre. Pour se familiariser avec la façon dont elles se propagent, il est conseillé de visiter les sites suivants:

Russel D., Acoustics Animations : Longitudinal and Transverse Wave Motion
Le site illustre la propagation 1D d'ondes longitudinales et transverses harmoniques ainsi que celle des ondes de Rayleigh guidées par une surface libre.

Baker S.A. and Wysession M.E., Seismic shear wave propagation : a synthetic animation
On y trouve des animations de la propagation et de la réflexion de fronts d'onde dans le manteau terrestre, et des exemples de tracés de rayons et d'enregistrements.

Tanimoto T., Southern California wavefield reconstruction
Visualisation de la propagation de champs d'ondes enregistrés par un réseau dense de stations sismologiques.

Caractéristiques

Une onde sismique se propage en ondulant (λ = VT), déplace les particules autour de leur position d'équilibre, déforme et contraint le milieu élastique de propagation. Les grandeurs physiques intervenant dans la propagation et leur ordre de grandeur sont les suivants (les vecteurs sont en caractères gras).

grandeuruniténotationcomposantesordres de grandeur
longueur d'onde m λ . 10-3 m (grains dans les roches) - 107 m (taille de la Terre)
période
fréquence
pulsation
s
Hz
rad/s
T
f = 1/T
ω = 2πf
. 106 Hz (ultrasons sur échantillons) - 10-3 Hz (vibrations propres de la Terre)
vitesses de propagation
ondes P, S
m/s VP
VS
stratification sedimentaire 102 m/s (VS dans les sols) - 104 m/s (VP dans le manteau)
VP/VS : 1.7 dans les roches consolidées, >10 dans les sols peu consolidés saturés en eau
densité
impédances
kg/m3
kg/m2/s
ρ
I = ρV
103 kg/m3 (eau) - 104 kg/m3 (noyau)
déplacement des particules m u ux , uy , uz U : 1 m (glissement sur faille) - 10-8m (bruit ambiant)
vitesse des particules m/s u/∂t ∂ux/∂t , ∂uy/∂t , ∂uz/∂t U/T - voir les définitions des différentes échelles de magnitudes à partir du rapport U/T
accélération des particules m/s2 2u/∂t2 2ux/∂t2 , ∂2uy/∂t2 , ∂2uz/∂t2 U/T2 : 10 m/s2≈1g (épicentres des grands tremblements de terre) - 10-11m/s2 (vibrations propres de la Terre)
déformations . ε εxx , εyy , εzz
εxy , εyz , εzx
U/λ : 10-6
contraintes Pa σ σxx , σyy , σzz
σxy , σyz , σzx
ρV2U/λ = IU/T : 104 - 1 Pa (pressions enregistrées dans l'eau en sismique réflexion, f = 10-100Hz)
direction de propagation . e ex , ey , ez
sinθ , 0 , cosθ
les sources sismiques émettent tous azimuts
vecteur d'onde m-1 k = (2π/λ)e kx , ky , kz caractérise une onde plane harmonique
facteur de qualité . Q . ∞ - 10 (sédiments peu consolidés)

Elasticité linéaire isotrope

Les roches sont des ressorts. Les déformations dues aux ondes sismiques sont suffisamment faibles (sauf au voisinage des sources) pour que les contraintes soient proportionnelles aux déformations. Les coefficients de proportionnalité sont les modules élastiques. Un milieu isotrope (comportement indépendant de l'orientation d'application des contraintes) est caractérisé par deux modules élastiques. Une compression uniaxiale fournit le module d'Young et le coefficient de Poisson de la roche. Les coefficients de Lamé λ et μ sont utilisés pour des états de contrainte-déformation quelconques. L'incompressibilité K est le rapport entre la pression moyenne et la déformation volumique.

. module élastique dilatation cisaillement variation volume
déformation . εxx = ∂ux/∂x εxy = (∂ux/∂y+∂uy/∂x)/2 divu = εxxyyzz
compression uniaxiale module Young E
coef. Poisson ν
σxx = Eεxx
εyy = -νεxx
. divu = (1-2ν)εxx
contrainte-déformation quelconque coef. Lamé λ , μ
incompressibilité K
σxx = λdivu+2μεxx σxy = 2μεxy divu = (σxxyyzz)/3K
relations entre modules E = 1010 - 1011 Pa
ν = 0.2 - 0.4 (roches) ; 0.5 (eau)
λ = Eν/((1+ν)(1-2ν)) μ = E/(2(1+ν)) K = λ+2μ/3 = E/(3(1-2ν))

Voici des valeurs de modules élastiques et densités pour l'air, l'eau et quelques minéraux :

air eau glace calcite
CaCO3
quartz
SiO2
olivine
(Mg,Fe)SiO4
K (GPa) 0.0001 2.4 8.4 70 37 130
μ (GPa) 0 0 3.6 29 44 80
ρ (kg/m3) 1 1000 920 2700 2700 3320

Lorsqu'on utilise des forages, on utilise les coordonnées cylindriques (r,θ,z). Les déformations et les contraintes s'écrivent :

déformations εrr = ∂ur/∂r εθθ = (1/r)(ur+∂uθ/∂θ) εzz = ∂uz/∂z εzr = (∂ur/∂z+∂uz/∂r)/2 ε = ((1/r)(∂ur/∂θ-uθ)+∂uθ/∂r)/2 εθz = ((1/r)∂uz/∂θ+∂uθ/∂z)/2
contraintes σrr = λdivu+2μεrr σθθ = λdivu+2μεθθ σzz = λdivu+2μεzz σzr = 2μεzr σ = 2με σθz = 2μεθz

avec divu = εrrθθzz = (1/r)(∂(rur)/∂r+∂uθ/∂θ)+∂uz/∂z

Equations des ondes élastiques en milieu homogène isotrope

Deux types d'onde peuvent se propager en milieu élastique isotrope avec des vitesses distinctes. Les ondes de compression, dites ondes P, sont polarisées longitudinalement (le déplacement des particules est orienté dans la direction de propagation). Elles sont les plus rapides. Elles sont l'équivalent en milieu élastique des ondes acoustiques (sons) se propageant dans l'air et l'eau. Les ondes de cisaillement, dites ondes S, sont polarisées transversalement ( le déplacement des particules se fait dans un plan orthogonal à la direction de propagation). Elles n'existent que dans les milieux pouvant être cisaillés élastiquement, donc pas dans les milieux acoustiques. Elles sont plus lentes que les ondes P. Pour des roches compactées (λ ≈ μ ou ν ≈ 0.25), VP/VS ≈ 1.7. Mathématiquement, la polarisation correcte des ondes s'obtient automatiquement en écrivant le déplacement des particules comme le gradient d'un potentiel scalaire φ pour les ondes P, ou comme le rotationnel d'un potentiel vecteur ψ pour les ondes S.

Ondes de compressiononde longitudinale 1Dondes acoustiquesondes P
déplacement ux(x,t) u = gradφ(x,t) u = gradφ(x,t)
déformations εxx = ∂ux/∂x divu = Δφ εxx = ∂2φ/∂x2 , εyy , εzz
εxy = ∂2φ/∂x∂y , εyz , εzx
contraintes σxx = E∂ux/∂x p = -Kdivu σxx = λΔφ+2μ∂2φ/∂x2 , σyy , σzz
σxy = 2μ∂2φ/∂x∂y , σyz , σzx
équations d'équilibre ρ∂2ux/∂t2 = ∂σxx/∂x ρ∂2u/∂t2 = -gradp ρ∂2ux/∂t2 = ∂σxx/∂x+∂σxy/∂y+∂σxz/∂z
ρ∂/∂x(∂2φ/∂t2) = (λ+2μ)∂/∂x(Δφ)
ρgrad(∂2φ/∂t2) = (λ+2μ)grad(Δφ)
équation des ondes 2ux/∂t2 = VL22ux/∂x2 2p/∂t2 = VA2Δp 2φ/∂t2 = VP2Δφ
vitesse de propagation VL2 = E/ρ VA2 = K/ρ VP2 = (λ+2μ)/ρ

Pour des ondes S se propageant dans un plan vertical (c'est le cas lorsque la vitesse de propagation est constante ou qu'elle ne dépend que de la profondeur, comme dans un milieu stratifié horizontement) , il est commode de distinguer deux types de polarisation transverses. Les ondes SH ont une polarisation horizontale perpendiculaire au plan vertical de propagation : le déplacement des particules a une seule composante non nulle. On peut donc écrire directement l'onde en fonction de cette composante. Les ondes SV ont une polarisation transverse dans le plan de propagation : les deux composantes du déplacement des particules s'obtiennent à partir d'un potentiel vecteur ayant une seule composante non-nulle, celle qui est orthogonale au plan de propagation.

Ondes de cisaillement (divu = 0)onde transverse 1Dondes SH (polarisation ⊥ plan de propagation)ondes SV(polarisation dans plan de propagation)
déplacement uy(x,t) uy(x,z,t) u(x,z,t) = rotψy(x,z,t)
ux = -∂ψy/∂z , uz = ∂ψy/∂x
déformations εxy = (1/2)∂uy/∂x εxy = (1/2)∂uy/∂x , εyz εxx = -∂2ψy/∂z∂x = -εzz
εxz = (-∂2ψy/∂z2+∂2ψy/∂x2)/2
contraintes σxy = μ∂uy/∂x σxy = μ∂uy/∂x , σyz σxx = -2μ∂2ψy/∂z∂x = -σzz
σxz = μ(-∂2ψy/∂z2+∂2ψy/∂x2)
équations d'équilibre ρ∂2uy/∂t2 = ∂σxy/∂x ρ∂2uy/∂t2 = ∂σxy/∂x+∂σyz/∂z ρ∂2ux/∂t2 = ∂σxx/∂x+∂σxz/∂z
ρ∂/∂z(∂2ψy/∂t2) = μ∂/∂z(∂2ψy/∂x2+∂2ψy/∂z2)
équation des ondes 2uy/∂t2 = VS22uy/∂x2 2uy/∂t2 = VS2Δuy 2ψy/∂t2 = VS2Δψy
vitesse de propagation VS2 = μ/ρ VS2 = μ/ρ VS2 = μ/ρ

Voici des valeurs de vitesses de propagation, coefficients de Poisson et densités pour quelques roches :

grés mal consolidé grés consolidé calcaire granite basalte granulite dunite
VP (km/s) 2.7 4.1 4.6 6.2 5.9 7.0 8.3
VS (km/s) 1.4 2.4 2.4 3.7 3.2 3.8 4.8
ν = (VP2/VS2-2)/(2VP2/VS2-2) 0.32 0.24 0.310.22 0.29 0.29 0.25
ρ (kg/m3) 2100 2400 2400 2650 2880 3000 3300

Ondes planes

Une onde plane se propage dans une direction e avec une vitesse VP ou VS. Le signe des composantes de e détermine le sens de propagation de l'onde dans les trois directions x, y et z. L'équation des ondes laisse toute liberté de choisir n'importe quelle fonction comme solution. Elle n'impose que la forme de la phase pour que l'onde se propage. Une onde plane harmonique correspond au choix d'une solution de type sinusoidal et au choix d'une période de vibration, donc d'une longueur d'onde. On vérifie que la polarisation des ondes P est longitudinale (u // e), celle des ondes S transverse (ue).

onde 1D ondes P ondes S
u→ = u(x/V-t)
u← = u(-x/V-t)
φ(e.x/VP-t)
ulongitudinal = gradφ = (e/VP)φ'
ψ(e.x/VS-t)
utransverse = rotψ = (e/VS)∧ψ'
cas harmonique : u = Ueiω(x/V-t) = Uei2π(x/λ-t/T) = Uei(kx-ωt)
u→ : k = ω/V
u← : k = -ω/V
φ = Φei(k.x-ωt) , k = (ω/VP)e
u = ikΦei(k.x-ωt)
ψ = Ψei(k.x-ωt) , k = (ω/VS)e
u = ik∧Ψei(k.x-ωt)

Pour les ondes harmoniques ∂φ/∂x = ikxφ , ∂φ/∂t = -iωφ . L'équation des ondes devient la relation de dispersion : ω2 = (kx2+ky2+kz2)V2 . C'est une sphère de rayon 1/V dans les coordonnées (kx/ω , ky/ω , kz/ω). En effet, le milieu étant homogène et isotrope, la vitesse de propagation est la même pour toutes les directions de propagation.
Une onde harmonique occupe tout l'espace avec des oscillations périodiques. Si on l'observe (ou on l'enregistre) dans sa direction de propagation e, on voit sa vraie longueur d'onde λ = 2π/k. Dans une direction différente e', on voit une
longueur d'onde apparente supérieure à la vraie longueur d'onde : λae' = λ/cos(angle(e,e')) Le long des axes on observe les longueurs d'onde apparentes λx = 2π/kx , λy = 2π/ky , λz = 2π/kz.

polarisation psvsh - déplacement P - déplacement SV - déplacement SH
ondePS.m

Ondes sphériques

Les ondes P sphériques émises par une source de compression ponctuelle à symétrie sphérique sont solution d'une équation des ondes 1D pour (rφ(r)).
Le déplacement radial des particules comporte un terme en 1/r, dominant loin de la source, et un terme en 1/r2, dominant près de la source.
Pour une onde harmonique, la limite champ proche/lointain est définie par la valeur de kr, c'est à dire par le rapport r/λ.
Pour une source sphérique où le déplacement U0 est imposé en r = r0<<λ, on obtient l'amplitude de l'onde en identifiant à U0 le terme de déplacement en champ proche en r0.

Δφ(r) = (1/r)∂2(rφ)/∂r2 car Δφ(r) = divgradφ(r) = div(∂φ/∂rer) = grad∂φ/∂r.er+∂φ/∂rdiver = ∂2φ/∂r2+2/r∂φ/∂r (diver = 2/r car diver = div(r/r) = grad(1/r).r+divr/r et divr = ∂x/∂x+∂y/∂y+∂z/∂z = 3)
2(rφ)/∂t2 = VP22(rφ)/∂r2 déplacement des particules champ lointain : kr>>1 champ proche : kr<<1
φ = (A/r)f(±r/VP-t) u = (A/r)er(-f/r±f'/VP) ucl = (±f'A/rVP)er ucp = (-Af/r2)er
φ = (A/r)ei(kr-ωt) u = (A/r)er(ik-1/r)ei(kr-ωt) ucl = (ikA/r)erei(kr-ωt) ucp = (-A/r2)erei(kr-ωt)
A = -U0r02 pour source de rayon r0<<λ

Fronts d'onde et rayons en milieu homogène isotrope

La propagation d'une onde à partir d'un point M quelconque atteint par l'onde à un temps t se fait dans une direction e à la vitesse de propagation V. Un front d'onde correspond à la surface T(x) sur laquelle la perturbation est localisée au temps t. Les rayons représentent les chemins le long duquel l'onde progresse avec la vitesse de propagation propre au milieu. Dans un milieu isotrope homogène, les rayons sont des droites orthogonales aux fronts d'onde. Un front d'onde balaye une direction e' différente de la direction de propagation e avec une vitesse apparente Vae' supérieure à la vitesse de propagation : Vae' = V/cos(angle(e,e')). Les vitesses apparentes dans trois directions mutuellement orthogonales sont liées géométriquement par l'équation iconale. La loi de Snell-Descartes exprime l'égalité des vitesses apparentes des fronts d'onde incidents, réfléchis et transmis sur une interface entre deux milieux homogènes (permettant ainsi la continuité des déplacements et contraintes de part et d'autre de l'interface au passage des ondes).

. équation front onde vitesses apparentes le long des axes équation iconale : Pythagore pour les ondes
plan T(x) = e.x/V gradT = e/V e.e=1
gradT.gradT = (∂T/∂x)2+(∂T/∂y)2+(∂T/∂z)2 = 1/V2
(1/Vax)2+(1/Vay)2+(1/Vaz)2 = 1/V2
sphérique T(x) = r/V gradT = er/V
quelconque T(x) gradT = eM/V
La propagation des fronts d'onde obéit au principe de Huygens : le front d'onde au temps t+dt est l'enveloppe des ondelettes sphériques de rayon Vdt centrées sur le front d'onde au temps t.
La géométrie des rayons découle du principe de Fermat : le temps de trajet est stationnaire par rapport à toute perturbation de la géométrie du rayon.
Ces deux principes sont directement utilisables pour déterminer la propagation d'une onde. Ils illustrent concrètement comment une onde trouve sa direction de propagation : elle les essaie toutes et choisit celles pour lesquelles il y a interférence constructive entre trajets voisins.

Source ponctuelle à distance d d'une interface plane : onde directe, réfléchie, conique

Un front d'onde sphérique incident sur une interface plane engendre une onde sphérique réfléchie et une onde transmise. Si la vitesse augmente au passage de l'interface (V2>V1), seule la partie du front d'onde sphérique incidente sous des angles inférieurs à l'incidence critique θc est transmise. La partie correspondant aux incidences supérieures à l'incidence critique est réfléchie totalement. Le front d'onde transmis, initialement tiré par le front d'onde incident sur l'interface, se propage alors horizontalement sous l'interface à la vitesse V2, supérieure à la vitesse apparente des ondes incidentes et réflechies sur l'interface. Il est en avance par rapport à l'onde incidente et engendre par diffraction un front d'onde conique qui remonte dans le milieu 1 sous l'incidence critique : l'onde conique est le "sillage" dans le milieu 1 de l'onde transmise. Pour une interface non pentée parallèle à la surface d'enregistrement, la vitesse apparente de l'onde conique en surface est la vitesse de propagation de l'onde transmise le long de l'interface. Pour une interface pentée faisant un angle α avec la surface d'enregistrement, la vitesse apparente de l'onde conique en surface est différente suivant que les récepteurs sont situés vers l'aval (Vav = V1/sin(θc+α)) ou vers l'amont (Vam = V1/sin(θc-α)) par rapport à la source ponctuelle.

Hodochrone T(X) des ondes directe, réfléchie, conique pour distance source-récepteurs SR=X et milieux homogènes de vitesses V1 et V2
pendage α onde directe onde réfléchie
réflexion totale pour incidence θ1c
onde conique (n'existe que si V2>V1)
incidence critique θc : sinθc = V1/V2 ;
interface horizontale
direco.m
T = X/V1 T = (4d2+X2)½/V1 distance et temps critique : Xc = 2dtgθc , Tc = 2d/cosθc/V1

vitesses apparentes horizontale et verticale : Vax = V1/sinθc = V2 ; Vaz = V1/cosθc

temps de trajet : T = X/Vax+2d/Vaz = X/V2+2dcosθc/V1

interface pentée T = X/V1 source virtuelle symétrique de la source par rapport à l'interface
xv = 2dsinα ; zv = 2dcosα

T = (4(dcosα)2+(X±2dsinα)2)½/V1
+ pour récepteurs vers l'aval
- pour récepteurs vers l'amont

distance critique (triangle source virtuelle, projeté vertical de la SV sur surface, récepteur) :
Xc = 2d(cosαtg(θc±α)-±sinα) = 2dsinθ/cos(θ±α)

vitesses apparentes suivant et perpendiculairement à l'interface :
V = V1/sinθc = V2 ; Va⊥α = V1/cosθc

distances suivant et perpendiculairement à l'interface :
xα = Xcosα ; d⊥α = 2d±Xsinα

temps de trajet : T = xα/V+d⊥α/Va⊥α
T = Xcosα/V2+(2d±Xsinα)cosθc/V1 = Xsin(θc±α)/V1+2dcosθc/V1

Si la droite source-récepteurs n'est pas parallèle à la direction de plus grande pente, la propagation se fait dans le plan orthogonal à l'interface passant par SR. α est le pendage apparent dans le plan de propagation, tel que sinα = sin(pendage)cos(azimut(RS)-azimut(pendage)).

Pour une interface entre deux milieux élastiques, il faut prendre en compte les ondes converties. Voici les figures dessinées par Cagniard en 1939 dans le cas où VP2>VS2>VP1>VS1.
ondes coniques converties

Pour une source ponctuelle située sur la surface libre d'une couche d'épaisseur d surmontant un demi-espace, les réflexions successives sur l'interface et la surface libre donnent naissance à des réverbérations. Les points de tangence commune entre les réflexions et les ondes coniques multiples sont alignés sur une droite de pente Tc / Xc = V2/V12, intermédiaire entre celle des ondes directes et des ondes coniques. (voir aussi ce sujet d'exam).

Ocean Acoustics Library : propagation des ondes directe, réfléchie, transmise, conique est une animation qui permet de comprendre la physique de l'interaction d'une onde sphérique avec une interface plane.
Heelan P.A., 1953, On the theory of head waves
Mitchell J.F. et R.J. Bolander, 1986, Structural interpretation using refraction velocities from marine seismic surveys montrent un bel enregistrement d'un point de tir en mer :
point de tir

Allen, R., Fratta D., Introduction to applied geophysics, mettent à disposition leur cours sur l'utilisation des ondes coniques en sismique réfraction.

Rayons, distances et temps de propagation pour V(z) - VRMS

Si la vitesse de propagation ne dépend que de la profondeur (stratification horizontale ou gradient vertical de vitesse), la géométrie des rayons est déterminée par la loi de Snell : le paramètre p d'un rayon (inverse de la vitesse apparente horizontale de l'onde) est constant. L'angle d'incidence θ(z) du rayon par rapport à la verticale augmente ou diminue avec l'augmentation ou la diminution de vitesse.
La distance horizontale X(p) parcourue le long d'un rayon et le temps de propagation T(p) s'expriment en fonction du paramètre p en sommant sur les intervalles de profondeur correspondant à l'épaisseur des couches dans un milieu stratifié.
En complément à la représentation de l'hodochrone T(X), il est utile de représenter le plan τ-p, où τ(p) = T(p)-pX(p) représente le temps vertical de propagation dans la stratification pour un rayon de paramètre p.
Dans un milieu à gradient de vitesse, la sommation discrète est remplacée par une intégration sur des intervalles d'épaisseur dz.
Si le gradient de vitesse a est constant (V(z) = V0+az), les rayons sont des arcs de cercle. La construction géométrique des rayons et le calcul de X(p) et T(p) dans ce cas se font simplement car les arcs de cercle correspondant aux différents paramètres sont tous centrés sur la droite V(z) = 0.
Dans ce cas, les fronts d'onde sont aussi des cercles orthogonaux à la famille de cercles des rayons. Ils sont tous centrés sur l'axe Oz à des profondeurs croissantes avec le temps de propagation.
Si on se limite à des angles d'incidence petits, on peut approximer T(X) entre deux profondeurs sous une forme hyperbolique analogue à celle dans un milieu de vitesse constante. La vitesse moyenne du milieu équivalent est VRMS, couramment utilisée dans le traitement des données de sismique réflexion.

loi de Snell : p = sinθ/V = (1/Vax) = cte pour un rayon paramètre du rayon distance horizontale X(p) entre deux points d'un même rayon temps de trajet T(p) entre deux points d'un même rayon
stratification plane horizontale (couches d'épaisseur di et vitesses constantes Vi)
VRMS - VRMS
hodostratif.m - hodostrat.m
p = sinθi/Vi X(p) = ∑ditgθi = p∑diVi/(1-p2Vi2)½ T(p) = ∑di/cosθiVi = ∑di/(1-p2Vi2)½Vi

τ(p) = T(p)-pX(p) = ∑dicosθi/Vi = ∑di(1-p2Vi2)½/Vi

gradient de vitesse V(z) (pour un gradient positif, un rayon atteint une profondeur maximum Zmax(p)) p = sinθ(z)/V(z) = 1/V(Zmax) X(p) = p∫dzV(z)/(1-p2V(z)2)½ T(p) = ∫dz/(1-p2V(z)2)½V(z)

τ(p) = ∫dz(1-p2V(z)2)½/V(z)

cas V(z) = V0+az : les rayons sismiques sont des cercles de rayon R(p)=1/(ap) centrés sur la droite V(z) = 0 sinθ(z) = (V0/a+z)/R = pV(z) X(p) = R[cosθ] = R[(1-p2V(z)2)½] T(p) = ∫Rdθ/V(z(θ)) = (1/a)∫dθ/sinθ = (1/a)[Log(tg(θ/2))] = (1/a)[Log(pV(z)/(1+(1-p2V(z)2)½))]
cas V(z) = V0+az : les fronts d'onde à T = -(1/a)[Log(tg(θ0/2))] sont des cercles orthogonaux aux rayons d'équation x2+(z-Zmax(T))2 = X(T)2
gradient vitesse
hodovz.m
tg(θ0/2) = e-aT X(T) = V0/atgθ0 = (V0/2a)(eaT-e-aT)

X(T) = (V0/a)sh(aT)

Zmax(T) = (V0/a)(1/sinθ0-1) = (V0/2a)(eaT+e-aT-2)

Zmax(T) = (V0/a)(ch(aT)-1)

approximation VRMS pour propagation subverticale
VRMS
vrms.m
θ ≈ 0
p ≈ θi/Vi
X(p) ≈ p∑diVi T(p) ≈ ∑di(1+p2Vi2/2)/Vi = ∑di/Vi+(p2/2)∑diVi = T(X=0)+X2/2∑diVi
T(X)2 ≈ T(0)2+X2/VRMS2
VRMS2 = (∑diVi)/(∑di/Vi) = (∫V(z)dz)/(∫dz/V(z))
VRMS(V0+az) = (a(V0z+az2/2)/Log(1+az/V0))½

En général, V augmente régulièrement avec z. Les rayons issus d'une source ponctuelle, courbés vers la surface, émergent à des distances horizontales croissantes pour des incidences θ0 décroissant de 90° à 0°. Le front d'onde arrive avec une vitesse apparente horizontale croissante en surface.
Si le gradient de vitesse augmente fortement dans un intervalle de profondeur, l'hodochrone présente une triplication
Si le gradient change de signe dans un intervalle de profondeur, il se forme une zone d'ombre : aucun rayon n'émerge dans l'intervalle de X correspondant.
Si la source de trouve dans la zone de faible vitesse, il se forme un guide d'onde pour les rayons quittant la source avec des incidences telles que le point bas (ou haut) du rayon est aussi dans la zone.
Un front d'onde plan incident sur une zone où V augmente avec z est totalement réfléchi si la zone contient le point bas des rayons.

triplication - zone d'ombre - guide d'onde
hodovzt.m - hodovzo.m

Lutter W.J., R. D. Catchings and C. M. Jarchow, 1994, An image of the Columbia Plateau from inversion of high-resolution seismic data montrent l'effet d'une stratification basalte (5 km/s) - sédiment (4 km/s) -socle (5.5 km/s) sur les premières arrivées.
Hornby B., 1993, Tomographic reconstruction of near-borehole slowness using refracted borehole sonic arrivals est un exemple d'utilisation d'ondes coniques et réfractées en diagraphie.
Russel D., Réfraction du son dans l'atmosphère
Garcés, M.A., Hansen, R. A. & Lindquist, K.G., 1998, Traveltimes for infrasonic waves propagating in a stratified atmosphere

Rayons, fronts d'onde, amplitude pour V(x,y,z)

Pour un milieu stratifié comportant des couches séparées par des interfaces pentées, on ne dispose plus d'une normale commune à toutes les interfaces. La loi de Snell s'écrit pour chaque interface en fonction de la direction du rayon incident e et de la normale à l'interface n. Les rayons réfléchis et transmis sont contenus dans le plan (e , n). Soit m le vecteur normal à n dans le plan (e , n), θ1 l'angle d'incidence par rapport à n), er et et les directions des rayons réfléchi et transmis, θ2 l'angle (et , n)

e = -cosθ1n+sinθ1m er = cosθ1n+sinθ1m er-e = 2cosθ1n = -2(e.n)n er = e-2(e.n)n
e/V1 = -cosθ1/V1n+sinθ1/V1m et/V2 = -cosθ2/V2n+sinθ2/V2m et/V2-e/V1 = (cosθ1/V1-cosθ2/V2)n et = V2/V1(e-(e.n+((V1/V2)2-(1-e.n2))½)n)

Voici un exemple pour un modèle comportant deux interfaces planes pentées. L'interface supérieure est la limite sédiments-socle. L'interface inférieure est le réflecteur dans le socle. Les trajets correspondent à des réflexions PP et PS sur le réflecteur pour une source et un récepteur placés à la limite sédiments-socle. La point de réflexion PP est l'intersection du segment source virtuelle-récepteur avec le réflecteur. Le point de réflexion PS est sur l'intersection du plan de propagation PP avec le réflecteur. Il correspond au temps minimum pour des trajets PS diffractés pour des points diffractants PS situés sur ce segment. Les rayons transmis dans les sédiments émergent en des points différents à la surface.

réflexion 3D - réflexion 3D -
refl3d.m

Si la vitesse de propagation varie continuement avec la profondeur et latéralement, les vitesses apparentes horizontales et verticales varient le long d'un même rayon. La variation le long d'un rayon de l'inverse d'une vitesse apparente dans une direction est égale à la composante du gradient de 1/V dans cette direction. Cette relation permet de déterminer la géométrie (tracer) des rayons par intégration en ajustant à chaque pas d'abscisse curviligne s la direction de propagation e. Dans la cas où V(x,z) varie linéairement avec x et z, les rayons sismiques sont des arcs de cercles centrés sur la droite V(x,z) = 0.

s abscisse curviligne d'un rayon x(s)
e = dx/ds
équation différentielle des rayons
e.e = 1 , e.e/∂x = 0 , gradT.e/∂x = 0
∂(dT/ds)/∂x = ∂(e.gradT)/∂x = e.grad∂T/∂x = d(∂T/∂x)/ds
courbure des rayons
de/ds = eN/R
torsion des rayons
b = e∧eN
deN/ds = -e/R-b
gradT = e/V
T = ∫ds/V(x(s))
d(e/V)/ds = d(gradT)/ds = grad(dT/ds) = grad(1/V) 1/R = VeN.grad(1/V) 1/τ = -Rb.d(Vgrad(1/V))/ds
cas V(x,z) = V0+axx+azz = V0+a(xsinα+zcosα)
par une rotation d'angle α, on se ramène au cas V(Z) = V0+aZ
les rayons sismiques sont des cercles centrés sur la droite V(x,z) = 0
d(sinθ/V)/ds = -ax/V2 = -asinα/V2
d(cosθ/V)/ds = -az/V2 = -acosα/V2
1/R = dθ/ds = -axcosθ/V+azsinθ/V = asin(θ-α)/V

En cherchant une solution à l'équation des ondes avec une vitesse V(x) sous forme d'une onde harmonique quasi-plane portée par un front d'onde d'équation T(x) et ayant une amplitude A(x), on obtient des termes dépendant des puissances successives de la pulsation. Le terme en ω2, dominant à haute fréquence, donne l'équation iconale des temps de trajet. Le terme en ω donne l'équation de transport permettant de calculer l'amplitude A(x).

onde quasi plane (solution haute fréquence de l'éq. des ondes) terme en ω2 → éq.iconale terme en ω → éq. transport tube de rayon de section dS(s)
φ = A(x)eiω(T(x)-t) gradT.gradT=1/V2 2gradA.gradT+AΔT = div(A2gradT) = 0 ∫∫∫div(A2gradT)dυ = ∫∫(A2/V)e.dσ = 0
A2dS/V = cte


Langan, R.T., I. Lerche, and R. T. Cutler, 1985, Tracing of rays through heterogeneous media: An accurate and efficient procedure
Bernasconi G. and G. Drufuca, 2001, 3-D traveltimes and amplitudes by gridded rays proposent une méthode simple de tracé de rayons 3D par une méthode de différences finies.
Sava P. and S. Fomel, M. 2001, 3-D traveltime computation using Huygens wavefront tracing montrent comment le principe de Huygens permet de déterminer des fronts d'onde et temps de trajets.

Continuité des déplacements et des contraintes sur une interface plane

On obtient les coefficients de réflexion-transmission des ondes élastiques sur une interface plane entre deux milieux élastiques en écrivant les conditions de continuité des composantes du déplacement et des contraintes qui s'appliquent sur l'interface de part et d'autre de l'interface. Les contraintes qui ne s'appliquent pas sur l'interface n'interviennent pas dans les conditions de continuité. Pour une surface libre, les contraintes sont nulles hors du milieu élastique et doivent donc s'annuler sur la surface libre. Par contre, le déplacement de la surface libre est quelconque (parfois destructeur).
Les calculs sont faits avec des ondes planes harmoniques pour lesquelles
la fréquence et la longueur d'onde apparente horizontale sont identiques pour l'onde incidente et les ondes réfléchies et transmises (loi de Snell). Les valeurs de kz pour chaque onde dépendent de sa nature (P ou S), de son angle d'incidence (les angles sont liés par la loi de Snell), et de son sens de propagation vertical (montante ou descendante) qui fixe le signe de kz.
Voici les différents cas possible suivant l'interface et la nature des ondes se propageant dans le plan vertical (x,z) et incidentes sur une interface plane horizontale :

ω et kx = ωp = ωsin(θ)/VP = ωsin(η)/VS de l'onde incidente sont les mêmes pour toutes les ondes réfléchies et transmises sur une interface plane horizontale
ondes descendantes : kzP = +ωcos(θ)/VP , kzS = +ωcos(η)/VS
ondes montantes :      kzP = -ωcos(θ)/VP , kzS = -ωcos(η)/VS
ondes SHondes P-SV
surface libre σzyIncidente SHzyRéfléchie SH = 0 σzzIncidente P ou SVzzRéfléchie PzzRéfléchie SV = 0
σzxIncidente P ou SVzxRéfléchie PzxRéfléchie SV = 0
interface solide - solide uIncidente SH+uRéfléchie SH = uTransmise SH
σzyIncidente SHzyRéfléchie SH = σzyTransmise SH
uxIncidente P ou SV+uxRéfléchie P+uxRéfléchie SV = uxTransmise P+uxTransmise SV
uzIncidente P ou SV+uzRéfléchie P+uzRéfléchie SV = uzTransmise P+uzTransmise SV
σzzIncidente P ou SVzzRéfléchie PzzRéfléchie SV = σzzTransmise PzzTransmise SV
σzxIncidente P ou SVzxRéfléchie PzxRéfléchie SV = σzxTransmise PzxTransmise SV
interface liquide - solide . uzIncidente P+uzRéfléchie P = uzTransmise P+uzTransmise SV
σzzIncidente PzzRéfléchie P = σzzTransmise PzzTransmise SV
σzxTransmise PzxTransmise SV = 0

Réflexion d'une onde plane harmonique SH sur une interface plane

Une onde plane harmonique SH se propageant dans le plan vertical (x,z), incidente sur une interface plane horizontale entre deux milieux élastiques, provoque sur cette interface une perturbation du déplacement uy et de la contrainte de cisaillement σzy. Cette perturbation balaye l'interface à la vitesse apparente horizontale de l'onde incidente VS1/sinη1, où η1 est l'angle du rayon incident par rapport à la verticale. Elle est équilibrée par les perturbations dues aux ondes SH réfléchie et transmise qui balayent l'interface à la même vitesse apparente horizontale (loi de Snell). Le déplacement uy et la contrainte σzy résultant des ondes SH incidente et réfléchie d'un coté de l'interface sont égaux à uy et σzy dus à l'onde SH transmise de l'autre coté de l'interface. On obtient ainsi les coefficients de réflexion et transmission, définis ici pour le déplacement des particules.
Dans le cas où VS1>VS2 et η1c, l'onde transmise devient évanescente : son amplitude décroit exponentiellement avec la distance à l'interface. La réflexion totale s'accompagne d'un déphasage χ à la réflexion donné par l'argument du coefficient de réflexion complexe. La somme de l'onde incidente et de l'onde réfléchie totalement produit une onde stationnaire verticalement et progressive horizontalement.

loi de Snell :
kx = ωp = ωsinη1/VS1 = ωsinη2/VS2
onde incidente SH descendante onde réfléchie SH montante onde transmise SH descendante
déplacement onde SH
uy(x,z,t) = Uyeikzzei(kxx-ωt)
UyI = 1
kzI = ωcosη1/VS1
UyR = R
kzR = -ωcosη1/VS1
UyT = T
kzT = ωcosη2/VS2
contrainte sur surface horizontale
σzy(x,z,t) = μ∂uy/∂z = ρVS2ikzuy = iωSzyeikzzei(kxx-ωt)
SzyI = ρ1VS1cosη1 SzyR = -ρ1VS1cosη1R SzyT = ρ2VS2cosη2T
continuité sur interface z =0
(uyI+uyR)(z=0) = uyT(z=0)
zyIzyR)(z=0) = σzyT(z=0)
1+R=T
(1-R)ρ1VS1cosη1 = Tρ2VS2cosη2
R = (ρ1VS1cosη12VS2cosη2)/(ρ1VS1cosη12VS2cosη2)
T = 2ρ1VS1cosη1/(ρ1VS1cosη12VS2cosη2)
faible contraste
ρ+δρ , V+δV
δη = tgηδV/V (δp=0)
R ≈ -δ(ρVcosη)/2ρVcosη = -(1/2)(δρ/ρ+(1-tg2η)δV/V) T ≈ 1
réflexion totale
sinηc=VS1/VS2
η1c , pVS2>1
kzR = -kzI = -ω(1-p2VS12)½/VS1
R = (ρ1VS1(1-p2VS12)½-iρ2VS2(p2VS22-1)½)/(ρ1VS1(1-p2VS12)½+iρ2VS2(p2VS22-1)½)
R = eiχ(p) , χ(p) = -2Arctg(ρ2VS2(p2VS22-1)½1VS1(1-p2VS12)½)
(uyI+uyR)(x,z,t) = (eikzIz+ei(kzRz+χ(p)))ei(kxx-ωt) = 2cos(kzIz-χ(p)/2)ei(kxx-ωt+χ(p)/2)
onde transmise évanescente : kzT = iω(p2VS22-1)½/VS2
T = 2ρ1VS1(1-p2VS12)½/(ρ1VS1(1-p2VS12)½+iρ2VS2(p2VS22-1)½)
T = ||T||eiχ(p)/2
uyT(x,z,t) = ||T||e-ωz(p2VS22-1)½/VS2ei(kxx-ωt+χ(p)/2)

Voici les coefficients de réflexion et transmission et le déplacement obtenus pour des angles d'incidence croissant dans le cas de la réflexion partielle et totale :
coefficients réflexion-transmission SH - réflexion totale SH
reflsh.m

Voici une animation dans le cas de la réflexion sur une surface libre (R=1) qui produit des noeuds de vibration (contrainte) pour z = nπ/kz et des ventres pour z = (n+1/2)π/kz.

Ocean Acoustics Library : réflexion totale d'une onde plane sur une interface plane

Réflexion d'une onde plane harmonique P ou SV sur une surface plane libre

Les ondes P et SV se propageant dans le plan vertical (x,z) exercent une compression σzz et un cisaillement σzx sur la surface libre horizontale d'un demi-espace. Elles sont donc couplées à la réflexion : une onde P ou SV incidente se réfléchit sous forme d'une onde P et d'une onde SV. Les contraintes σzz et σzx résultant de l'onde incidente et des deux ondes réfléchies doivent s'annuler sur la surface libre. On obtient ainsi les coefficients de réflexion sur la surface libre, définis ici pour les potentiels de déplacement des particules.
Dans le cas où une onde SV est incidente avec un angle η supérieur à l'incidence critique correspondant au rapport VS/VP, l'onde P réfléchie par la surface libre devient évanescente. L'onde SV est totalement réfléchie avec un déphasage χ' à la réflexion donné par l'argument du coefficient de réflexion complexe.

loi de Snell :
kx = ωp = ωsinθ/VP = ωsinη/VS
onde incidente P montante onde réfléchie P descendante onde réfléchie SV descendante
potentiels de déplacement P et SV
φ(x,z,t) = Φ eikzPzei(kxx-ωt)
ψy(x,z,t) = Ψ eikzSzei(kxx-ωt)
ΦI = 1
kzP = -ωcosθ/VP
ΦR = Rφ
kzP = ωcosθ/VP
ΨR = Rψ
kzS = ωcosη/VS
contraintes exercées par onde P sur une surface horizontale
σzzP = λΔφ+2μ∂2φ/∂z2 = -ρω2(1-2p2VS2)φ = ρω2SzzeikzPzei(kxx-ωt)
σzxP = 2μ∂2φ/∂z∂x = -2ρVS2ωpkzφ = ρω2SzxeikzPzei(kxx-ωt)
Szz = -(1-2p2VS2)
Szx = 2VS2pcosθ/VP
Szz = -(1-2p2VS2)Rφ
Szx = -(2VS2pcosθ/VP)Rφ
.
contraintes exercées par onde SV sur une surface horizontale
σzzS = 2μ∂2ψy/∂z∂x = -2ρVS2ωpkzψ = ρω2SzzeikzSzei(kxx-ωt)
σzxS = μ(-∂2ψy/∂z2+∂2ψy/∂x2) = ρω2(1-2p2VS2)ψ = ρω2SzxeikzSzei(kxx-ωt)
. . Szz = -(2VS2pcosη/VS)Rψ
Szx = (1-2p2VS2)Rψ
onde P incidente, contraintes nulles sur surface libre z = 0
zzPI + σzzPR + σzzSR)(z=0) = 0
zxPI + σzxPR + σzxSR)(z=0) = 0
(1-2p2VS2)(1+Rφ)+(2VS2pcosη/VS)Rψ = 0
(2VS2pcosθ/VP)(1-Rφ)+(1-2p2VS2)Rψ = 0
réflexion P-SV surface libre D = (1-2p2VS2)2+4VS4p2cosθ/VPcosη/VS
Rφ = (-(1-2p2VS2)2+4VS4p2cosθ/VPcosη/VS)/D
Rψ = (-4(1-2p2VS2)VS2pcosθ/VP)/D
onde SV incidente, contraintes nulles sur surface libre z = 0
zzSI + σzzPR + σzzSR)(z=0) = 0
zxSI + σzxPR + σzxSR)(z=0) = 0
(1-2p2VS2)R'φ+(2VS2pcosη/VS)(-1+R'ψ) = 0
-(2VS2pcosθ/VP)R'φ+(1-2p2VS2)(1+R'ψ) = 0
réflexion P-SV surface libre R'φ = (4(1-2p2VS2)VS2pcosη/VS)/D
R'ψ = (-(1-2p2VS2)2+4VS4p2cosθ/VPcosη/VS))/D
réflexion totale de l'onde SV , sinη > VS/VP (η > 36°)
onde P réfléchie évanescente : kzP = iω(p2VP2-1)½/VP
R'ψ = eiχ'(p)
χ'(p) = -2Arctg((4VS4p2(p2VP2-1)½/VPcosη/VS)/(1-2p2VS2)2)

Réflexion d'une onde plane harmonique P ou SV sur une interface plane

Sur une interface plane horizontale entre deux demi-espaces élastiques, les ondes P et SV se propageant dans le plan vertical (x,z) sont reliées par la loi de Snell et par la continuité des composantes du déplacement ux et uz et des contraintes σzz et σzx de part et d'autre de l'interface. Une onde incidente P ou SV génère deux ondes P et SV réfléchies et deux ondes P ou SV transmises. Les quatre équations de continuité s'écrivent sous forme d'une matrice multipliée par le vecteur des coefficients de réflexion-transmission, définis ici pour les potentiels de déplacement des particules.
Pour une interface liquide-solide, il y a continuité de uz et σzz de part et d'autre de l'interface et le cisaillement σzx dans le milieu élastique doit s'annuler à l'interface.
Dés que l'angle d'incidence est supérieur à une valeur d'incidence critique, l'onde pour laquelle la valeur d'incidence est supérieure à 90° devient évanescente. Les coefficients de réflexion-transmission deviennent complexes. Le déphasage χ' à la réflexion totale de l'onde SV est donné par l'argument du coefficient de réflexion.

loi de Snell : kx = ωp = ωsinθ1/VP1 = ωsinη1/VS1 = ωsinθ2/VP2 = ωsinη2/VS2
P incidente descendante
ΦI = 1
kzI = ωcosθ1/VP1
P réfléchie montante
ΦR = Rφ
kzP = -ωcosθ1/VP1
SV réfléchie montante
ΨR = Rψ
kzS = -ωcosη1/VS1
P transmise descendante
ΦT = Tφ
kzP = ωcosθ2/VP2
SV transmise descendante
ΨT = Tψ
kzS = ωcosη2/VS2
Ux = p + pRφ + (cosη1/VS1)Rψ = pTφ -(cosη2/VS2)Tψ
Uz = cosθ1/VP1 -(cosθ1/VP1)Rφ + pRψ = (cosθ2/VP2)Tφ + pTψ
Szz = ρ1(1-2p2VS12) + ρ1(1-2p2VS12)Rφ - (2ρ1VS12pcosη1/VS1)Rψ = ρ2(1-2p2VS22)Tφ + (2ρ2VS22pcosη2/VS2)Tψ
Szx = -2ρ1VS12pcosθ1/VP1 + (2ρ1VS12pcosθ1/VP1)Rφ + ρ1(1-2p2VS12)Rψ = -(2ρ2VS22pcosθ2/VP2)Tφ + ρ2(1-2p2VS22)Tψ
écriture matricielle
cas liquide-solide
AR = B ; R = [Rφ , Rψ , Tφ , Tψ]t ; B = [-p , -cosθ1/VP1 , -ρ1(1-2p2VS12) , 2ρ1VS12pcosθ1/VP1]t
ALRL = BL ; RL = [Rφ , Tφ , Tψ]t ; BL = [-cosθ1/VP1 , -ρ1 , 0]t
A[4,4]
AL[3,3]
p cosη1/VS1 -p cosη2/VS2
-cosθ1/VP1 p -cosθ2/VP2 -p
ρ1(1-2p2VS12) 1) -2ρ1VS12pcosη1/VS1 2(1-2p2VS22) -2ρ2VS22pcosη2/VS2
1VS12pcosθ1/VP1 (0) ρ1(1-2p2VS12) 2VS22pcosθ2/VP2 2(1-2p2VS22)
P incidente descendante ou montante
coefficients réflexion-transmission P-SV
reflpsv.m
Rφ = det([B,A[:,2:4]])/det(A)

tφ = det([b,A[:,2:4]])/det(A)

Rψ = det([A[:,1],B,A[:,3:4]])/det(A)

tψ = det([A[:,1],b,A[:,3:4]])/det(A)

Tφ = det([A[:,1:2],B,A[:,4]])/det(A)

rφ = det([A[:,1:2],b,A[:,4]])/det(A)

Tψ = det([A[:,1:3],B])/det(A)

rψ = det([A[:,1:3],b])/det(A)

P incidente descendante
D = det(A) = P-Q
Rφ = (-P-Q)/D
Rψ =2S/D
AR = B ; R = [Rφ , Rψ , Tφ , Tψ]t ;
B = [-p , -cosθ1/VP1 , -ρ1(1-2p2VS12) , 2ρ1VS12pcosθ1/VP1]t = [-a11 , a21 , -a31 , a41]t ;
D = (a11A11+a31A31)-(a21A21+a41A41) ;
DRφ = -(a11A11+a31A31)-(a21A21+a41A41)
DRψ/2 = a11a21(a33a44-a43a34) +a11a41(a23a34-a33a24) +a31a41(a13a24-a23a14)
DTφ/2 = -(a11a21(a32a44-a42a34) +a11a41(a22a34-a32a24) +a31a41(a12a24-a22a14))
DTψ/2 = (a11a21(a32a43-a42a33) +a11a41(a22a33-a32a23) +a31a41(a12a23-a22a13)
P incidente montante
kzI = -ωcosθ2/VP2
Ar = b ; r = [tφ,tψ,rφ,rψ]t
b = [p , -cosθ2/VP2 , ρ2(1-2p2VS22) , 2ρ2VS22pcosθ2/VP2]t = [-a13 , a23 , -a33 , a43]t
D = (a13A13+a33A33)-(a23A23+a43A43) ;
Drφ = -(a13A13+a33A33)-(a23A23+a43A43)
Drψ/2 = a13a23(a31a42-a41a32) +a13a43(a21a32-a31a22) +a33a43(a11a22-a21a12)
Dtφ/2 = -(a13a23(a34a42-a44a32) +a13a43(a24a32-a34a22) +a33a43(a14a22-a24a12))
Dtψ/2 = a13a23(a34a41-a44a31) +a13a43(a24a31-a34a21) +a33a43(a14a21-a24a11)
cas liquide-solide
coefficients réflexion-transmission liquide-solide
reflpsl.m
ALRL = BL ; RL = [RφL , TφL , TψL]t ; BL = [-cosθ1/VP1 , -ρ1 , 0]t
ALrL = bL ; rL = [tφL , rφL , rψL]t ; bL = [-cosθ2/VP2 , ρ2(1-2p2VS22) , 2ρ2VS22pcosθ2/VP2]t
DL = -ρ22cosθ1/VP1((1-2p2VS22)2+4VS24p2cosθ2/VP2cosη2/VS2) -ρ1ρ2cosθ2/VP2
DLRφL = -ρ22cosθ1/VP1((1-2p2VS22)2+4VS24p2cosθ2/VP2cosη2/VS2) 1ρ2cosθ2/VP2
DLTφL = -2ρ1ρ2cosθ1/VP1(1-2p2VS22) ; DLtφL = -2ρ22cosθ2/VP2(1-2p2VS22)
DLTψL = -4ρ1ρ2cosθ1/VP1VS22pcosθ2/VP2
SV incidente descendante ou montante
coefficients réflexion-transmission SV-P
reflpsv.m
R'φ = det([B',A[:,2:4]])/D

t'φ = det([b',A[:,2:4]])/D

R'ψ = det([A[:,1],B',A[:,3:4]])/D

t'ψ = det([A[:,1],b',A[:,3:4]])/D

T'φ = det([A[:,1:2],B',A[:,4]])/D

r'φ = det([A[:,1:2],b',A[:,4]])/D

T'ψ = det([A[:,1:3],B'])/D

r'ψ = det([A[:,1:3],b'])/D

SV incidente descendante
ΨI = 1 ; kzI = ωcosη1/VS1
D = det(A) = -P'+Q'
R'ψ = (-P'-Q')/D
AR'=B' ; R' = [R'φ , R'ψ , T'φ , T'ψ]t ;
B' = [cosη1/VS1 , -p , -2ρ1VS12pcosη1/VS1 , -ρ1(1-2p2VS12)]t = [a12 , -a22 , a32 , -a42]t ;
D = -(a12A12+a32A32)+(a22A22+a42A42) ;
DR'ψ = -(a12A12+a32A32)-(a22A22+a42A42)
SV incidente montante
kzI = -ωcosη2/VS2
Ar'=b' ; r' = [t'φ , t'ψ , r'φ , r'ψ]t
b' = [cosη2/VS2 , p , -2ρ2VS22pcosη2/VS2 , ρ2(1-2p2VS22)]t = [a14 , -a24 , a34 , -a44]t ;
D = -(a14A14+a34A34)+(a24A24+a44A44)
Dr'ψ = -(a14A14+a34A34)-(a24A24+a44A44)
réflexion totale de l'onde SV descendante
p = sinη1/VS1
1/VS1 > p > 1/VP1
p > 1/VS2 > 1/VP2
R'ψ = eiχ'(p)
termes de A imaginaires : a14 = i(p2VS22-1)½/VS2 ; a21 = -i(p2VP12-1)½/VP1 ; a23 = -i(p2VP22-1)½/VP2 ;
a34 = -2iρ2VS22p(p2VS22-1)½/VS2 ; a41 = 2iρ1VS12p(p2VP12-1)½/VP1 ; a43 = 2iρ2VS22p(p2VP22-1)½/VP2 ;
⇒ A12, A32 et P' sont imaginaires purs, A12, A32 et Q' sont réels
R'ψ = (-P'-Q')/(-P'+Q') = eiχ'(p)
χ'(p) = 2Arctg(P'/iQ') = 2Arctg((1/i)(a12A12+a32A32)/(a22A22+a42A42))

Voici les coefficients de réflexion-transmission P-SV et SV-P calculés pour des angles d'incidence 0°, 15° et 30° dans l'eau pour un modèle de vitesses-densité stratifié type Mer du Nord :
Modèle vitesses-densité Mer du Nord - Réflexion-transmission Mer du Nord

Approximations des coefficients de réflexion P-P et P-SV pour faibles contrastes

Lorsque les paramètres élastiques et la densité varient peu au passage d'une interface plane, les coefficients de réflexion s'expriment linéairement en fonction des variations relatives de ces paramètres. Ces expressions approximées sont utilisées dans l'analyse des amplitudes en sismique réflexion.

faibles contrastes ρ2 = ρ+δρ ; VP2 = VP+δVP ; θ2 = θ+δθ ; δθ = tgθδVP/VP ; VS2 = VS+δVS ; η2 = η + δη ; δη = tgηδVS/VS
a11 = -a13 = a22 = -a24 = p ; a12 = cosη/VS ; a21 = -cosθ/VP ; a31 = a42 = ρ(1-2p2VS2) ; a32 = -2ρVS2pcosη/VS ; a41 = 2ρVS2pcosθ/VP
a14 = (a12+δa12) ; a23 = (a21+δa21) ; a33 = a44 = -(a31+δa31) ; a34 = (a32+δa32) ; a43 = (a41+δa41)
D = det(A) D = (p2-a21a12)(a312-a41a32) - (2pa21)(-2a31a32) + (-p2-a21a12)(a32a41+a312) + (a12a21+p2)(-a312-a41a32) - (-2pa12)(2a31a41) + (p2-a21a12)(a312-a41a32)
D = 2(p2-a21a12)(a312-a41a32) -2(p2+a21a12)(a32a41+a312) +4pa31(a21a32+a12a41) = -4p2a41a32 -4a21a12a312 +4pa31(a21a32+a12a41)

D = 4cosθ/VPcosη/VS [p2(2ρVS2p)2 +(ρ(1-2p2VS2))2 +4ρ2p2VS2(1-2p2VS2)] = 4ρ2cosθ/VPcosη/VS

Rφ DRφ = (-p2-a21a12)((a31+δa31)2-(a41+δa41)(a32+δa32)) - (-p(a21+δa21)+pa21)(-a32(a31+δa31)-a31(a32+δa32)) + (p2-a21(a12+δa12))(a32(a41+δa41)+a31(a31+δa31)) + (a12(a21+δa21)+p2)(a31(a31+δa31)-a41(a32+δa32)) - (-pa12-p(a12+δa12))(-a31(a41+δa41)+a41(a31+δa31)) + (p2-(a21+δa21)(a12+δa12))(-a312-a41a32))
DRφ = (-p2-a21a12)(2a31δa31-a41δa32-a32δa41) - 2pa32a31δa21 + (p2-a21a12)(a32δa41+a31δa31)-a21δa12(a32a41+a312) + (p2+a12a21)(a31δa31-a41δa32)+a12δa21(a312-a41a32) + 2pa12(-a31δa41+a41δa31) + (a21δa12+a12δa21)(a312+a41a32)
DRφ = 2(p(pa32-a12a31)δa41 + a12(pa41-a21a31)δa31 + a31(a31a12-pa32)δa21)

DRφ = 2ρ(-2p2cosη/VSδ(ρVS2cosθ/VP) +cosη/VScosθ/VPδ(ρ(1-2p2VS2) + ρ(1-2p2VS2)cosη/VSδ(-cosθ/VP))
DRφ = 2ρcosη/VS(cosθ/VP((1-4p2VS2)δρ - 4ρp2δ(VS2)) - ρδ(cosθ/VP))

Rφ = (1/2)((1-4p2VS2)δρ/ρ - 8p2VS2δVS/VS - δ(cosθ/VP)/(cosθ/VP))
Rφ = (1/2)((1-4p2VS2)δρ/ρ - 8p2VS2δVS/VS + (1/cosθ)2δVP/VP)
Rφ = (1/2)((δρ/ρ+δVP/VP) + (δVP/VP -4VS2/VP2(δρ/ρ+2δVS/VS)sin2θ + (tg2θ-sin2θ)δVP/VP) ≈ R0+Gsin2θ

Rψ DRψ = (2pa21)((a31+δa31)2-(a41+δa41)(a32+δa32)) - p((a21+δa21)+a21)(a31(a31+δa31)-a41(a32+δa32)) + (-p2-a21(a12+δa12))(-a31(a41+δa41)+a41(a31+δa31)) + p(-(a21+δa21)+a21)(-a31(a31+δa31)-a41(a32+δa32)) - (p2-a21(a12+δa12))(a31(a41+δa41)+a41(a31+δa31)) + (p2-(a21+δa21)(a12+δa12))(2a31a41)
DRψ = 2pa21(2a31δa31-a41δa32-a32δa41) - 2pa21(a31δa31-a41δa32)-pδa21(a312-a41a32) + (-p2-a21a12)(-a31δa41)+a41δa31)) + (-pδa21)(-a312-a41a32) - (p2-a21a12)(a31δa41+a41δa31) +a21δa122a31a41 + (-2a31a41(a21δa12+a12δa21))
DRψ = 2(p(a21a31-pa41)δa31 + a21(-pa32+a12a31)δa41 + a41(pa32-a31a12)δa21)

DRψ = -2ρpcosθ/VP(δ(ρ(1-2p2VS2)) + cosη/VSδ(2ρVS2cosθ/VP) + 2ρVS2cosη/VSδ(-cosθ/VP))
DRψ = -2ρpcosθ/VP ((1-2p2VS2+2VS2cosθ/VPcosη/VS)δρ + 2ρ(-p2+cosη/VScosθ/VP)δ(VS2)

Rψ = -(1/2)pVS/cosη ((1-2p2VS2+2VS2cosθ/VPcosη/VS)δρ/ρ + 4VS2(-p2+cosη/VScosθ/VP)δVS/VS)
Rψ ≈ -(1/2)sinη ((1+2VS/VP)δρ/ρ + 4VS/VPδVS/VS) pour incidence faible

Onde de Rayleigh guidée par une surface libre et de Stoneley guidée par une interface

L' onde de Rayleigh correspond à la propagation le long de la surface libre d'un demi-espace (z ≥ 0 ) de deux ondes P et SV évanescentes de même vitesse de propagation horizontale VR < VS. L'annulation des contraintes sur la surface libre fixe VR. VR est fonction du rapport VP/VS et est indépendante de la fréquence (VR ≈ 0.9VS). Le déplacement des particules en surface est elliptique rétrograde car il existe un déphasage de π/2 entre les composantes ux et uz.

1/VR = p = kx/ω > 1/VS > 1/VP onde P évanescente
kzP = iω(p2VP2-1)½/VP
onde S évanescente
kzS = iω(p2VS2-1)½/VS
σzz = σzx = 0 sur surface libre z = 0 (1-2p2VS2)Φ+i2pVS(p2VS2-1)½Ψ = 0
-i2VS2p(p2VP2-1)½/VPΦ+(1-2p2VS2)Ψ = 0
équation de Rayleigh déterminant = 0
(1-2p2VS2)2-4p2VS4(p2VS2-1)½/VS(p2VP2-1)½/VP = 0
(2-VR2/VS2)2-4(1-VR2/VP2)½(1-VR2/VS2)½ = 0

(ux , uz) = (Ux(z) , Uz(z))eiω(x/VR-t)
Ux = iωpΦe-ωz(p2VP2-1)½/VP+ω(p2VS2-1)½/VSΨe-ωz(p2VS2-1)½/VS
Uz = -ω(p2VP2-1)½/VPΦe-ωz(p2VP2-1)½/VP+iωpΨe-ωz(p2VS2-1)½/VS

Ux = eiπ/2ωpΦ(e-ωz(p2VP2-1)½/VP+(1/(2p2VS2)-1)e-ωz(p2VS2-1)½/VS)
Uz = eω(p2VP2-1)½/VPΦ(e-ωz(p2VP2-1)½/VP+(1/(2p2VS2)-1)-1e-ωz(p2VS2-1)½/VS)

De même, une interface (z = 0) entre deux demi-espaces élastiques peut guider une onde de Stoneley (seulement pour certaines valeurs des rapports de vitesses et densités). Elle est constituée de deux ondes de Rayleigh couplées (deux ondes P1, SV1 évanescentes pour z ≤ 0 couplées à deux ondes P2 , SV2 évanescentes pour z ≥ 0) se propageant à la même vitesse horizontale VST < VS1. VST est déterminée par l'annulation du déterminant de la matrice A obtenue lors du calcul des coefficients de réflexion-transmission pour les ondes P-SV. Il n'existe pas toujours une solution, sauf pour une interface liquide-solide. Dans ce cas, lorque la vitesse des ondes S dans le solide est inférieure à la vitesse des ondes P dans le liquide ("slow formation"), la vitesse de l'onde de Stoneley est proche de celle de VS dans le solide.

eq. de Stoneley
interf. liq-sol.
équation Stoneley
DL = -ρ22cosθ1/VP1((1-2p2VS22)2+4VS24p2cosθ2/VP2cosη2/VS2) -ρ1ρ2cosθ2/VP2 = 0

(1-2p2VS22)2-4VS24p2(p2VP22-1)½/VP2(p2VS22-1)½/VS2 + (ρ12)((p2VP22-1)½/VP2)/((p2VP12-1)½/VP1) = 0

Biot, M.A., 1952, Propagation of elastic waves in a cylindrical bore containing a fluid calcule les vitesses de propagation des ondes de Rayleigh et Stoneley dans le cas d'un forage.
Stephen, R.A., F. Cardo-Casas, and C. H. Cheng , 1985, Finite-difference synthetic acoustic logs montrent la propagation d'une onde de Stoneley dans un puits pour une diagraphie acoustique.

Ondes de Love et de Rayleigh guidées dans une couche d'épaisseur H

L'onde de Love est une onde guidée formée par interférence constructive d'ondes SH réfléchies totalement dans une couche de vitesse VS1 surmontant un demi-espace de vitesse VS2>VS1. La réflexion sur la surface libre se fait sans déphasage, celle sur l'interface avec un déphasage χ , obtenu lors du calcul du coefficient de réflexion d'une onde SH pour un angle d'incidence dans la couche supérieur à ηc. La vitesse de propagation VL de l'onde de Love est la vitesse apparente horizontale des ondes SH dans la couche, telle que VS1≤VL≤VS2. Elle dépend de la fréquence (dispersion). Les différents modes propres s'établissant verticalement dans la couche (mode fondamental et harmoniques) correspondent à une surface libre placée sur les noeuds (contrainte) ou les ventres (déplacement) de la vibration stationnaire imposée par la réflexion totale sur l'interface. Ils correspondent à des incidences différentes des ondes SH rebondissant dans la couche, donc à des courbes de dispersion distinctes.

réflexion totale SH onde stationnaire en z
kz pour modes propres
onde progressive en x , kx = (ω2/VS2-kz2)½
vitesse de phase V(ω) = Vax = ω/kx
ω ≥ fréquence de coupure ωc
kx réel pour onde progressive en x
modes propres 1D
do, do, sol, do, mi, sol...
eikz+e-ikz = 2cos(kz)
Hkn = nπ , nλn/2 = H , n = 1,2,3...
plaque à bords libres
Rsup = Rinf = 1
eikzz+e-ikzz = 2cos(kzz)
Hkz n = nπ , nλz n/2 = H , n = 1,2,3...
kx n = (ω2/VS2-(nπ/H)2)½
Vn(ω) = ω/kx n
ωc n = nπVS/H
couche (-H ≤ z ≤ 0)
sur 1/2 espace (z≥ 0)
Rsurf = 1 ; Rinterf = eiχ(p)
eik1zz+eiχ(p)e-ik1zz = 2eiχ(p)/2cos(k1zz-χ(p)/2)
Hk1z n+χ(p)/2 = nπ , n = 0,1,2,3...
kx n = (ω2/VS12-(nπ-χ(p)/2)2/H2)½
VL n(ω) = ω/kx n = 1/pn
kx nc = ωc n/VS2 ; χ(1/VS2) = 0
ωc n = nπ/(H(1/VS12-1/VS22)½)
relation de dispersion Love
displov.m
k1z = (ω/VS1)(1-p2VS12)½
χ(p) = -2Arctg(ρ2VS2(p2VS22-1)½1VS1(1-p2VS12)½) = -2Arctg(I2/I1)

relation de dispersion de l'onde de Love :
ω(pn) = (-χ(pn)/2+nπ)VS1/(H(1-pn2VS12)½) = (Arctg(I2/I1)+nπ)VS1/(H(1-pn2VS12)½)

De même, une onde SV totalement réfléchie dans une couche telle que VS1<VS2<VP1<VP2 engendre une onde guidée de Rayleigh se propageant à la vitesse VR telle que VS1≤VR≤VS2. La réflexion totale sur la surface libre produit un déphasage χ'1(p) obtenu lors du calcul des coefficients de réflexion P-SV sur une surface libre. Celle à l'interface avec le demi-espace produit un déphasage χ'2(p) obtenu lors du calcul des coefficients à l'interface. On établit la relation de dispersion en tenant compte de la somme χ'12(p) de ces deux déphasages dans l'écriture du mode propre vertical dans la couche. On obtient ainsi les courbes de dispersion des harmoniques (n≥1 , VS1≤VRn≤VS2, VRn=VS2 à la fréquence de coupure fcn) de l'onde de Rayleigh guidée par la couche.

conditions aux limites modes propres en z onde progressive en x fréquence de coupure
Rsurf = eiχ'1(p)
Rinterf = eiχ'2(p)
χ'12(p) = χ'1(p)+χ'2(p)
(n-χ'12(p)/2π)λ1z n/2 = H
Hk1z n+χ'12(p)/2 = nπ
kx n = (ω2/VS12-(nπ-χ'12(p)/2)2/H2)½
VR n(ω) = ω/kx n = 1/pn
ωc n = (nπ-χ'12(1/VS2))/2)/(H(1/VS12-1/VS22)½)
rel. dispersion Rayleigh
harmoniques n≥1
VS1≤VR n≤VS2
χ'1(p) = -2Arctg((4VS14p2(p2VP12-1)½/VP1(1-p2VS12)½/VS1)/(1-2p2VS12)2) = -2Arctg(S)
χ'2(p) = -2Arctg(Im(D)/Re(D)) = -2Arctg(I)
tg(Hω(1-p2VS12)½/VS1) = -tg((χ'1(p)+χ'2(p))/2) = (I+S)/(1-IS)
ω(pn) = (-(χ'1(pn)+χ'2(pn))/2+nπ)VS1/(H(1-pn2VS12)½) = (Arctg((I+S)/(1-IS))+nπ)VS1/(H(1-pn2VS12)½)

La courbe de dispersion du mode fondamental de l'onde de Rayleigh est obtenue en considérant quatre ondes P et SV montantes et descendantes dans la couche et deux ondes P et SV évanescentes dans le demi-espace. La vitesse de propagation VR0 correspond à l'annulation du déterminant de la matrice des conditions aux limites sur la surface libre en z = -H et sur l'interface en z = 0. Elle est telle que VRS1≤VR0≤VRS2 où VRS1 est la vitesse de l'onde de Rayleigh guidée à la surface libre d'un demi-espace de vitesse VS1 et VRS2 celle pour un demi-espace de vitesse VS2. A basse fréquence, VR0 = VRS2 (l'onde de longueur d'onde >> H ne voit pas la couche).

cond. limites descendante P1 , kzP1 = ω(1-p2VP12)½/VP1
ZP1 = e-ikzP1H
descendante S1 , kzS1 = ω(1-p2VS12)½/VS1
ZS1 = e-ikzS1H
montante P1 , -kzP1 montante S1 , -kzS1 évanescente P2 , kzP2 = iω(p2VP22-1)½/VP2 évanescente S2 , kzS2 = iω(p2VS22-1)½/VS2
z = -H σzz(kzP12)ZP1 -2μ1kxkzS1ZS1 σzz(kzP12)/ZP1 1kxkzS1/ZS1 0 0
-2μ1kxkzP1ZP1 σzx(kzS12)ZS1 1kxkzP1/ZP1 σzx(kzS12)/ZS1 0 0
z = 0 kx -kzS1 kx kzS1 -kx kzS2
kzP1 kx -kzP1 kx -kzP2 -kx
σzz(kzP12) -2μ1kxkzS1 σzz(kzP12) 1kxkzS1 zz(kzP22) 2kxkzS2
-2μ1kxkzP1 σzx(kzS12) 1kxkzP1 σzx(kzS12) 2kxkzP2 zx(kzS22)

La partie haute fréquence de la courbe de dispersion du mode fondamental (VRS1≤VR0≤VS1) correspond à l'onde de Rayleigh guidée par la surface libre de la couche couplée à une onde de Stoneley guidée par l'interface (ou deux ondes P et SV évanescentes dans la couche depuis la surface couplées à quatre ondes évanescentes depuis l'interface). A haute fréquence, VR0 = VRS1 (l'onde de longueur d'onde << H ne voit que la couche).

cond. lim. ondes guidées par surf.libre en z = -H et prolongées en z = 0 par
ZP1 = e-ω(p2VP12-1)½/VP1H (0) ; ZS1 = e-ω(p2VS12-1)½/VS1H (0)
ondes guidées par interf. en z = 0 et prolongées en z = -H par ZP1 , ZS1
z = -H Rayleigh surf. libre (1-2p2VS12)ZP1 (0) -i2pVS1(p2VS12-1)½ZS1 (0) 0 0
+i2VS12p(p2VP12-1)½/VP1ZP1 (0) +(1-2p2VS12)ZS1 (0) 0 0
z = 0 pZP1 (0) -i(p2VS12-1)½/VS1ZS1 (0) Stoneley à l'interface
i(p2VP12-1)½/VP1ZP1 (0) pZS1 (0)
ρ1(1-2p2VS12)ZP1 (0) i2ρ1VS12p(p2VS12-1)½/VS1ZS1 (0)
-i2ρ1VS12p(p2VP12-1)½/VP1ZP1 (0) ρ1(1-2p2VS12)ZS1 (0)

Russel D., Guide d'onde

Vitesse de groupe

superposition d'ondes harmoniques de vitesses de phase V(ω)=ω/k interférence constructive en (x,t) ⇔ phase stationnaire pour ondes harmoniques (la variation de U est négligeable) ⇒ vitesse de groupe Vg
u(x,t) = ∫U(ω)eiω(x/V(ω)-t)dω = ∫U(ω)ei(k(ω)x-ωt)dω = ∫U(k)ei(kx-ω(k)t)dk d(k(ω)x-ωt)/dω = (dk/dω)x-t = 0 ⇒ Vg = dω/dk = V+kdV/dk = 1/(dk/dω) = 1/(p+ωdp/dω)
vitesse groupe onde de Love - groupe onde de Love
displov.m

relation de dispersion de l'onde de Love : tg(Hω(pn)(1-pn2VS12)½/VS1) = (ρ2VS2(p2VS22-1)½)/(ρ1VS1(1-p2VS12)½)

dérivation de la relation de dispersion :
(H/VS1)((1-pn2VS12)½-ωpnVS12(1-pn2VS12)dpn/dω)/cos2(Hω(1-pn2VS12)½/VS1)) = (ρ2VS2)/(ρ1VS1)(1-pn2VS12)(pnVS22(pn2VS22-1)+pnVS12(pn2VS22-1)½/(1-pn2VS12))dpn/dω

dpn/dω = (H/VS1)(1-pn2VS12)/
(HωpnVS1+(ρ2VS2)/(ρ1VS1)cos2(Hω(1-pn2VS12)½/VS1)(pnVS22(pn2VS22-1)+pnVS12(pn2VS22-1)½/(1-pn2VS12)))

Kelly K.R. , 1983, Numerical study of Love wave propagation construit des sismogrammes synthétiques par superposition d'ondes de Love dans une couche et illustre l'effet des variations latérales de vitesses.
Roth M. and K. Holliger, 1999, Inversion of source-generated noise in high-resolution seismic data montrent l'effet de la dispersion des ondes guidées dans une couche de faible épaisseur.

Onde de Love guidée en milieu stratifié

Les relations de dispersion s'obtiennent en empilant les couches d'épaisseur hi au dessus de l'interface avec le demi-espace z > 0 où se produit la réflexion totale avec un déphasage χ. Dans chaque couche les ondes harmoniques montante et descendante sont prolongées vers le haut par un déphasage vertical ±kizhi. Si la face supérieure de la couche est la surface libre, la condition de réflexion totale donne la relation de dispersion de l'onde de Love. Si la face supérieure de la couche est une interface avec une couche au dessus, les conditions de continuité sur le déplacement uy et la contrainte σyz donnent les amplitudes des ondes montantes et descendantes dans la couche supérieure (le tableau se lit de bas en haut).

z
Ii = ρiVSi2kiz
onde descendante ekizz
kiz = ω(1/VSi2-p2)½
onde montante e-ikizz continuité uy et σyz à l'interface
réflexion totale si surface libre (montante = descendante) : relation de dispersion
-(h3+h2+h1)
couche 3
D3e-ik3zh3 = D3(a3-ib3) M3eik3zh3 = M3(a3+ib3) A3b3+a3B3 = 0
tg(k3zh3)+(I2/I3)tg(k2zh2)+((I1/I3)-(I1/I2)tg(k3zh3)tg(k2zh2))tg(k1zh1+χ/2) = 0
-(h2+h1)
couche 3
D3 = A3-iB3 M3 = A3+iB3 A3 = A2a2-B2b2 = a1a2-(I1/I2)b1b2
B3 = (I2/I3)(A2b2+B2a2) = (I2/I3)a1b2+(I1/I3)b1a2
-(h2+h1)
couche 2
D2e-ik2zh2 = D2(a2-ib2) M2eik2zh2 = M2(a2+ib2) A2b2+a2B2 = 0
tg(k2zh2)+(I1/I2)tg(k1zh1+χ/2) = 0
-h1
couche 2
D2 = A2-iB2 M2 = A2+iB2 A2 = a1
B2 = (I1/I2)b1
-h1
couche 1
D1e-ik1zh1 = a1-ib1 M1eik1zh1 = a1+ib1 b1 = 0
sin(k1zh1+χ/2) = 0
0
couche 1
D1 = e-iχ/2 M1 = eiχ/2 χ(p) = -2Arctg(ρVS(p2VS2-1)½1VS1(1-p2VS12)½) = -2Arctg(I/I1)
> 0
1/2 espace
onde évanescente e-kzz
kz = |ω|(p2-1/V2)½

Ondes cylindriques

Pour une symétrie cylindrique, les solutions harmoniques de l'équation des ondes font intervenir les fonctions de Bessel et Hankel d'ordre n = 0,1,2... où n correspond à la variation azimutale einθ.
Pour krr>>1, elles sont similaires à des sinusoïdes d'amplitude décroissant comme (1/r)½ :

fonctions Bessel

eq. de Besselsolutiondérivation
r2d2f/dr2+rdf/dr+ (kr2r2-n2)f(r) = 0 f(r) = Hn(krr) = Jn(krr)+iYn(krr) ≈ (2/πkrr)½ei(krr-π/4-nπ/2)

Hn(ikrr) ≈ i-(n+1)(2/πkrr)½e-krr

dHn(krr)/dr = krHn-1(krr)-(n/r)Hn(krr) = (n/r)Hn(krr)-krHn+1(krr)

En l'absence de variation azimutale, les potentiels de déplacement des ondes P et S sont :

Δφ(r,z) = ∂2φ/∂r2+(1/r)∂φ/∂r+∂2φ/∂z2 = (1/VP2)∂2φ/∂t2 φ(r,z) = f(r)ei(kzz-ωt) kp2 = ω2/VP2-kz2 f(r) = AH0(kpr)
Δψθ(r,z) = (1/VS2)∂2ψθ/∂t2 ψθ(r,z) = g(r)ei(kzz-ωt) ks2 = ω2/VS2-kz2 g(r) = BH1(ksr)

Dans le cas où il y a une variation angulaire avec θ, on a :

Δφ(r,θ,z) = Δφ(r,z) + (1/r)22φ/∂θ2 = (1/VP2)∂2φ/∂t2 φ(r,θ,z) =f(r)cos(nθ)ei(kzz-ωt) kp2 = ω2/VP2-kz2 f(r) = AHn(kpr)
Δψr(r,θ,z) = (1/VS2)∂2ψr/∂t2 ψr(r,z) = g(r)sin(nθ)ei(kzz-ωt) ks2 = ω2/VS2-kz2 g(r) = BHn+1(ksr)
Δψθ(r,θ,z) = (1/VS2)∂2ψθ/∂t2 ψθ(r,z) = -g(r)cos(nθ)ei(kzz-ωt)
Δψz(r,θ,z) = (1/VS2)∂2ψz/∂t2 ψz(r,z) = h(r)sin(nθ)ei(kzz-ωt) h(r) = CHn(ksr)

Déplacements et déformations sur la paroi d'un forage :

ur = cos(nθ) uθ = sin(nθ) uz = cos(nθ) εrr = cos(nθ) ε = sin(nθ) εrz = cos(nθ)
Aei(kzz-ωt) n/rHn(kpr)-kpHn+1(kpr) -n/rHn(kpr) ikzHn(kpr) (n(n-1)/r2-kp2)Hn(kpr)+kp/rHn+1(kpr) -n(n-0.5)/r2Hn(kpr)+nkp/rHn+1(kpr) ikz(n/rHn(kpr)-kpHn+1(kpr))
Bei(kzz-ωt) ikzHn+1(ksr) ikzHn+1(ksr) -ksHn(ksr) ikz(ksHn(ksr)-(n+1)/rHn+1(ksr)) ikz(ks/2Hn(ksr)-(n+1)/rHn+1(ksr)) (-nks/rHn(ksr)+(ks2-kz2)Hn+1(ksr))/2
Cei(kzz-ωt) n/rHn(ksr) -n/rHn(ksr)+ksHn+1(ksr) 0 n(n-1)/r2Hn(ksr)-nks/rHn+1(ksr) (-n(n-1)/r2+ks2/2)Hn(ksr)-ks/rHn+1(ksr) (ikzn/2r)Hn(ksr)

Les solutions P et S en coordonnées cylindriques sont développées en détail dans : Sinha, B.K., E. Simsek, and S. Asvadurov, 2009 Influence of a pipe tool on borehole modes

Ondes P et S dues à une force ponctuelle

En l'absence de forces de volume, les 3 équations d'équilibre élastiques écrites en termes des composantes du déplacement des particules s'écrivent :
ρ∂2u/∂t2 = (λ+μ)graddivuΔuΔu = (Δux, Δuy, Δuz) = grad(divu)-rot(rot(u))

On obtient l'équation des ondes P, ∂2φ/∂t2 = VP2Δφ, à partir de ces 3 équations en posant u = gradφ. En effet Δgradφ = gradΔφ car rot(gradφ) = 0. Donc ρgrad2φ/∂t2 = (λ+2μ)grad(Δφ).
On obtient l'équation des ondes S, ∂2ψ/∂t2 = VS2Δψ, en posant u = rotψ avec divψ = 0. En effet Δrotψ = -rotrotrotψ = rotΔψ car divrotψ = 0 et divψ = 0. Donc ρrot2ψ/∂t2 = μrotΔψ.

Pour une source ponctuelle due à une force de module F(t) orientée dans une direction e, la force par unité de volume s'écrit :
f = F(t)eδ(x) = -F(t)/(4π)eΔ(1/r) = -F(t)/(4π)Δ(e/r) = -(F(t)/(4π))(graddiv(e/r)-rotrot(e/r)).
En effet, ∫δ(x)dx = 1 et ∫Δ(1/r)dx = ∫grad(1/r).dS = (-1/r2)er.4πr2er = -4π impliquent δ(x) = -1/(4π)Δ(1/r). Comme e est constant, eΔ(1/r) = Δ(e/r).
On considère le cas harmonique du à une force oscillant sinusoïdalement avec une pulsation ω : F(t) = Fe-iωt.

θ étant l'angle entre la direction e d'application de la force et la direction radiale er de propagation de l'onde depuis le point d'application de la force, on a e = cosθer-sinθeθ.
On utilise l'expression des
opérateurs différentiels en coordonnées sphériques (r,&theta,ξ) : gradφ(r,θ) = &partφ/∂rer+1/r&partφ/∂θeθ et rot(ψ(r,θ)eξ) = (1/r)(1/sinθ∂(sinθψ)/∂θer-∂(rψ)/∂reθ).

Les équations d'équilibre avec la partie du terme source générant les ondes P s'écrivent ρgrad2φ/∂t2 = (λ+2μ)grad(Δφ)-(F(t)/(4π)(graddiv(e/r)
L'équation des ondes P devient ∂2φ/∂t2 = VP2Δφ-(F(t)/4πρ)div(e/r). Elle comporte un terme source proportionnel à (dive/r) = grad(1/r).e = -1/r2er.e = -cosθ/r2
L'amplitude du potentiel de déplacement des particules et du déplacement en champ lointain est ainsi proportionnelle à cosθ.
En écrivant le potentiel du déplacement des particules sous la forme de la divergence d'un vecteur Φ(r)e, on sépare les dépendances en r et θ.

u = gradφ(r,θ) 2φ/∂t2 = VP2Δφ-(F(t)/4πρ)div(e/r) φ(r,θ) = div(Φ(r)e) 2Φ/∂t2 = VP2ΔΦ-F(t)/(4πρr) 2(rΦ)/∂t2 = VP22(rΦ)/∂r2-F(t)/(4πρ)
F(t) = Fe-iωt, k = ω/VP 2(rΦ)/∂r2+k2(rΦ) = F/(4πρVP2) Φ = F/(4πρVP2k2)(1-eikr)/r φ = gradΦ(r).e = cosθ∂Φ/∂r = -Fcosθ/(4πρVP2k2)((ik-1/r)eikr/r+1/r2)
u = cosθ∂2Φ/∂r2er-1/r∂Φ/∂rsinθeθ = cosθ∂2Φ/∂r2er+1/r∂Φ/∂r(e-cosθer)
u = ( -F/(4πρVP2k2)) ((cosθ(2/r2-2ik/r-k2)er+1/r(ik-1/r)(e-cosθer))eikr/r) -(1/r3)(3cosθer-e))
ucl = Fcosθer/(4πρVP2)eikr/r ucp = -F/(4πρVP2)(e-3cosθer)((i/kr-1/(kr)2)eikr/r+1/r3)

Les équations d'équilibre avec la partie du terme source générant les ondes S s'écrivent ρrot2ψ/∂t2 = μrotΔψ+(F(t)/(4π))rotrot(e/r).
L'équation des ondes S devient ∂2ψ/∂t2 =VS2Δψ+(F(t)/(4πρ))rot(e/r). Elle comporte un terme source proportionnel à rot(e/r) = grad(1/r)∧e = -(1/r2)er∧e = sinθeξ/r2
L'amplitude du potentiel de déplacement des particules et du déplacement en champ lointain est ainsi proportionnelle à sinθ.
En écrivant le potentiel du déplacement des particules sous la forme du rotationnel d'un vecteur Ψ(r)e, on sépare les dépendances en r et θ.

u = rotψ(r,θ) 2ψ/∂t2 = VS2Δψ+(F(t)/4πρ)rot(e/r) ψ(r,θ) = -rot(Ψ(r)e) 2Ψ/∂t2 = VS2ΔΨ-F(t)/(4πρr) 2(rΨ)/∂t2 = VS22(rΨ)/∂r2-F(t)/(4πρ)
F(t) = Fe-iωt, k = ω/VS 2(rΨ)/∂r2+k2(rΨ) = F/(4πρVS2) Ψ = F/(4πρVS2k2)(1-eikr)/r ψ(r,θ) = -gradΨ(r)∧e = sinθ∂Ψ/∂reξ = -Fsinθeξ/(4πρVS2k2)((ik-1/r)eikr/r+1/r2)
u = (1/r)(2cosθer∂Ψ/∂r-sinθeθ∂(r∂Ψ/∂r)/∂r) = (1/r)(2(e+sinθeθ)∂Ψ/∂r-sinθeθ(∂Ψ/∂r+r∂2Ψ/∂r2))
u = -F/(4πρVS2k2) (1/r)((2(e+sinθeθ)(ik-1/r)-sinθeθ(-rk2-ik+1/r))eikr/r +(1/r2)(2e+3sinθeθ))
ucl = -Fsinθeθ/(4πρVS2)eikr/r ucp = F/(4πρVS2)(e-3cosθer) ((i/kr-1/(kr)2)eikr/r+1/r3)

Voici comment les termes de déplacement en champ proche et lointains contribuent au déplacement des particules et à la polarisation des ondes P et S aux distances de la source voisines de la longueur d'onde :
force ponctuelle : champ proche et lointain - force ponctuelle : onde P - force ponctuelle : onde S -
force.m

Anisotropie : cas isotrope transverse

Lorsque le milieu de propagation n'est pas isotrope, les caractéristiques des ondes sismiques sont plus complexes : les vitesses de propagation dépendent de la direction de propagation, les polarisations des ondes ne sont plus parfaitement longitudinales ou transverses, les vitesses des ondes SH et SV ne sont plus égales (biréfringence).
Un milieu stratifié horizontalement composé de couches isotropes fines devant la longueur d'onde a un comportement anisotrope du type millefeuilles : si on le comprime verticalement, les couches élastiques sont en série ; si on le comprime horizontalement, les couches élastiques sont en parallèle. Un tel milieu composite peut être remplacé par un milieu équivalent homogène au comportement anisotrope caractérisé par une isotropie transverse. Les relations contrainte-déformation font intervenir cinq coefficients élastiques au lieu de deux pour le milieu isotrope.
Pour établir les coefficients élastiques anisotropes équivalents dans ce cas, on considère des états de contrainte-déformation qui permettent de calculer un coefficient anisotrope à partir de moyennes des coefficients isotropes.

cas isotrope : 2 coef. élastiques λ et μ isotropie transverse (axe principal Oz) : 5 coef. élastiques exemple pour des argiles (x109Pa)
σxx = c11εxx+c12εyy+c12εzz
σxy = (c11-c12xy
c11 = λ+2μ , c12 = λ
σxx = c11εxx+c12εyy+c13εzz , σyy = c12εxx+c11εyy+c13εzz
σzz = c13εxx+c13εyy+c33εzz
σyz = 2c44εyz , σzx = 2c44εzx
σxy = 2c66εxy , c66 = (c11-c12)/2
c11 = 45 , c12 = 10 , c13 = 8
c33 = 28
c44 = 11
c66 = (c11-c12)/2 = 17.5

σ , ε appliquées sur milieu finement stratifié d'épaisseur H couches horizontales d'épaisseur hi, pi = hi/H modules anisotropes équivalents
couches comprimées en série verticalement : σzz ≠ 0
sans déformation latérale : εxx = εyy = 0
σzzi = σzz , εzzi = σzz/(λi+2μi)
εzz = ∑piεzzi = σzz∑pi/(λi+2μi)
σxxi = σyyi = λiεzzi
σzz = c33εzz
c33 = 1/∑pi/(λi+2μi)
couches comprimées en parallèle horizontalement selon x : σxx ≠ 0
sans déformation selon y : εyy = 0
sans déformation verticale moyenne : εzz = ∑piεzzi = 0
déformation selon x uniforme : εxxi = εxx
σxxi = (λi+2μixxiεzzi
σxx = ∑piσxxi = εxx∑pii+2μi)+∑piλiεzzi [1]
σzzi = σzz , σzz = λiεxx+(λi+2μizzi [2]
[2] ⇒ σzz∑pi/(λi+2μi) = εxx∑piλi/(λi+2μi) + ∑piεzzi = εxx∑piλi/(λi+2μi) [3]
[2,3] ⇒ ∑piλiεzzi = εxx((∑piλi/(λi+2μi))2/∑pi/(λi+2μi)-∑piλi2/(λi+2μi)) [4]
σxx = c11εxx
[1,4] ⇒ c11 = c33(∑piλi/(λi+2μi))2+4∑piμiii)/(λi+2μi)
couches cisaillées en série horizontalement : σyz ≠ 0 σyz = 2μiεyzi
εyz = ∑piεyzi = (σyz/2)∑pii
σyz = 2c44εyz
c44 = 1/∑pii
couches cisaillées en parallèle horizontalement : σxy ≠ 0
déformation uniforme : εxyi = εxy
σxyi = 2μiεxy
σxy = ∑piσxyi = 2∑piμiεxy
σxy = 2c66εxy
c66 = ∑piμi , c12 = c11-2c66
couches comprimées en parallèle horizontalement selon x : σxx ≠ 0
sans contrainte verticale : σzz = 0
sans déformation selon y : εyy = 0
déformation selon x uniforme : εxxi = εxx
σxxi = (λi+2μixxiεzzi
σzzi = σzz = 0 , λiεxx+(λi+2μizzi = 0
εzz = ∑piεzzi = -εxx∑piλi/(λi+2μi)
σzz = c13εxx+c33εzz = 0
c13 = c33(∑piλi/(λi+2μi))

Voici l'illustration de ces différents états de déformations et contraintes et les modules anisotropes équivalents obtenus pour une stratification alternant en proportion égale des couches sédimentaires rapides et lentes :
modules anisotropes équivalents
anisomod.m

ondes planes SH (polarisation horizontale) ondes planes quasi-P (polarisation α quasi-longitudinale : α ≈ θ)
ondes planes quasi-SV (polarisation β quasi-transverse : β ≈ η-π/2)
uy(x,z,t) = Uyei(kxx+kzz-ωt)
ρω2 = (c66kx2+c44kz2)
la relation de dispersion est une ellipse de demi-axes (c66/ρ) et (c44/ρ)
u(x,z,t) = (Ux,0,Uz)ei(kxx+kzz-ωt)
ρω2Ux = (c11kx2+c44kz2)Ux+(c13+c44)kxkzUz
ρω2Uz = (c13+c44)kxkzUx+(c44kx2+c33kz2)Uz
(ρω2-c11kx2-c44kz2)(ρω2-c44kx2-c33kz2)-((c13+c44)kxkz)2 = 0
VSH2 = (1/ρ)(c66sin2η+c44cos2η) VQP2 = (1/2ρ)(c44+c11sin2θ+c33cos2θ+D)
VQS2 = (1/2ρ)(c44+c11sin2θ+c33cos2θ-D)
D2 = ((c11-c44)sin2θ-(c33-c44)cos2θ)2+(c13+c44)2sin2
tgα = Ux/Uz = (c13+c44)sin2θ/2(ρVQP2-c11sin2θ-c44cos2θ) = (c13+c44)sin2θ/((c44-c11)sin2θ+(c33-c44)cos2θ+D)
tgβ = Ux/Uz = (c13+c44)sin2θ/2(ρVQS2-c11sin2θ-c44cos2θ) = (c13+c44)sin2θ/((c44-c11)sin2θ+(c33-c44)cos2θ-D)
tgα.tgβ = ((c13+c44)sin2θ)2/(((c44-c11)sin2θ+(c33-c44)cos2θ)2-D2) = -1 ⇒ β = α-π/2
vitesse de groupe ⇔ phase stationnaire dans ∫∫U(kx,kz)ei(kxx+kzz-ω(kx,kz)t)dkxdkz pour avoir interférence constructive des ondes harmoniques en (x,z,t)
∂(kxx+kzz-ω(kx,kz)t)/∂kx = 0 , x-(∂ω/∂kx)t = 0 , VGx = ∂ω/∂kx
∂(kxx+kzz-ω(kx,kz)t)/∂kz = 0 , z-(∂ω/∂kz)t = 0 , VGz = ∂ω/∂kz
VG = gradkω(kx,kz) ⇒ ek.VG = ek.gradkω(kx,kz) = dω/dk = Vphase ⇒ VG ≥ Vphase
VGx = (c66/ρω)kx , VGz = (c44/ρω)kz
VG2 = (1/ρ2ω2)((c66kx)2+(c44kz)2)
sinηg = VGx/VG = (c66kx)/((c66kx)2+(c44kz)2)½
cosηg = VGz/VG = (c44kz)/((c66kx)2+(c44kz)2)½
VG2 = 1/(ρ(sin2ηg/c66+cos2ηg/c44))
VGx = ∂(kV)/∂kx = V∂k/∂kx+k∂V/∂kx
k = (kx2+kz2)½ , ∂k/∂kx = kx/k = sinθ
kx = ksinθ , 1 = ∂k/∂kxsinθ+kcosθ∂θ/∂kx , ∂θ/∂kx = cosθ/k
VGx = Vsinθ+cosθdV/dθ , VGz = Vcosθ-sinθdV/dθ
VG2 = V2+(dV/dθ)2
tgθg = VGx/VGz = (tgθ+(1/V)dV/dθ)/(1-(tgθ/V)dV/dθ)
paramètres d'anisotropie de Thomsen utilisés dans le traitement des données de sismique réflexion : αT , βT , γT , δT , εT
η = 0 , VSH2 = c44/ρ = βT2 , η = 90° , VSH2 = c66
γT = (VSH2(90)-VSH2(0))/2VSH2(0) = (c66-c44)/2c44
VSH2 = βT2((c66/c44)sin2η+cos2η) = βT2(1+2γTsin2η)
VG2 = βT2((1+2γT)/(1+2γTcos2ηg))
θ = 0 , D = (c33-c44) , VQP2 = c33/ρ = αT2 , VQS2 = c44/ρ = βT2
θ = 90° , D = (c11-c44) , VQP2 = c11/ρ , VQS2 = c44/ρ = βT2
εT = (VQP2(90)-VQP2(0))/2VQP2(0) = (c11-c33)/2c33
δT = (1/2αT)(d2VQP/dθ2)(0) = ((c13+c44)2-(c33-c44)2)/(2c33(c33-c44))
pour l'exemple des argiles (ρ=2380 kg/m3) : αT = 3.43 km/s, βT = 2.15 km/s , γT = 0.3 , δT = 0.07 , εT = 0.3
faible anisotropie
γT<<1 , γT ≈ (VSH(90)-VSH(0))/VSH(0)
VSH ≈ βT(1+γTsin2η)
εT<<1 , δT<<1 , εT ≈ (VQP(90)-VQP(0))/VQP(0)
VQP ≈ αT(1+δTsin2θcos2θ+εTsin4θ) = αT(1+δTsin2θ+(εTT)sin4θ)
VQS ≈ βT(1+(αT2T2)(εTT)sin2θcos2θ)

Voici les vitesses de phase et de groupe et les polarisations obtenues avec les valeurs du tableau ci-dessus, ainsi que les groupes d'ondes SH et QP pour différents éventails d'angle de phase :
vitesses anisotropes - onde SH anisotrope - onde QP anisotrope
anisopro.m

Pour le calcul des coefficients de réflexion-transmission QP-QSV, on écrit les composantes du déplacement des particules en fonction de l'angle de polarisation α pour les ondes QP et de l'angle β = angle de polarisation QSV + π/2 pour les ondes QS (comme dans le cas isotrope). La loi de Snell s'applique toujours mais ne donne pas explicitement les différents angles d'incidence puisque les vitesses dépendent de ces angles. La relation de dispersion permet d'obtenir les angles θ(p) et η(p), ou plus généralement kzP = ωqP(p) et kzS = ωqS(p), imaginaires pour les ondes évanescentes.

Loi de Snell : p = sinθ1/VQP11) = sinθ2/VQP22) = sinη1/VQS11) = sinη2/VQS22) ;
qP1(p) = (cosθ1/VQP1)(p) ; qP2(p) = (cosθ2/VQP2)(p) ; qS1(p) = (cosη1/VQS1)(p) ; qS2(p) = (cosη1/VQS2)(p) ;
Relation de dispersion : (ρ-c11p2-c44q2)(ρ-c44p2-c33q2)-((c13+c44)pq)2 = 0
Aq4+Bq2+C = 0 ; A = c33c44 ; B = -c33(ρ-c11p2)-c44(ρ-c44p2)-(c13+c44)2p2 ; C = (ρ-c11p2)(ρ-c44p2) ; Δ = (B2-4AC)½
qP = ((-B-Δ)/2A)½ ; qS = ((-B+Δ)/2A)½ ;
Polarisation : α = Arctg((c13+c44)pqP/(ρ-c11p2-c44(qP)2)) ; β = π/2+Arctg((c13+c44)pqS/(ρ-c11p2-c44(qS)2))
QP incidente descendante + réfléchie montante
UI = (sinα1,cosα1) , URP = RP(sinα1,-cosα1)
kzI = ωqP1(p) , kzRP = -ωqP1(p)
QSV réfléchie montante
URS = RS(cosβ1,sinβ1)
kzRS = -ωqS1(p)
QP transmise descendante
UTP = TP(sinα2,cosα2)
kzTP = ωqP2(p)
QSV transmise descendante
UTS = TS(-cosβ2,sinβ2)
kzTS = ωqS2(p)
Ux = sinα1(1+RP)+ cosβ1RS = sinα2TP-cosβ2TS
Uz = cosα1(1-RP)+ sinβ1RS = cosα2TP+ sinβ2TS
Σzz = (c13sinα1p+c33cosα1qP1)(1+RP) + (c13cosβ1p-c33sinβ1qS1)RS = (d13sinα2p+d33cosα2qP2)TP + (-d13cosβ2p+d33sinβ2qS2)TS
Σzx = c44(sinα1qP1+cosα1p)(1-RP) + c44(-cosβ1qS1+sinβ1p)RS = d44(sinα2qP2+cosα2p)TP + d44(-cosβ2qS2+sinβ2p)TS
écriture matricielle : AR = B ; R = [RP , RS , TP , TS]t ; B = [-sinα1 , -cosα1 , -(c13sinα1p+c33cosα1qP1) , c44(sinα1qP1+cosα1p)]t
Asinα1cosβ1-sinα2cosβ2
-cosα1sinβ1-cosα2-sinβ2
(c13sinα1p+c33cosα1qP1) (c13cosβ1p-c33sinβ1qS1) - (d13sinα2p+d33cosα2qP2) (d13cosβ2p-d33sinβ2qS2)
c44(sinα1qP1+cosα1p) c44(cosβ1qS1-sinβ1p) d44(sinα2qP2+cosα2p) d44(-cosβ2qS2+sinβ2p)

Voici les coefficients obtenus pour une interface entre un milieu isotrope lent et un milieu isotrope transverse rapide dans le cas où une onde P est incidente dans le milieu isotrope :
coefficients de réflexion-transmission anisotropes
anisopro.m

Winterstein D.F. and G. S. De, 2001, VTI documented montrent l'effet de l'isotropie transverse sur des ondes S enregistrées en sismique de puits 9 composantes.
Johnston, J. E. and Christensen, N. I. 1995, Seismic anisotropy of shales mesurent les vitesses de phase et de groupe sur des argiles.
Godfrey, N.J. ; Christensen, N. I. and Okaya, D.A. 2000, Anisotropy of schists: Contribution of crustal anisotropy to active source seismic experiments and shear wave splitting observations montrent l'importance de la contribution des schistes pour l'anisotropie dans la croûte.
Picotti, S., A. Vuan, J. M. Carcione, H. J. Horgan, and S. Anandakrishnan (2015), Anisotropy and crystalline fabric of Whillans Ice Stream (West Antarctica) inferred from multicomponent seismic data
Graebner M., 1992, Plane-wave reflection and transmission coefficients for a transversely isotropic solid

Atténuation

déplacement
eq. équilibre
densité énergie cinétique (Jm-3) densité énergie élastique (Jm-3) densité puissance transportée Pt (Wm-2) intensité I (Wm-2)
u(x,t)
ρ∂2u/∂t2 = ∂σ/∂x
wc = (ρ/2)(∂u/∂t)2
∂wc/∂t = (∂σ/∂x)(∂u/∂t)
we = ∫σdε = (E/2)ε2
∂we/∂t = σ∂/∂x(∂u/∂t)
∂(wc+we)/∂t = ∂/∂x(σ∂u/∂t)
Pt = -σ∂u/∂t
(1/T)∫0TPtdt
u(x,t) = Ucos(kx-ωt) wc = (ρ/2)ω2U2sin2(kx-ωt) we = (ρ/2)V2k2U2sin2(kx-ωt) Pt = ρVω2U2sin2(kx-ωt) = V(wc+we) (ρ/2)Vω2U2
u = gradφ(x,t)
ρ∂2u/∂t2 = -gradp
wc = (ρ/2)(∂u/∂t)2
∂wc/∂t = -gradp.∂u/∂t
we = (K/2)(divu)2
∂we/∂t = -pdiv∂u/∂t
∫dx∂(wc+we)/∂t = -∫dxdiv(p∂u/∂t) = -∫(p∂u/∂t).dS
Pt = p∂u/∂t
(1/2)ℜ(p∂u*/∂t)

2π/facteur de qualité = énergie dissipée par cycle / énergie élastique maximum cycle = λ amplitude déplacement vitesse complexe VQ = V(1-i/2Q)
2π/Q = -(δw)cycle/wmax π/Q = -(dU/dx)λ/U U(x) = Ue-αx
α = π/Qλ = k/2Q = ω/2QV
u(x,t) = Ue-αxcos(kx-ωt) u(x,t) = Ue-αxeiω(x/V-t) = Ueiω(x(1+i/2Q)/V-t)
u(x,t) = Ueiω(x/VQ-t)

viscoélasticité linéaire énergie dissipée par cycle énergie élastique maximum facteur de qualité
σ = σmcos(ωt+χ)
ε = εmcosωt
(δw)cycle = ∫0Tσdε = πσmεmsinχ wmax = ∫0T/4σdε ≈ -(1/2)σmεmcosχ 1/Q = tgχ ≈ χ

modèle viscoélastique , module élastique E , viscosité η nombre d'onde complexe k+iα faible atténuation , 1/Q << 1 1/Q < 1 , dispersion V(ω) = ω/k
E et η en parallèle, temps relaxation τ = η/E
σ = Eε+η∂ε/∂t = E(1-iωτ)ε = Mε
M = E(1+ω2τ2)½e-iχ , 1/Q = tgχ = ωτ
ρ∂2u/∂t2 = E(1-i/Q)∂2u/∂x2
ρω2 = E(1-i/Q)(k+iα)2 , k02 = ρω2/E
k02 = (1-i/Q)(k+iα)2
k+iα = k0(1+i/2Q)
k = k0 , V = V0 = ω/k0
α = k0/(2Q) = ω2τ/(2V0)
k22 = k02/(1+1/Q2)½ ≈ k02(1-1/2Q2)
k22 = k02/(1+1/Q2) ≈ k02(1-1/Q2)
V(ω) ≈ V0(1+3/(8Q2)) = V0(1+3ω2τ2/8) , α = k0/(2Q)
(E et η en parallèle) en série avec E' , σ = M''ε
1/M'' = 1/M+1/E' = 1/E(1-iωτ)+1/E' = (1/E+1/E'-iωτ/E')/(1-iωτ)
1/E'' = 1/E+1/E' , τ'' = τE''/E' , M'' = E''(1-iωτ)/(1-iωτ'') = ||M''||e-i(χ-χ'')
1/Q = tg(χ-χ'') = ω(τ-τ'')/(1+ω2ττ'')
d(1/Q)/dω = 0 pour ωr = (ττ'') et (1/Q)max = (τ-τ'')/2(ττ'')½
ρω2 = M''(k+iα)2 , k''02 = ρω2/E''
k22 = ρω2/||M''|| = k''02(1+ω2τ''2)½/(1+ω2τ2)½
k22 = ρω2ℜ(1/M'') = k''02(1+ω2ττ'')/(1+ω2τ2)
k2 = ω2/V(ω)2 = (k''02/2)((1+ω2τ''2)½/(1+ω2τ2)½+(1+ω2ττ'')/(1+ω2τ2))

modèle à Q constant (Kjartansson, 1979, JGR, 84, 4737)
M = M0(iω/ω0) = M0(|ω/ω0|)eiπγsign(ω)
1/Q = tg(πγ)
ρω2 = M0(|ω/ω0|)eiπγsign(ω)(k+iα)2
k22 = k02(|ω/ω0|)-2γ
k22 = k02(|ω/ω0|)-2γcos(πγ)
k = k0(|ω/ω0|)cos(πγ/2)
V(ω) = V0(|ω/ω0|)γ/cos(πγ/2)
α = k0(|ω/ω0|)sin(πγ/2) = ktg(πγ/2)
dispersion pour 1/Q << 1
V(ω)/V0 ≈ (|ω/ω0|)1/πQ = eLog(|ω/ω0|)/πQ
V(ω)/V0 ≈ 1+Log(|ω/ω0|)/πQ

Sun, L.F., B. Milkereit, and D. R. Schmitt, 2009 Measuring velocity dispersion and attenuation in the exploration seismic frequency band mettent en évidence la dispersion et l'atténuation sur des données vibrosismiques.

Poroélasticité et vitesses dans les roches saturées en fluide

La porosité Φ est la proportion volumique d'espace poreux ΩΦ dans un volume élémentaire Ω. Seuls les pores interconnectés dans lesquels un fluide peut s'écouler interviennent dans Φ, qui atteint 30% pour certains grès. Pour une porosité saturée par un fluide, le volume de fluide Ωf est égal au volume des pores. On considère une déformation volumique ε de la roche due à une variation de la pression de confinement -σ. L'augmentation ζ du contenu volumique en fluide est due à l'augmentation du volume des pores et à la compression volumique du fluide.

Φ = ΩΦ/Ω = Ωfε = δΩ/Ω ζ = (δΩΦ-δΩf)/Ω = Φ(δΩΦΦ-δΩff) = Φ(εΦf)

La théorie de la poroélasticité linéaire isotrope exprime ε et ζ comme une combinaison linéaire de σ et pf , la variation de pression de fluide dans les pores. La signification des coefficients de proportionnalité α (coefficient de Biot-Willis) et B (coefficient de Skempton) apparaît en considérant les cas d'une roche drainée (pression de fluide constante, incompressibilité K), non-drainée (contenu en fluide constant, incompressibilité Ksat) et d'une pression différentielle nulle (pressions de confinement et de fluide égales). Dans ce dernier cas, la déformation homogène est celle de la roche sans pore d'incompressibilité K0 et la déformation volumique εΦ des pores est égale à la déformation volumique ε de la roche : on obtient ainsi la relation de Gassmann entre Ksat, K0 , K , Kf (l'incompressibilité du fluide) et Φ. Cette relation permet de déterminer la vitesse de propagation des ondes de compression dans les roches poreuses saturées en fluide à basse fréquence (f < 100Hz).

poroélasticité linéaire ε = (1/K)(σ+αpf)
ζ = (α/K)(σ+pf/B)
incompressibilité coef. poroélastiques
pf = 0 ε = (1/K)σ K α = ζ/ε
ζ = 0 ε = ((1-αB)/K))σ K = (1-αB)Ksat B = -pf
σ = -pf ε = ((1-α)/K))σ = σ/K0 K = (1-α)K0 .
εΦ = ζ/Φ+εf = ε
εf = -pf/Kf
(α/ΦK)(1-1/B)+1/Kf = 1/K0 1/B = 1-(ΦK/α)(1/K0-1/Kf) = α/(1-K/Ksat)
relation de Gassmann 1/Ksat = (1-α)/K+α/K(1-1/(1-(ΦK/α)(1/K0-1/Kf)))
1/Ksat = 1/K0+Φ/(ΦK/α+1/(1/Kf-1/K0))
1/(1/Ksat-1/K0) = K/α+1/(Φ(1/Kf-1/K0)) = 1/(1/K-1/K0)+1/(Φ(1/Kf-1/K0))
densité et vitesses ρ = (1-Φ)ρ0+Φρf VP2 = (Ksat+4μ/3)/ρ VS2 = μ/ρ

Wang Z., 2001, Fundamentals of seismic rock physics montre comment la relation de Gassmann est utilisée dans l'étude sismique des réservoirs.
Smith, T.M., C. H. Sondergeld, and C. S. Rai, 2003, Gassmann fluid substitutions: A tutorial

Réflexion 1D sur gradient de vitesse (ρ=cte)

couchesdéplacementscontinuité u et σ en z=0 et z=H
z<0 V=V0 u(z) = eiωz/V0+Re-iωz/V0 1+R = AV0n+BV0m
(iω/V0)(1-R) = a(AnV0n-1+BmV0m-1)

TeiωH/V1 = AV1n+BV1m
(iω/V1)TeiωH/V1 = a(AnV1n-1+BmV1m-1)

⇒ B/A = -V1n-m(n-iω/a)/(m-iω/a)

0≤z≤H V(z)=V0+az -ρω2u(z) = ∂/∂z(ρV(z)2∂u/∂z)
u(z) = A(V0+az)n+B(V0+az)m
n,m = -1/2±(1/4-ω2/a2)½
n-m = 2(1/4-ω2/a2)½
z>H V=V1 u(z) = Teiωz/V1
réflexion sur gradient vitesse R = -(V0n-m-V1n-m)/((n+iω/a)/(n-iω/a)V0n-m-(m+iω/a)/(m-iω/a)V1n-m) = 1/(2iω/a+(n-m)(V0n-m+V1n-m)/(V0n-m-V1n-m))

Liner, C. and Bodmann, B. (2010). The Wolf ramp: Reflection characteristics of a transition layer
Chapman, N.R., J. F. Gettrust, R. Walia, D. Hannay, G. D. Spence, W. T. Wood, and R.D. Hyndman , 2002, High-resolution, deep-towed, multichannel seismic survey of deep-sea gas hydrates off western Canada
Wolf, A., 1937, The reflection of elastic waves from transition layers of variable velocity Read More: http://library.seg.org/doi/ref/10.1190/1.3476312

Réflexion-transmission d'une onde monochromatique sur une couche ou une stratification

La réflexion (transmission) d'une onde plane monochromatique SH sur (à travers) une couche d'épaisseur H incluse dans un espace homogène est formée par interférence entre la réflexion (transmission) primaire en entrant dans (à travers) la couche et les réflexions multiples qui sortent de la couche dans la direction de réflexion (transmission). Si Re et Te=(1+Re) sont les coefficients de réflexion et transmission en entrant dans la couche, Rs=(-Re) et Ts=(1-Re) en sortant de la couche, et kz le nombre d'onde vertical dans la couche, les coefficients RH et TH sont donnés par la somme des termes de la série des multiples successifs. Pour Re<<1, il suffit de prendre en compte le premier terme de la série :

RH1 = Re+1 réflexion dans la couche , Re<<1 , TeTs = 1-Re2 ≈ 1 RHn = Re+n réflexions dans la couche , RH = RHn (n→∞) RH1 , RH , TH1, TH
RH1 = Re+TeeikzHRseikzHTs = Re(1-(1-Re2)e2ikzH) ≈ Re(1-e2ikzH)
RH1 = 0 pour kzH=nπ ou λz = 2H/n
RH1 = 2Re pour kzH=(2n+1)π/2 ou λz = 4H/(2n+1)
RHn = Re+TeRsTse2ikzH(1+∑(RseikzH)2n)
RH = Re(1-(1-Re2)e2ikzH/(1-Re2e2ikzH)) = Re(1-e2ikzH)/(1-Re2e2ikzH)
réflexion sur une couche
TH1 = TeeikzHTs = (1-Re2)eikzH ≈ eikzH THn = TeeikzHTs(1+∑(RseikzH)2n)
TH = (1-Re2)eikzH/(1-Re2e2ikzH)

Alternativement, on peut considérer qu'il y a une onde montante (déplacement uy pour une onde SH, potentiel φ pour une onde P acoustique) et une onde descendante dans la couche, somme respectivement de toutes les ondes montantes et descendantes dans cette couche. Pour une stratification, les ondes montantes et descendantes dans deux couches ayant une interface commune sont liées par des relations de continuité. On obtient les coefficients de réflexion-transmission sur la stratification en remontant depuis l'interface avec le demi-espace inférieur et en utilisant itérativement ces relations (le tableau se lit de bas en haut). Inversement, on obtient les coefficients de réflexion sur les interfaces à partir des rapports onde montante/descendante dans chaque couche. Le calcul des coefficients de réflexion-transmission sur la stratification pour différentes fréquences permet d'obtenir la trace sismique en temps par transformée de Fourier inverse (sommation des réponses harmoniques, sujette aux problèmes numériques d'échantillonnage et de périodicité). Pour tenir compte de la surface libre (R = 1), il suffit d'utiliser l'onde réfléchie par la stratification et par la surface libre comme onde descendante, et ainsi de suite pour obtenir les réflexions multiples successives dues à la surface libre.

couche
Rinterface
ondes descendantes
Dj = Tj+1Dj+1eikzjhj-Rj+1Mje2ikzjhj
Dj+1 = (Dje-ikzjhj+Rj+1Mjeikzjhj)/Tj+1
ondes montantes
Mj+1 = Rj+1Dj+1+(1-Rj+1)Mjeikzjhj
Mj+1 = (Rj+1Dje-ikzjhj+Mjeikzjhj)/Tj+1
coef. réflexion sur stratification
Rh1..j = Mj+1/Dj+1
Rj+1 = (Rh1..je-ikzjhj-Rh1..j-1eikzjhj)/(e-ikzjhj-Rh1..j-1Rh1..jeikzjhj)
coef. transmission par stratification
Th1..j = D0/Dj+1
couche 3
R3
D3 = (e-ikz2h2+R3Rh1eikz2h2)D2/T3 M3 = (R3e-ikz2h2+Rh1eikz2h2)D2/T3 Rh12 = M3/D3
R3 = (Rh12e-ikz2h2-Rh1eikz2h2)/(e-ikz2h2-Rh1Rh12eikz2h2)
Th12 = D0/D3 = Th1D2/D3
couche 2
R2
D2 = (e-ikz1h1+R2R1eikz1h1)D1/T2 M2 = (R2e-ikz1h1+R1eikz1h1)D1/T2 Rh1 = M2/D2
R2 = (Rh1e-ikz1h1-R1eikz1h1)/(e-ikz1h1-R1Rh1eikz1h1)
Th1 = D0/D2 = T1D1/D2
couche 1
R1
D1 M1 R1 = M1/D1 T1 = D0/D1
1/2 espace D0 M0 = 0 montante+descendante
stratif.m

Pour des ondes P-SV, il y a deux ondes (potentiels φ et &psiy) montantes et deux ondes descendantes dans chaque couche. On obtient les coefficients de réflexion-transmission sur la stratification en calculant la matrice de transfert de la stratification, produit des matrices de transfert dans chaque couche. Sur l'interface inférieure de la stratification avec le demi-espace sous-jacent, les ondes montantes sont les ondes descendantes réfléchies. Dans le demi-espace sus-jacent les ondes descendantes sont les ondes incidentes choisies (P et/ou SV) et les rapports montantes/descendantes donnent les coefficients de réflexion. Si la stratification est sous l'eau (sismique marine), il faut multiplier la matrice de transfert précédente par la matrice calculée à l'interface liquide-solide.

ondes descendantes P et SV
DPj = (TPPj+1DPj+1+TSPj+1DSj+1+RPPj+1MPjeikzPjhj+RSPj+1MSjeikzSjhj)eikzPjhj
DSj = (TPSj+1DPj+1+TSSj+1DSj+1+RPSj+1MPjeikzPjhj+RSSj+1MSjeikzSjhj)eikzSjhj
ondes montantes P et SV
MPj+1 = RPPj+1DPj+1+RSPj+1DSj+1+TPPj+1MPjeikzPjhj+TSPj+1MSjeikzSjhj
MSj+1 = RPSj+1DPj+1+RSSj+1DSj+1+TPSj+1MPjeikzPjhj+TSSj+1MSjeikzSjhj
(TPPj+1TSSj+1-TPSj+1TSPj+1)DPj+1

(TSPj+1TPSj+1-TSSj+1TPPj+1)DSj+1

TSSj+1e-ikzPjhj -TSPj+1e-ikzSjhj (TSPj+1RPSj+1-TSSj+1RPPj+1)eikzPjhj (TSPj+1RSSj+1-TSSj+1RSPj+1)eikzSjhj) DPj
DSj
MPj
MSj
TPSj+1e-ikzPjhj -TPPj+1e-ikzSjhj (TPPj+1RPSj+1-TPSj+1RPPj+1)eikzPjhj (TPPj+1RSSj+1-TPSj+1RSPj+1)eikzSjhj)
MPj+1

MSj+1

RPPj+1 RSPj+1 TPPj+1eikzPjhj TSPj+1eikzSjhj DPj+1
DSj+1
MPj
MSj
RPSj+1 RSSj+1 TPSj+1eikzPjhj TSSj+1eikzSjhj
1 (0)
0 (1)
RPPh1.n (RSPh1.n)
RPSh1.n (RSSh1.n)
DM(4,4) = ∏ matrices transfert dans couches 1 à n 1
0
RPP1
RPS1
0
1
RSP1
RSS1
DP1
DS1
cas liquide-solide
DPj = (TPPj+1DPj+1+RPPj+1MPjeikzPjhj+RSPj+1MSjeikzSjhj)eikzPjhj
MPj+1 = RPPj+1DPj+1+TPPj+1MPjeikzPjhj+TSPj+1MSjeikzSjhj
1
RPPeau.h1.n
DMliq-sol(2,4) DM(4,4) 1
0
RPP1
RPS1
0
1
RSP1
RSS1
DP1
DS1

Voici les potentiels de déplacement calculés avec cette méthode en utilisant un modèle de vitesses-densité stratifié type Mer du Nord pour des paramètres d'ondes planes correspondant à des angles d'incidence dans l'eau θ = 0°, 15° et 30° dans les cas correspondant à la sismique réflexion terrestre (réflexion P-P et SV-SV sans eau) , à la sismique réflexion marine (réflexion P-P dans l'eau) et à la sismique réflexion marine avec enregistrement au fond de l'eau (réflexion P et SV sous l'eau pour une onde P incidente dans l'eau).
Modèle vitesses-densité Mer du Nord - Réflexion-transmission Mer du Nord - Stratification P-SV incidence 0° - Stratification P-SV incidence 15° - Stratification P-SV incidence 30°
stratifpsv.m

Liu Y. and D. R. Schmitt, 2003, Amplitude and AVO responses of a single thin bed
Chapman, C. H., 2003, Yet another elastic plane-wave, layer-matrix algorithm donne un programme plus élaboré, utilisable pour des incidences supérieures à l'incidence critique.

Monitoring sismique en milieu élastique stratifié

La sismique réflexion est utilisable pour détecter les changements de propriétés élastiques (vitesses, densité, atténuation) dans le sous-sol. Il faut répéter des expériences identiques et détecter les changements dans les enregistrements (temps de trajet, amplitudes, forme d'ondes, fréquences ....) dus aux modifications du milieu sous-jacent. La méthode ci-dessus est utilisable pour simuler les changements correspondant à des modèles différents. Même si l'on ne modifie les propriétés élastiques que dans une couche fine, correspondant par exemple à un réservoir en production, toutes les arrivées ayant traversé cette couche sont affectées. L'exemple suivant montre ce qui se passe pour des ondes planes P et SV se propageant dans un modèle simple comportant une couche de 40 m d'épaisseur à 1000 m de profondeur dans laquelle la vitesse VP diminue de 3200 à 2800 m/s, les autres paramètres étant inchangés.

Modèle vitesses-densité monitoring - Réflexion-transmission monitoring - Monitoring P-SV incidence 0° - Monitoring P-SV incidence 15° - Monitoring P-SV incidence 30°
monitorpsv.m

Réflexion-transmission d'une onde SH sur une stratification cyclique

Si l'onde est incidente sur une stratification cyclique constituée d'une alternance régulière de 2 couches telles que les coefficients de réflexion-transmission à chaque interface soient alternativement (Re,Te) et (Rs,Ts) , l'effet des réflexions multiples de même ordre sur les différentes interfaces est cumulatif. Pour un grand nombre d'interfaces, les modules des termes du coefficient de transmission tenant compte des premiers multiples dans chaque couche sont supérieurs à celui de la transmission directe. En supposant que kzH est le même dans les 2 couches et que la stratification comporte 2m+1 couches (2m+2 interfaces), les premiers termes du coefficient de transmission sont :

1 trajet à (TeTs)m+1 TS0 = TeeikzHTs∏(eikzHTeeikzHTs) = (1-Re2)m+1e(2m+1)ikzH transmission par stratification
transtrat.m
N1(2m+1) = 2m+1 trajets à R2 TS1 = TS0N1(2m+1)Re2e2ikzH
N2(2m+1) = N1(2m+1)+...+N1(1) trajets à R4
N1(2m) trajets à -R2,TeTs
TS2 = TS0(N2(2m+1)Re4
-N1(2m)(1-Re2)Re2)e4ikzH
N3(2m+1) = N2(2m+1)+...+N2(1) trajets à R6
2N2(2m+1)-2 trajets à -R4,TeTs
N1(2m-1) trajets à R2,(TeTs)2
TS3 = TS0(N3(2m+1)Re6
-(2N2(2m+1)-2)(1-Re2)Re4
+N1(2m-1)(1-Re2)2Re2)e6ikzH
N4(2m+1) = N3(2m+1)+...+N3(1) trajets à R8
3N3(2m+1)+N2(2m)-3 trajets à -R6,TeTs
3N2(2m+1)-N2(2)-2 trajets à R4,(TeTs)2
N1(2m-2) trajets à -R2,(TeTs)3
TS4 = TS0(N4(2m+1)Re8
-(3N3(2m+1)+N2(2m)-3)(1-Re2)Re6
+(3N2(2m+1)-N2(2)-2)(1-Re2)2Re4
+N1(2m-2)(1-Re2)3Re2)e8ikzH

Les termes successifs des coefficients de transmission et réflexion sont calculables sans avoir à dénombrer tous les trajets possibles en utilisant les récurrences reliant les ondes montantes Mj,n+1 et descendantes Dj,n+1 de part et d'autre de chaque interface j après n+1 réflexion-transmissions à celles provenant des interfaces immédiatement au dessus Dj-1,n et en dessous Mj+1,n après n réflexion-transmissions. On obtient ainsi les réponses impulsionelle et harmonique complètes de la stratification. Les réponses impulsionnelles en réflexion et transmission sont caractérisées par des codas avec une amplitude décroissant exponentiellement. Les réponses harmoniques sont caractérisées par des pics de réflectivité (trous de transmissivité) pour H=(2n+1)λz/4 .

initialisation réflexion-transmission sur interface impaire réflexion-transmission sur interface paire réponses impulsionelle et harmonique stratification
refltranstrat.m
M1,1 = Re , D1,1 = Te
M2,2 = RsTeeikzH , D2,2 = TsTeeikzH
R2n+1 = M1,2n+1 = TsM2,2neikzH
D1,2n+1 = RsM2,2neikzH
M2j+1,2n+1 = (ReD2j,2n+TsM2j+2,2n)eikzH
D2j+1,2n+1 = (TeD2j,2n+RsM2j+2,2n)eikzH
M2j,2n+2 = (RsD2j-1,2n+1+TeM2j+1,2n+1)eikzH
D2j,2n+2 = (TsD2j-1,2n+1+ReM2j+1,2n+1)eikzH
M2m+2,2n+2 = RsD2m+1,2n+1eikzH
T2n+2 = D2m+2,2n+2 = TsD2m+1,2n+1eikzH

Anstey N.A. and R. F. O'Doherty, 2002, Cycles, layers, and reflections: Part 1 expliquent l'origine des réflexions observées dans les bassins sédimentaires caractérisés par une sédimentation cyclique.
Stewart, R.R., P. D. Huddleston, and T. K. Kan, 1984, Seismic versus sonic velocities: A vertical seismic profiling study expliquent la différence entre les vitesses mesurées localement et depuis la surface dans un puits par l'effet de stratification cyclique.

Relations réflexion-transmission par une stratification

On obtient des relations entre les coefficients de réflexion-transmission pour des ondes incidentes par en haut ou par en bas sur la stratification en remarquant qu'une inversion du signe du temps inverse le sens de propagation de toutes les ondes harmoniques et remplace leurs amplitudes complexes par leurs conjuguées (les retards de phase correspondant au temps de propagation vertical dans la stratification deviennent des avances). La propagation dans la stratification des ondes réfléchies et transmises auxquelles ce retournement temporel a été préalablement appliqué reconstitue l'onde incidente initiale. En effet, les ondes refont le même chemin en sens inverse et arrivent en phase sur les plans limitant la stratification.

cas SH onde incidente en haut (a) onde incidente en bas (b) retournement (a) + propagation retournement (b) + propagation
sens propagation
haut 1 Rh 0 th Rh* 1 = RhRh*+thTh* th* 0 = Rhth*+thrh*
bas Th 0 rh 1 0 = ThRh*+rhTh* Th* 1 = Thth*+rhrh* rh*

Pour le cas d'une couche d'épaisseur H dans un milieu homogène, Re étant le coefficient de réflexion SH sur l'interface à l'entrée de la couche, on vérifie :

couche onde incidente en haut (a) onde incidente en bas (b) retournement (a) + propagation
haut 1 Rh = Re(1-e2ikzH)/(1-Re2e2ikzH) 0 th = Th RhRh*+thTh* = (Re2(1-e2ikzH)(1-e-2ikzH)+ (1-Re2)2) /((1-Re2e2ikzH)(1-Re2e-2ikzH)) = 1
bas Th = (1-Re2)eikzH/(1-Re2e2ikzH) 0 rh = Rh 1 ThRh*+rhTh* = Re(1-Re2)(eikzH(1-e-2ikzH)+e-ikzH(1-e2ikzH)) /((1-Re2e2ikzH)(1-Re2e-2ikzH)) = 0

Par ailleurs, les flux d'énergie entrant et sortant de la stratification doivent être égaux.

impédance SH x cos(incidence) onde incidente en haut onde incidente en bas
flux énergie entrant sortant entrant sortant
haut : Ia = (ρVScosη)haut Ia IaRhRh* 0 Iathth*
bas : Ib = (ρVScosη)bas 0 IbThTh* Ib Ibrhrh*
conservation flux Ia(1-RhRh*) = IbThTh* Ib(1-rhrh*) = Iathth*
relations transmission/réflexion thTh* = 1-RhRh* = 1-rhrh* = (Ib/Ia)ThTh* th = (Ib/Ia)Th RhRh* = rhrh*

Si la stratification est limitée par une surface libre, il y a réflexion totale sur celle-ci. Pour une onde incidente descendante, la réponse en réflexion de la stratification est renvoyée vers le bas. L'autocorrélation de la réponse en transmission descendante D dans le demi-espace sous la stratification est proportionnelle à la somme de l'onde incidente, de la réponse en réflexion et de sa retournée temporelle. De même, pour une onde incidente montante, l'autocorrélation de la réponse en transmission montante M dans la couche sous la surface libre est proportionnelle à la somme de l'onde incidente, de la réponse en réflexion et de sa retournée temporelle.

onde incidente transmission + 1 réflexion sur surface libre autocorrélation réflexion + retournée temporelle
descendante dans couche D = Th(1+Rh) DD* = ThTh*(1+Rh)(1+Rh)* = (Ia/Ib)Thth*(1+Rh)(1+Rh)* = (Ia/Ib)(1-RhRh*)(1+Rh)(1+Rh)* DD* ≈ (Ia/Ib)(1+Rh+Rh*)
montante dans demi-espace M = th(1+Rh) MM* = thth*(1+Rh)(1+Rh)* = (Ib/Ia)Thth*(1+Rh)(1+Rh)* = (Ib/Ia)(1-RhRh*)(1+Rh)(1+Rh)* MM* ≈ (Ib/Ia)(1+Rh+Rh*)

Réflexion-diffusion d'une onde acoustique plane sur une surface libre faiblement ondulée

Une surface libre présentant une ondulation sinusoïdale z=h(x) de la topographie autour du plan z=0 perturbe la réflexion d'une onde plane. Pour une ondulation de faible amplitude, on considère que la surface engendre une onde réfléchie sur un plan horizontal z=0 selon la loi de Snell plus une onde diffusée par les ondulations. On obtient une approximation de l'onde harmonique diffusée ps(x,z) à partir de l'onde harmonique p0(x,z) se propageant en l'absence de perturbation en supposant que Ps<<P0 et en faisant un développement limité de la condition au limite (pression nulle) qui s'applique sur l'interface ondulée.

onde harmonique incidente+réfléchie sur surface libre z=0 pression nulle sur surface libre
z = h(x) = Heikhx
onde = p0 + onde diffusée Ps<<P0
kz0H<<1
onde harmonique diffusée dans la direction θs
p0(x,z) = P0eikx0x(e-ikz0z-eikz0z)

p0(x,0) = 0

p(x,h(x)) = 0

p(x,0)+h(x)∂p/∂z(x,0) ≈ 0

p(x,z) = (p0+ps)(x,z)

ps(x,0) = -h(x)∂(p0+ps)/∂z(x,0)

ps(x,0) ≈ -h(x)∂p0/∂z(x,0)

ps(x,0) =Pseikxsx

Ps = -2ikz0HP0

kxs = (kx0+kh)

ps(x,z) = Pseikxsxeikzsz

sinθs = sinθ0+Vkh

kzs = (ω2/V2-kxs2)½

Bibliographie