Seismic P-wave modelling in an inhomogeneous transversely isotropic medium with a tilted symmetry axis
View Patent ↗An improved method for P-wave modeling in inhomogeneous transversely isotropic media with tilted symmetry axis (TTI media), suitable for anisotropic reverse-time migration, is based on an acoustic TI approximation. The resulting wave equations (2.20) & (2.21) are derived directly from first principles, Hooke's law and the equations of motion, and therefore make no assumptions on spatial variation of medium parameters. Like in the acoustic VTI case, the wave equations are written as a set of two second-order partial differential equations. However, unlike in the acoustic VTI case, the acoustic TTI wave equations contain mixed second-order derivatives. The discretization scheme uses centered finite-difference operators for first- and second-order derivative operators to approximate the mixed and non-mixed second-order derivatives in the wave equation. The discretization scheme is stabilized by slightly weighing down the mixed derivatives, with almost negligible effect on the wave field kinematics.
1. A method for seismic P-wave modelling to generate a seismic image of a subsurface formation that is represented in a Cartesian coordinate frame as an inhomogeneous transversely isotropic (TI) acoustic medium with a tilted symmetry axis of variable non-vertical direction, the method comprising the steps of:
a) measuring seismic P-waves, excited by a seismic source, and propagated through the subsurface formation;
b) generating a pseudo-acoustic stress-strain relationship of a stress tensor and a strain tensor of modelled P-waves, in the Cartesian coordinate frame by means of a TI elastic tensor, wherein each of the tensors is expressed in a rotated local Cartesian coordinate frame which rotated local Cartesian coordinate frame is rotated relative to the Cartesian coordinate frame in which the subsurface formation is imaged and aligned with said tilted symmetry axis of the variable non-vertical direction of the TI medium, in which stress-strain relationship a shear velocity along the axis of symmetry is set to zero, and which stress-strain relationship comprises an axial scalar stress component that is co-axial to the tilted symmetry axis and a lateral scalar stress component in a plane perpendicular to the tilted symmetry axis;
c) combining the pseudo-acoustic stress-strain relationship of step b with an equation of motion to generate a coupled wave equation for the axial and lateral scalar stress components, which equation contains mixed and non-mixed second-order spatial derivatives;
d) discretizing first-order spatial derivatives and non-mixed second-order spatial derivatives by centered finite-differences, with dedicated selection of coefficients;
e) using bi-directional combinations of the discretized first-order derivatives for the mixed second-order derivatives in the coupled wave equation of step c, and using discretized second-order derivatives for the non-mixed second-order derivatives, while stability of an explicit time-stepping method is established by weighing down the mixed second-order derivatives; and
f) forward propagating a simulated shot, and backward propagating seismic P-waves measured in step a), through an anisotropic migration model in accordance with steps a-e to generate the seismic image of the subsurface formation.
2. The method of claim 1 , wherein the elastic TI tensor, the stress tensor, and the strain tensor are expressed in a rotated local Cartesian coordinate system, which is aligned with the variable direction of the axis of symmetry in the TI acoustic medium, and where a pseudo-acoustic approximation is applied by setting the shear velocity along the axis of symmetry to zero, V S =0, thereby generating two scalar wave fields, comprising the axial scalar stress component σ′ V for the axial stress component and the lateral scalar stress component σ′ H for the stress component in the plane perpendicular to the axis of symmetry, and the pseudo-acoustic stress-strain relationship is expressed by the formulas:
σ
H
′
=
ρ
V
P
2
{
(
1
+
2
ɛ
)
(
ɛ
11
′
+
ɛ
22
′
)
+
1
+
2
δ
ɛ
33
′
}
,
σ
V
′
=
ρ
V
P
2
{
1
+
2
δ
(
ɛ
11
′
+
ɛ
22
′
)
+
ɛ
33
′
}
,
where ρ is a medium density, V P is a P-velocity along the axis of symmetry, δ and ε are anisotropy parameters, known as Thomsen parameters, and where ε′ 33 is an axial component of strain along the axis of symmetry, and ε′ 11 and ε′ 22 are lateral strain components in the plane perpendicular to the tilted symmetry axis.
3. The method of claim 2 , wherein the pseudo-acoustic approximation implies a physical stability constraint ε−δ≧0, and the coupled wave equation for the two scalar stress components σ′ V and σ′ H is derived from first principles, which principles comprise the equation of motion and the pseudo-acoustic stress-strain relationship, thereby taking into account spatial variation of medium parameters, including the variable direction of the symmetry axis of the TI medium, and the coupled wave equation is expressed by the formulas:
∂
2
σ
H
′
∂
t
2
=
V
P
2
[
(
1
+
2
ɛ
)
(
ɛ
¨
11
′
+
ɛ
¨
22
′
)
+
1
+
2
δ
ɛ
¨
33
′
]
∂
2
σ
V
′
∂
t
2
=
V
P
2
[
1
+
2
δ
(
ɛ
¨
11
′
+
ɛ
¨
22
′
)
+
ɛ
¨
33
′
]
,
where t is a time-variable, and {umlaut over (ε)}′ 11 , {umlaut over (ε)}′ 22 , {umlaut over (ε)}′ 33 are second-order time-derivatives of ε′ 11 , ε′ 22 , ε′ 33 , which are expressed by
(
ɛ
¨
11
′
+
ɛ
¨
22
′
)
=
∑
k
∑
l
(
R
1
k
R
1
l
+
R
2
k
R
2
l
)
∑
j
∂
2
∂
x
l
∂
x
j
[
(
R
1
k
R
1
j
+
R
2
k
R
2
j
)
σ
H
′
+
R
3
k
R
3
j
σ
V
′
]
,
ɛ
¨
33
′
=
∑
k
∑
l
R
3
k
R
3
l
∑
j
∂
2
∂
x
l
∂
x
j
[
(
R
1
k
R
1
j
+
R
2
k
R
2
j
)
σ
H
′
+
R
3
k
R
3
j
σ
V
′
]
,
wherein the latter two equations specify a spatial differential operator consisting of both mixed second-order derivatives and non-mixed second-order derivatives, in which x l and x j are spatial Cartesian coordinates and R ij are entries of a rotation matrix that transforms vectors from the global to the local rotated coordinate system.
4. The method of claim 3 , further comprising designing:
a discrete first-order derivative operator D and
a discrete second-order derivative operator Δ, which are designed as high-order centered finite-difference operators, with tuned coefficients, according to principles of spectral approximation, in which operator coefficients are chosen with the aim of obtaining good approximations for derivatives of Fourier components, while Δ−D 2 is aimed to be negative definite.
5. The method of claim 4 wherein the second-order temporal derivatives
∂
2
σ
H
′
∂
t
2
and
∂
2
σ
V
′
∂
t
2
are treated by a second-order divided difference, leading to a three-term time-stepping scheme for numerical evolution of the wave fields after discretization of the spatial differential operator.
6. The method of claim 1 wherein a stability analysis is used to show that the evolution of σ′ V and σ′ H is stable for a spatial discretization of the coupled wave equation expressed by the formulas according to claim 3 , in which all second-order derivatives (both mixed and non-mixed) are approximated by combinations D l D j , where D l and D j are discrete first-order derivative operators for the x l - and x j -coordinates, designed according to the method of claim 4 , provided that:
a. ε−δ≧0 throughout the medium,
b. homogeneous boundary conditions are applied,
c. an appropriate ratio of the two scalar stress components is honoured in areas of ellipticity (nodes with ε=δ), and
d. a source term and/or initial condition is consistent with a valid strain tensor.
7. The method of claim 6 wherein to prevent numerical artifacts, which would be generated for the spatial discretization according to claim 6 , in particular for high spatial frequency components of the wave field, discrete second-order derivative operators Δ j , designed in accordance with claim 4 , are used for the non-mixed second-order derivatives with respect to the x j -coordinates, whereupon, to avoid losing stability, each of the differences Δ j −D j 2 is a negative definite operator.
8. The method of claim 7 , wherein the requirement that the differences Δ j −D j 2 are negative definite operators is enforced by weighing down the discrete first-order derivative operators in the remaining mixed derivatives by a factor, which is slightly smaller than unity, thereby retaining stability.
9. The method according to claim 8 , wherein the steps according to claim 8 relieve the specific requirement d. of claim 6 of a strain-consistent source term or initial condition.
10. The method of claim 1 wherein absorbing boundary conditions are implemented by slightly tapering off the wave field at any time step in the proximity of the boundaries of the computational domain.
11. The method of claim 1 wherein the computational grid has uniform grid spacings for each of the x j -coordinates.
12. The method of claim 1 wherein the computational grid has variable grid spacings.
13. The method of claim 3 wherein the second-order temporal derivatives
∂
2
σ
H
′
∂
t
2
and
∂
2
σ
V
′
∂
t
2
are treated by a second-order divided difference, leading to a three-term time-stepping scheme for numerical evolution of the wave fields after discretization of the spatial differential operator.
14. The method of claim 4 wherein to prevent numerical artifacts, which would be generated for the spatial discretization according to claim 6 , in particular for high spatial frequency components of the wave field, discrete second-order derivative operators Δ j , designed in accordance with claim 4 , are used for the non-mixed second-order derivatives with respect to the x j -coordinates, whereupon, to avoid losing stability, each of the differences Δ j −D j 2 is a negative definite operator.
15. The method of claim 5 wherein to prevent numerical artifacts, which would be generated for the spatial discretization according to claim 6 , in particular for high spatial frequency components of the wave field, discrete second-order derivative operators Δ j , designed in accordance with claim 4 , are used for the non-mixed second-order derivatives with respect to the x j -coordinates, whereupon, to avoid losing stability, each of the differences Δ j −D j 2 is a negative definite operator.