Algorytm migracji MG(F-K) w monoklinalnym ośrodku model TTI

Andrzej Kostecki

The Oil and Gas Institute, Krakow

Algorithm MG(F-K) of migration in model TTI

anisotropy

Introduction

One of the most frequently encountered models in anisotropic medium is the model of Transverse Isotropy (TI). If the symmetry axis of model TI coincides with the vertical axis of right-angle coordinates, then we receive horizontal thin-layered system proposed by G. Postma [12] and called VTI (Vertical Transversely Isotropic). The application of algorithms of propaga-tion and migrapropaga-tion of „isotropic” waves in anisotropic medium VTI results in deformations and relocation of reproduced structures [6, 10], and the higher the

anisotropic parameters are, the greater their extent is. (parameters by L. Thomsen [14]).

In the domain of wavenumbers and frequencies, the al-gorithm of migration MG(F-K) was presented in anisotropic medium type VTI [11]. This algorithm uses approximative version of vertical wavenumber [5] with reference to me-dium VTI. This article presents the algorithm of migration MG(F-K) in monoclinal medium marked as TTI (Tilted Transversely Isotropic) model, whose symmetry axis is tilted at θ angle to the vertical axis (Fig. 1).

Algorithmic solutions

In case of a thin-layered arrangement, arbitrarily oriented in relation to the Cartesian coordinate system x, y, z, it is appropriate to use general law of tensor rotation in Bond`s formulation [2, 3, 13] which allows to obtain relation between matrix Dφθ _{of elastic modules in} measur-ing coordinate system x, y, z, and appropriate matrix C of these modules in coordinate system x’, y’, z’.

The is the ensuing relation: ( )CR( )T

R

Dϕθ = ϕθ ϕθ (1) where φ denotes the angle of rotation of system x’, y’, z’ with regard to axis z, whereas θ is the tilt angle of the symmetry plane of isotropy TI.

The dimensions of matrices Rφθ and RϕθT are 6 × 6 and they transform the matrix of elastic modules C into sym-metric matrix Dφθ_{. Matrices R}

φθ and RϕθT transpose the

vectors of stress and strain from system x’, y’, z’ to system x, y, z by means of rotational matrix:

θ

θ

ϕ

θ

ϕ

ϕ

θ

ϕ

θ

ϕ

ϕ

θ

cos sin cos sin sin cos cos sin sin cos sin cos A 0 − − − = (2)

Fig. 1. Geometrical model of monoclinal thin-layered system, angle θ is the tilt angle between axis x’ and horizontal plane. Rotation angle φ of system x’, y’, z’

The subject of discussion will be the situation when the rotation angle φ = 90o_{, i.e. when the coincidence of axis} y and y’ occurs (Fig. 1). Then, similarly as for the medium VTI [4] the shear waves of type SH may be separated from longitudinal waves P and shear waves SV. It means that displacements Uy of oscillating particles of the medium

towards axis y are independent from displacements Ux and

Uz in directions x and z respectively. Thus, only

compo-nents Ux and Uz can be discussed, assuming that Uy and

its derivatives equal zero.

Symmetrical matrix C, dimensions 6 × 6, represents components of tensor Cijlk in the medium of transverse

isotropy TI. In abbreviated notation by Voigt [14] this matrix can be presented in this way:

66 44 44 33 13 13 13 11 12 13 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 C C C C C C C C C C C C C = (3) Matrix Dφ=90o_,θ

= D (omitting indices φ and θ) in

dis-cussed case may be presented in the following way:

66 64 55 53 52 51 46 44 35 33 32 31 25 23 22 21 15 13 12 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 d d d d d d d d d d d d d d d d d d d d D = (4)

and the elements can be expressed by means of components of tensor Cijkl (in Voigt`s notation) and the tilt angle θ as

follows:

d11 = C11 cos4 θ + 2C13 cos2 θ sin2 θ + C33 sin4 θ +

+ 4C44 (sin θ cos θ)2 (5)

d21 = d12 = C12 cos2 θ + C13 sin2 θ (6) d31 = d13 = (C11 cos2 θ + C13 sin2 θ) sin2 θ + (C13 cos2 θ +

+ C33 sin2 θ) cos2 θ – 4C44 sin2 θ cos2 θ (7) d51 = d15 = [C13 cos2 θ + C33 sin2 θ – C11 cos2 θ –

– C13 sin2 θ + 2C44 (cos2 θ – sin2 θ –

– sin2 _{θ)] sin θ cos θ (8)}

