Inversion, migration, and imaging related to isotropic wave-mode- independent attenuation
A method and apparatus for identifying features of a subsurface region, including: obtaining an initial physical property model and survey data for the subsurface region; identifying a current model to be the initial physical property model; and executing one or more iterations of: generating synthetic data and forward wavefields with the current model and the survey data by forward modeling with forward wave equations representing isotropic wave-mode-independent attenuation; generating adjoint wavefields with the synthetic data and the survey data by adjoint modeling with adjoint wave equations representing isotropic wave-mode-independent attenuation; computing an objective function gradient with the forward wavefields and the adjoint wavefields by solving gradient equations with the corresponding wave equations representing isotropic wave-mode-independent attenuation; computing a search direction of the objective function; searching for a possible improved model along the search direction; and updating the current model to be the possible improved model.
1. A method for identifying features of a subsurface region, comprising:
obtaining an initial physical property model and survey data for the subsurface region;
identifying a current model to be the initial physical property model; and
executing one or more iterations of:
generating synthetic data and forward wavefields with the current model and the survey data by forward modeling with forward wave equations representing isotropic wave-mode-independent attenuation;
generating adjoint wavefields with the synthetic data and the survey data by adjoint modeling with adjoint wave equations representing isotropic wave-mode-independent attenuation;
computing an objective function gradient with the forward wavefields and the adjoint wavefields by solving gradient equations associated with the corresponding wave equations representing isotropic wave-mode-independent attenuation;
computing a search direction of the objective function;
searching for a possible improved model along the search direction;
if the possible improved model more closely matches the survey data than the current model, updating the current model to be the possible improved model; and
if a difference between the synthetic data and the survey data has not converged, executing another iteration; and
after the difference between the synthetic data and the survey data has converged:
identify a final model to be the last current model; and
generate an image of the subsurface region based on features of the final model.
2. The method of claim 1 , further comprising generating an image of the subsurface region with at least one of the current model and the possible improved model.
3. The method of claim 2 , wherein executing one or more iterations comprises at least one of Wave Equation Migration (WEM) and Reverse Time Migration (RTM).
4. The method of claim 1 , wherein at least two iterations are executed.
5. The method of claim 4 , wherein executing one or more iterations comprises at least one of Full Wavefield Inversion (FWI) and Least-Squares Reverse Time Migration (LSRTM).
6. The method of claim 1 , wherein the initial physical property model comprises a property representing arbitrary velocity anisotropy throughout at least a portion of the subsurface region.
7. The method of claim 1 , further comprising, after the difference between the synthetic data and the survey data has converged:
identifying a final model to be the last current model; and
managing hydrocarbons in the subsurface region based on features of the final model.
8. The method of claim 1 , wherein:
the forward wave equations comprise velocity-attenuation forward wave equations,
the adjoint wave equations comprise velocity-attenuation adjoint wave equations, and
the gradient equations comprise velocity-attenuation gradient equations.
9. The method of claim 8 , wherein:
the velocity-attenuation forward wave equations comprise:
ρ
∂
v
→
∂
t
=
Π
σ
→
-
∑
l
=
1
L
ω
l
m
→
l
+
S
→
v
,
∂
σ
→
∂
t
=
C
Π
T
v
→
+
S
→
σ
,
and
∂
m
→
l
∂
t
=
-
ω
l
m
→
l
+
(
a
l
T
Π
T
)
T
σ
→
,
where
{right arrow over (σ)}=(σ xx σ yy σ zz σ yz σ xz σ xy ) T is a stress vector,
{right arrow over (v)}=(v x v y v z ) T is a velocity vector,
C is a stiffness matrix representing the medium properties,
L represents a number of modeled attenuation mechanisms,
π a matrix of spatial derivatives,
p is density,
ω l is an l th relaxation frequency and l=1, 2, . . . , L,
α l is a matrix representing the l th relaxation modulus,
{right arrow over (m)} l =(m x l m y l m z l ) T is a memory variable vector, and
{right arrow over (S)} σ is a source of the stress vector;
the velocity-attenuation adjoint wave equations comprise:
∂
σ
_
→
∂
t
=
Π
T
(
1
ρ
v
_
→
)
-
∑
l
=
1
L
α
l
Π
T
m
_
→
l
+
∂
F
∂
σ
→
,
∂
v
_
→
∂
t
=
Π
C
σ
_
→
+
∂
F
∂
v
→
,
and
∂
m
_
→
l
∂
t
=
-
(
1
ρ
ω
l
v
_
→
+
ω
l
m
_
→
l
)
,
where
{right arrow over (σ)} =( σ xx σ yy σ zz σ yz σ xz σ xy ) T is an adjoint stress vector,
{right arrow over (v)} ( v x v y v z ) T is an adjoint velocity vector,
{right arrow over (m)} l =( m x l m y l m z l ) T is an adjoint memory variable vector, and
F is an objective function; and
the velocity-attenuation gradient equations comprise:
∂
F
∂
C
=
∫
0
T
σ
_
→
(
Π
T
v
→
)
T
dt
,
∂
F
∂
ρ
=
∫
0
T
-
〈
1
ρ
2
v
_
→
,
Π
σ
→
-
∑
l
=
1
L
ω
l
m
→
l
〉
dt
,
and
∂
F
∂
α
l
=
-
∫
0
T
σ
→
(
Π
T
m
_
→
l
)
T
dt
.
10. The method of claim 1 , wherein:
the forward wave equations comprise auxiliary velocity-attenuation forward wave equations,
the adjoint wave equations comprise auxiliary velocity-attenuation adjoint wave equations, and
the gradient equations comprise auxiliary velocity-attenuation gradient equations.
11. The method of claim 10 , wherein:
the auxiliary velocity-attenuation forward wave equations comprise:
∂
v
→
∂
t
=
1
ρ
Π
σ
→
-
η
→
,
∂
σ
→
∂
t
=
C
Π
T
v
→
+
S
→
σ
,
∂
ξ
→
l
∂
t
=
-
η
→
+
ω
l
(
v
→
-
ξ
→
l
)
,
and
η
→
=
∑
l
=
1
L
α
l
ω
l
(
v
→
-
ξ
→
l
)
-
1
ρ
S
→
v
;
where,
{right arrow over (σ)}=(σ xx σ yy σ zz σ yz σ xz σ xy ) T is a stress vector,
{right arrow over (v)}=(v x v y v z ) T is a velocity vector,
C is a stiffness matrix representing the medium properties,
L represents a number of modeled attenuation mechanisms,
π a matrix of spatial derivatives,
p is density,
ω l is an l th relaxation frequency and l=1, 2, . . . , L,
α l is a matrix representing the l th relaxation modulus,
{right arrow over (S)} σ is a source of the stress vector;
ξ
→
l
=
v
→
-
1
ρ
α
l
-
1
m
→
l
.
is an l th auxiliary variable vector, and
the auxiliary velocity-attenuation adjoint wave equations comprise:
∂
σ
_
→
∂
t
=
Π
T
(
1
ρ
v
_
→
)
+
∂
F
∂
σ
→
,
∂
v
_
→
∂
t
=
Π
C
σ
_
→
-
∑
l
=
1
L
ω
l
(
η
_
→
l
-
ξ
_
→
l
)
+
∂
F
∂
v
→
,
∂
ξ
_
→
l
∂
t
=
ω
l
(
η
_
→
l
-
ξ
_
→
l
)
,
and
η
_
→
l
=
α
l
(
v
_
→
+
∑
l
=
1
L
ξ
_
→
l
)
;
and
where
{right arrow over (σ)} =( σ xx σ yy σ zz σ yz σ xz σ xy ) T is an adjoint stress vector,
{right arrow over (v)} ( v x v y v z ) T is an adjoint velocity vector,
{right arrow over (ξ)} l is the l th adjoint auxiliary variable vector, and
F is an objective function for Full Wavefield Inversion (FWI), and
the auxiliary velocity-attenuation gradient equations comprise:
∂
F
∂
C
=
∫
0
T
σ
_
→
(
Π
T
v
→
)
T
dt
,
∂
F
∂
ρ
=
-
∫
0
T
〈
1
ρ
2
v
_
→
,
Π
σ
→
〉
dt
,
and
∂
F
∂
α
l
=
-
∫
0
T
〈
v
_
→
+
∑
l
=
1
L
ξ
_
→
l
,
ω
l
(
v
→
-
ξ
→
l
)
〉
dt
.
12. The method of claim 1 , wherein the forward wave equations are based on an auxiliary variable that is a function of velocity and a memory variable.
13. The method of claim 12 , further comprising using the auxiliary variable to cause spatial derivatives of memory variables to drop out of the adjoint wave equations.
14. The method of claim 1 , further comprising:
defining one or more checkpoints comprising a forward simulation wavefield state; and
utilizing the one or more checkpoints to perform a cross-correlation.
15. The method of claim 1 , further comprising:
defining a plurality checkpoints, each comprising a forward simulation wavefield state; and
utilizing the plurality of checkpoints to access the respective forward simulation wavefield states in time-reverse order.
16. The method of claim 1 , wherein the forward modeling is implemented in the time domain using memory variables.
17. The method of claim 1 , wherein the forward modeling is implemented in the frequency domain.
18. The method of claim 1 , further comprising, after the difference between the synthetic data and the survey data has converged:
identifying an initial attenuation earth model to be the last current model; and
refining the initial attenuation earth model by inversion into a general attenuation earth model.
19. The method of claim 18 , wherein the general attenuation earth model comprises differing attenuation parameters for differing wave propagation modes.
20. The method of claim 18 , wherein the general attenuation earth model comprises anisotropic attenuation parameters for differing wave propagation modes.
21. A system for generating an image of a subsurface region comprising:
a data module comprising survey data for the subsurface region;
a model module comprising an initial physical property model of the subsurface region; and
an equation module comprising wave equations representing isotropic wave-mode-independent attenuation, comprising:
a representation of velocity in the subsurface region comprising a vector that accounts for attenuation effects; and
a representation of stress in the subsurface region comprising a vector that does not account for attenuation effects.
22. The system of claim 21 , wherein the initial physical property model comprises a property representing arbitrary velocity anisotropy throughout at least a portion of the subsurface region.
23. The system of claim 21 , wherein the wave equations comprise:
ρ
∂
v
→
∂
t
=
Π
σ
→
-
∑
l
=
1
L
ω
l
m
→
l
+
S
→
v
,
∂
σ
→
∂
t
=
C
Π
T
v
→
+
S
→
σ
,
∂
m
→
l
∂
t
=
-
ω
l
m
→
l
+
(
a
l
T
Π
T
)
T
σ
→
,
where
{right arrow over (σ)}=(σ xx σ yy σ zz σ yz σ xz σ xy ) T is the representation of stress,
{right arrow over (m)} l =(m x l m y l m z l ) T is the representation of veloity,
{right arrow over (v)}=(v x v y v z ) T is a velocity vector,
C is a stiffness matrix representing the medium properties,
L represents a number of modeled attenuation mechanisms,
π a matrix of spatial derivatives,
p is density,
ω l is an l th relaxation frequency and l=1, 2, . . . , L,
α l is a matrix representing the l th relaxation modulus, and
{right arrow over (S)} σ is a source of the stress.
24. A method for generating an image of a subsurface region, comprising:
obtaining an initial physical property model and survey data for the subsurface region;
identifying a current model to be the initial physical property model;
iteratively executing:
generating synthetic data and forward wavefields with the current model and the survey data by forward modeling with forward wave equations representing isotropic wave-mode-independent attenuation;
generating adjoint wavefields with the synthetic data and the survey data by adjoint modeling with adjoint wave equations representing isotropic wave-mode-independent attenuation;
computing an objective function gradient with the forward wavefields and the adjoint wavefields by solving gradient equations associated with the corresponding wave equations representing isotropic wave-mode-independent attenuation;
computing a search direction of the objective function;
searching for a possible improved model along the search direction;
if the possible improved model more closely matches the survey data than the current model, updating the current model to be the possible improved model; and
if a difference between the synthetic data and the survey data has not converged, executing another iteration; and
after the difference between the synthetic data and the survey data has converged;
identifying a final model to be the last current model; and
generating an image of the subsurface region based on features of the final model.
25. The method of claim 24 , wherein the initial physical property model comprises a property representing arbitrary velocity anisotropy throughout at least a portion of the subsurface region.
26. The method of claim 24 , wherein:
the forward wave equations comprise velocity-attenuation forward wave equations,
the adjoint wave equations comprise velocity-attenuation adjoint wave equations, and
the gradient equations comprise velocity-attenuation gradient equations.
27. A geophysical data analysis system comprising:
a processor; and
a display configured to display graphical representations of a geophysical data set, wherein the geophysical data analysis system is configured to:
obtain an initial physical property model and survey data for the subsurface region;
identify a current model to be the initial physical property model;
iteratively execute:
generate synthetic data and forward wavefields with the current model and the survey data by forward modeling with forward wave equations representing isotropic wave-mode-independent attenuation;
generate adjoint wavefields with the synthetic data and the survey data by adjoint modeling with adjoint wave equations representing isotropic wave-mode-independent attenuation;
compute an objective function gradient with the forward wavefields and the adjoint wavefields by solving gradient equations associated with the corresponding wave equations representing isotropic wave-mode-independent attenuation;
compute a search direction of the objective function;
search for a possible improved model along the search direction;
if the possible improved model more closely matches the survey data than the current model, update the current model to be the possible improved model; and
if a difference between the synthetic data and the survey data has not converged, execute another iteration; and
after the difference between the synthetic data and the survey data has converged:
identify a final model to be the last current model; and
generate an image of the subsurface region based on features of the final model.
28. The system of claim 27 , wherein the initial physical property model comprises a property representing arbitrary velocity anisotropy throughout at least a portion of the subsurface region.
29. The method of claim 27 , wherein:
the forward wave equations comprise velocity-attenuation forward wave equations,
the adjoint wave equations comprise velocity-attenuation adjoint wave equations, and
the gradient equations comprise velocity-attenuation gradient equations.