d22 = C11 (9)

d32 = d23 = C12 sin2 θ + C13 cos2 θ (10) d25 = d52 = (C13 – C12)sin θ cos θ (11) d33 = (C11 sin2 θ + C13 cos2 θ) sin2 θ + (C13 sin2 θ +

+ C33 cos2 θ) cos2 θ + 4C44 sin2 θ cos2 θ (12) d35 = d53 = [C13 sin2 θ +C33 cos2 θ – C11 sin2 θ –

– C13 cos2 θ – 2C44 (cos2 θ – sin2 θ)] sin θ cos θ (13) d44 = C44 cos2 θ + C66 sin2 θ (14) d46 = d64 = (C44 – C66)sin θ cos θ (15) d55 = (C11 – 2C13 + C33) sin2 θ cos2 θ + C44 (cos2 θ – sin2 θ)2 (16) d66 = C44 sin2 θ + C66 cos2 θ (17) and other remaining components are equal zero.

For small tilt angles, when θ → o_{0 we have:} d11 → C11 d12 = d21 → C12 d13 = d31 → C13 d51 = d15 → 0 d23 = d32 → C13 d25 = d52 → 0 d33 → C33 (18) d35 = d53 → 0 d44 → C44 d46 = d64 → 0 d55 → C44 d66 → C66

So for small tilt angles θ → 0, as expected, matrix C

D → →0

θ ₍₁₉₎ The initial point in the discussion on construction of algorithmic solutions will be Hook`s law – basic relation between the tensor of stress Tij and strain Eij, which results

in a conclusion that each stress component is a linear func-tion of strain, i.e.:

Tij = dijkl Ekl = dijkl Elk (20)

While the tensor of strain is

(

l,k kl,

)

lk U U

E = +

2

1 ₍₂₁₎

Tij = dij11E11+dij22E22+dij33E33+2dij23E23+

+2dij13E13+2dij12E12 (22)

Substituting i,j = 1, 2, 3, we obtain all the components of stress tensor: T11, T22, T33, T23 = T32, T13 = T31, T12 = T21. In the matrix notation this is as follows:

0 2 2 0 2 2 0 12 31 23 33 22 11 21 12 13 31 23 32 33 22 11 = = + = = = = = = = = = = = = = = = y , x z , x x , z z , y yz z , z zz y , y yy x , x xx U E U U E U E E U E E U E E U E E D T T T T T T T T T (23)

In relation (23) it has been considered that derivatives of the field with regard to coordinate y equal zero. Start-ing with the general law of movement (disregardStart-ing the external force) 2 2 t U T i j , ij _∂ ∂ =

ρ

(24) where ρ is the medium density, and t denotes time, let us write equations for the horizontal Ux(U1) and vertical Uz(U3) component 21 2 3 13 1 11 t U T T , + , =

ρ

∂_∂ (25a) 23 2 3 33 1 31 t U T T , + , =

ρ

∂_∂ (25b)

Using the matrix equation (23) with relation to equa-tions (25), the following relaequa-tions are received:

+ + + + _{x,zz x,xz z,xx} xx x, d U d U d U U d11 55 2 15 15

(

_{13 55}

)

₅₃ 2 ₂ t U U d U d d x zz z, zx z, + = ∂_∂ + +

ρ

_(26a)

(

+

)

+ + + + _{x,xz x,zz z,zz} xx x, d d U d U d U U d51 31 55 35 33 2 2 55 35 2 t U U d U d z xx z, xz z, _∂ ∂ = + +

ρ

(26b) Adopting for small tilt angles θ d15 = d51 ≈ 0 and d53 = d35 ≈ 0 and applying Fourier`s transformation

(x → kx,z→ kz,t → ω) z for equations (26) we obtain matrix

equation analogical to Christoffel`s relation:

(

)

(

)

2 2 0 55 2 33 55 31 55 13 2 2 55 2 11 ₌ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − + + + − + z x x z z x z x z x U U k d k d k k d d k k d d k d k d

ρω

ρω

₍₂₇₎

In equation (27), kx and kz denote wavenumbers in

horizontal and vertical direction, whereas ω is frequency. It should be noted that matrix equation (27) for medium VTI, i.e. in the case when the tilt angle θ = 0 transforms into analogical equation derived by Q. Han and R.S. Wu [5]. Equation (27) results in dispersion relation

0 2 2 1 4 0k +bk +b = b z z (28) where:

(

)

(

)

(

)

2 2 4 55 11 2 55 11 4 2 55 33 2 55 11 2 13 33 11 2 1 55 33 0 2

ω

ρ

ρω

ρω

+ + − = + − − − = = x x x k d d d d k b d d d d d d d k b d d b (29) In general, equation (28) has four solutions correspond-ing to longitudinal wave qP and shear wave qSV (polarized in plane x-z) forward and backward propagation. Q. Han and R.S. Wu [5] concluded on the basis of numerous experi-ments that wave velocity qSV does not provide significant contribution in the quantity of vertical wavenumber of longitudinal wave qP, therefore it can be assumed that velocity of wave qSV equals zero [1]. On assumption that b0 = 0, we receive:

(

)

2 11 2 4 2 2 33 2 2 13 33 11 2 1 x ' x ' k d b d d d d k b

ρω

ω

ρ

ρω

− = − − = (30) Let expressions _b'

1 and b2' be represented by parameters analogical to Thomsen`s, i.e.:

33 33 11 2d d d − =

ε

(31)

(

) (

)

(

33 55

)

33 2 55 33 2 55 13 2d d d d d d d − − − + = ∂ (32) hence 33 11 2 1 d d q= +

ε

= and for d55 = 0 2 33 2 13 2 1 dd = ∂ +

Assuming that the velocity of longitudinal wave 2 1 33 ⎟⎟ ⎠ ⎞ ⎜⎜ ⎝ ⎛ = ρ C Vp a d33 = C33cos2θ we receive: 2 2 2 33 Vpcos Vpp d =

ρ

θ

=

ρ

(33)

Similarly, using parameters ε and q = 1 + 2ε and η = 2(ε – δ) we have:

(

)

2 2 4 4 2 4 33 33 11d 1 2 d q Vpcos q Vpp d = +

ε

=

ρ

θ

=

ρ

₍₃₄₎

where: Vpp = Vp cos θ denotes the velocity of longitudinal

wave along the vertical axis

(

)

2

(

)

2 4 33 2 13 1 2 d 1 2 Vpp d = + ∂ = + ∂

ρ

₍₃₅₎

(

)

2 11 1 2 Vpp d = +

ε

ρ

(36) From relation (29) and relations (30)-(35) we obtain the expression for vertical component of wavenumber kz

2 1 2 2 2 2 2 2 4 4 ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ − − = x p x p p z S _S qS _k k k

ω

_ω

_η

ω

(37) where: pp p _V

S = 1 denotes vertical slowness.

It is not difficult to notice that parameters q and η which can be seen in expression (37) are represented by their values in model VTI, therefore for horizontally laminated medium

q = qVTI � cos2 θ

η = ηVTI � cos2 θ (38)

where: qVTI and ηVTI denote the values of these parameters

in VTI medium.

Thus relation (37) for medium TTI is analogical to the relation derived by Q. Han and R.S. Wu [5] for the horizontally laminated medium VTI. In order to use it in the process of propagation and migration of compressional waves one should know parameters q and η, the velocity of P waves along the axis and the tilt angle – θ, the plane of isotropy. We will apply received wavenumber kz for

migration MG(F-K) in the domain of wavenumbers and frequencies and spatial coordinates x and time t [8, 11]. Migration process performed in this way occurs in two stages. At the first stage, the relocation of the wave field takes place

(

k ,z z

)

e U

(

k ,z ,ω

)

U' x j+∆ ,

ω

= −ikz0∆z x j (39)

from the level of zj to zj + ∆z by means of exponential

operator with vertical wavenumber kz0 corresponding to a homogeneous medium. At the second stage, correc-tion of the wave field follows U′(x, zj + ∆z, ω) – Fourier

transforms (kx → x) of the field U′(kx, zj + ∆z, ω) by way

of spatial filter Fj(x, ω) = [1 – i/2∆zMj(x)]–1, which is the

sum of Neumann’s power series

( )x z

(

z z

)

ikx x j k k k e dk M 1 2 2 x 0 0 − Σ = − (40)

This relation may be represented in this way:

Uzj + ∆z + = Fj (x, ω)U′(x, zj + ∆z, ω) (41)

The correction positions the wave field in the function of spatial coordinates, taking into account the differences between parameters of a homogenous medium and para-meters of heterogeneous medium in the function of lateral coordinates. For the prestack migration, the algorithm of extrapolation will be the product of corrective functions F related to the sources and receivers, while corrected field U′ will be a function of coordinates of sources and receivers. With the zero-offset migration in relation (37) slowness Sp must be multiplied by 2. Further steps follow

in an analogical way as in model VTI and it was discussed in detail in the article by A. Kostecki [11]. It should be noted that even when the elastic parameters are indepen-dent from the spatial coordinates, it is essential to take into consideration the differences in vertical wavenumbers as a result of different tilt angles of the laminated medium. Artykuł nadesłano do Redakcji 17.09.2009. Przyjęto do druku 29.10.2009.

Recenzent: dr Anna Półchłopek

Literature

[1] Alkhalifah T.: Acoustic approximation for processing in

transversely isotropic media. Geophysics, 63, 623-631,

1998.

[2] Auld B.A.: Acoustic fields and waves in solids. Krieger Publishing Company, vol. 1, 2, 1990.

[3] Bansal R., Sen M.: Finite – difference modelling of S-wave

splitting in anisotropic media. Geophysical Prospecting

56, 293-312, 2008.

[4] Crampin S., Chesnokov E, Hipkin R.: Seismic anisotropy

– the state of the art. Geophys. J. Roy. Astr. Soc., 76, 1-16,

1984.

[5] Han Q., Wu R.S.: A one-way dual domain propagator

for scalar qP waves in VTI medium. Geophysics, vol. 70,

D9-D-17, 2005.

[6] Issack I., Lavton D.: A practical method for estimating

effective parameters of anisotropy from reflection seismic data. Geophysics, vol. 69, 681-689, 2004.

[7] Jianlin Zhu, Ji Dorman: Two-dimensional, three-component

wave propagation in a transversely isotropic medium with arbitrary orientation – finite – element modeling.

Geophy-sics, vol. 65, no 3, pp. 934-942, 2000.

sum-ming up the wavefield in a medium with lateral velocity heterogeneities. Studies by Oil and Gas Institute, no 94,

1998.

[9] Kostecki A., Półchłopek A.: Stable depth extrapolation of

seismic wavefields by Neumann series. Geophysics, 63,

2063-2071, 1998.

[10] Kostecki A., Półchłopek A.: A study of structural

repro-duction in anisotropic medium VTI. Prace INiG, no 137,

pp. 39-60, 2006.

[11] Kostecki A.: Algorithms of depth migration in anisotropic

medium VTI. Nafta-Gaz, no 11, 661-666, 2007.

[12] Postma G.M.: Wave propagation in a stratified media. Geophysics, vol. 20, no 4, 1955.

[13] Schoenberg M.A.: Seismic characterization of reservoirs

containing multiple fracture sets. Geophysical Prospecting,

vol. 57, no 2, pp. 169-186, 2009.

[14] Thomsen L.: Weak elastic anisotropy. Geophysics, vol. 51, 1954-1966, 1986.

Andrzej KOSTECKI – Professor of geophysics. The main subject of interest – electromagnetic and seismic wave propagation, reproduction of deep geological structures by means of seismic migra-tion, the analysis of migration velocities, seismic anisotropy. The author of 130 publications.

ZAKŁAD TECHNIKI STRZELNICZEJ INiG w KRAKowIE

OFERTA NA WYKONYWANIE BADAŃ ODPORNOŚCI CIŚNIENIOWEJ I TERMICZNEJ URZĄDZEŃ

Oferujemy wykonywanie badań odporności ciśnieniowej urządzeń w warunkach podwyższo-nej temperatury i w temperaturze otoczenia na stanowisku termobarycznym INiG w Krakowie.

– maksymalne wymiary gabarytowe:

◦ długość

◦ średnica 1850 mm 140 mm

– maksymalny ciężar: 100 kG

– ciecz robocza: olej Iterm 5Mb

– maksymalne ciśnienie badania: 120 MPa

– maksymalna temperatura badania: 180°C

– ogrzewanie pośrednie w płaszczu olejowym, – rejestracja ciągła ciśnienia i temperatury,

– możliwe wykonywanie badań „do zniszczenia” w przypadku urządzeń o niewielkiej objętości, – możliwość wyprowadzenia sygnału elektrycznego z badanego urządzenia linią 2-przewodową.

Kierownik: mgr inż. Antoni Frodyma

Adres: ul. Bagrowa 1, 30-733 Kraków

Telefon: 12 653-25-12 wew. 137

Faks: 12 653-16-65