Device and method for calculating 3D angle gathers from reverse time migration
View Patent ↗A method for calculating angle domain common image gathers (ADCIGs). The method includes calculating a source wavefield p F of a seismic source; calculating a receiver wavefield p B of a seismic receiver; applying an algorithm of anti-leakage Fourier transform (ALFT) to transform the source wavefield p F to a wavenumber domain; applying the ALFT algorithm to the receiver wavefield to transform the receiver wavefield in the wavenumber domain; determining an imaging condition to the ALFT source and receiver wavefields in the wavenumber domain; computing a reflection angle θ and an azimuth angle φ of the source wavefield p F and receiver wavefield p B in the wavenumber domain; calculating the ADCIGs in the wavenumber domain; and applying an inverse fast Fourier transform (FFT) to determine the ADCIGs in the space domain.
1. A method for calculating, in a computing device, angle domain common image gathers (ADCIGs) for use in generating a final image of a subsurface of the earth, the method comprising:
calculating with wavefield forward propagation, in a time-space domain, a source wavefield p F of a seismic source;
calculating with wavefield back propagation, in the time-space domain, a receiver wavefield p B of a seismic receiver configured to detect a seismic wave generated by the seismic source, wherein the seismic wave emitted by the seismic source enters the subsurface of the earth and gets reflected on a reflector of the subsurface prior to arriving at the seismic receiver;
applying an algorithm of anti-leakage Fourier transform (ALFT) to transform the source wavefield p F to a wavenumber domain;
applying the ALFT algorithm to the receiver wavefield to transform the receiver wavefield in the wavenumber domain;
determining an imaging condition to the ALFT source and receiver wavefields in the wavenumber domain;
computing a reflection angle θ and an azimuth angle φ of the source wavefield p F and receiver wavefield p B in the wavenumber domain;
calculating the ADCIGs in the wavenumber domain; and
applying an inverse fast Fourier transform (FFT) to determine the ADCIGs in the space domain.
2. The method of claim 1 , further comprising:
applying local windows for calculating Fourier coefficients of the ALFT algorithm.
3. The method of claim 2 , further comprising:
using an oversampling scheme in the wavenumber domain.
4. The method of claim 3 , further comprising:
selecting those wavenumbers that correspond to energetic components and removing the remaining wavenumbers prior to the determining step.
5. The method of claim 1 , further comprising:
calculating the source wavefield p F as a forward propagation of a source signature.
6. The method of claim 5 , wherein the source wavefield p F is calculated based on equations:
{
(
1
v
2
∂
2
∂
t
2
-
∇
2
)
p
F
(
x
;
t
)
=
0
,
p
F
(
x
,
y
,
z
=
0
;
t
)
=
δ
(
x
-
s
)
∫
-
∞
t
f
(
t
′
)
ⅆ
t
′
,
where v is the velocity of p F , ∇ 2 is the Laplacian operator, x is a position where the Laplacian operator is applied to the source wavefield p F , t is the time at which the source wavefield p F is present at position x, z is a depth component of the x position, δ(x−s)f(t) is the source signature, δ is a Dirac function, and s is a location of the seismic source.
7. The method of claim 1 , further comprising:
calculating the receiver wavefield p B as a backward propagation of data recorded by the seismic receiver.
8. The method of claim 7 , wherein the receiver wavefield p B is calculated based on equations:
(
1
v
2
∂
2
∂
t
2
-
∇
2
)
p
B
(
x
;
t
)
=
0
,
p
B
(
x
,
y
,
z
=
0
;
t
)
=
δ
G
obs
(
x
,
y
;
s
;
t
)
,
where v is the velocity of p B , ∇ 2 is the Laplacian operator, x is a position where the Laplacian operator is applied to the reflector wavefield p R , t is the time at which the receiver wavefield p B is present at position x, z is a depth of the position x, δG obs (x, y;s;t) is the recorded data, and s is the location of the seismic source.
9. The method of claim 1 , further comprising:
determining a true amplitude migration asymptotic equation in both the reflection angle θ and the azimuth angle φ, wherein the equation is:
R
(
x
,
θ
0
,
φ
0
)
=
∫
∫
∫
δ
(
θ
-
θ
0
)
δ
(
φ
-
φ
0
)
ⅆ
s
ⅆ
r
ⅆ
ω
32
π
ⅈ
ω
cos
θ
A
(
s
,
x
,
r
)
p
z
s
p
z
r
ⅇ
-
ⅈ
ω
T
(
r
,
x
,
s
)
δ
G
obs
sin
θ
,
where R is a reflectivity of the subsurface, x is a position where the reflectivity is calculated, s is a position of the seismic source, r is a position of the seismic receiver, δG obs is data recorded by the seismic receiver, δ is a Dirac function, w is a frequency of the wavefields, A is a product between (i) a ray amplitude from location s to location x and (ii) a ray amplitude from location x to location r, p z s is related to p F and p z r is related to p B .
10. The method of claim 9 , further comprising:
transforming the true amplitude migration asymptotic equation into a true amplitude reverse time migration (RTM) equation which is given by
R
_
(
x
,
θ
0
,
φ
0
)
=
16
π
∫
∫
∫
v
(
x
)
sin
θ
δ
(
θ
-
θ
0
)
δ
(
φ
-
φ
0
)
p
F
p
B
ⅆ
s
ⅆ
ω
,
where the common image gathers (CIGs) are given by R (x,θ 0 ,φ 0 )=R(x,θ 0 ,φ 0 )ν(x)/2 cos θ 0 , p F =G*(x,s)p z s , and p B =∫drG*(x,r)p z r δG obs , with G being a Green function.
11. The method of claim 10 , wherein the ADCIGs are given by:
R
_
s
(
x
,
θ
0
,
φ
0
)
=
16
π
v
(
x
)
sin
θ
0
∑
ω
∑
k
∑
(
k
r
-
k
s
)
δ
(
θ
-
θ
0
)
δ
(
φ
-
φ
0
)
p
F
(
k
s
,
ω
)
p
B
(
k
r
,
ω
)
ⅇ
ⅈ
k
·
x
,
with k being the wavenumber.
12. The method of claim 11 , wherein the ADCIGs are rewritten as
R
_
s
(
x
,
θ
0
,
φ
0
)
=
16
π
v
(
x
)
sin
θ
0
∑
k
I
s
(
k
,
θ
0
,
φ
0
)
ⅇ
ⅈ
k
·
x
,
based on the imaging condition
I
s
(
k
,
θ
0
,
φ
0
)
=
∑
ω
∑
(
k
r
-
k
s
)
δ
(
θ
-
θ
0
)
δ
(
φ
-
φ
0
)
p
F
(
k
s
,
ω
)
p
B
(
k
r
,
ω
)
expressed in the wavenumber domain.
13. The method of claim 12 , further comprising:
applying the inverse FFT transform to I s (k,θ 0 ,φ 0 ) to obtain the ADCIGs.
14. The method of claim 13 , wherein the ADCIGs are 3-dimensional.
15. The method of claim 1 , wherein the computing device is a specific computer that generates the final image and displays the final image on a display.
16. A non-transitory computer readable medium including computer executable instructions, wherein the instructions, when executed, implement a method for calculating angle domain common image gathers (ADCIGs) for use in generating a final image of a subsurface of the earth, the method comprising:
calculating with wavefield forward propagation, in a time-space domain, a source wavefield p F of a seismic source;
calculating with wavefield back propagation, in the time-space domain, a receiver wavefield p B of a seismic receiver configured to detect a seismic wave generated by the seismic source, wherein the seismic wave emitted by the seismic source enters the subsurface of the earth and gets reflected on a reflector of the subsurface prior to arriving at the seismic receiver;
applying an algorithm of anti-leakage Fourier transform (ALFT) to transform the source wavefield p F to a wavenumber domain;
applying the ALFT algorithm to the receiver wavefield to transform the receiver wavefield in the wavenumber domain;
determining an imaging condition to the ALFT source and receiver wavefields in the wavenumber domain;
computing a reflection angle θ and an azimuth angle φ of the source wavefield p F and receiver wavefield p B in the wavenumber domain;
calculating the ADCIGs in the wavenumber domain; and
applying an inverse fast Fourier transform (FFT) to determine the ADCIGs in the space domain.
17. The medium of claim 16 , further comprising:
applying local windows for calculating Fourier coefficients of the ALFT algorithm.
18. The medium of claim 17 , further comprising:
using an oversampling scheme in the wavenumber domain.
19. The medium of claim 18 , further comprising:
selecting those wavenumbers that correspond to energetic components and removing the remaining wavenumbers prior to the determining step.
20. A processing device configured to calculate angle domain common image gathers (ADCIGs) for use in generating a final image of a subsurface of the earth, the processing device comprising:
an interface configured to receive seismic data; and
a processor connected to the interface and configured to,
calculate with wavefield forward propagation, in a time-space domain, a source wavefield p F of a seismic source,
calculate with wavefield back propagation, in the time-space domain, a receiver wavefield p B of a seismic receiver configured to detect a seismic wave generated by the seismic source, wherein the seismic wave emitted by the seismic source enters the subsurface of the earth and gets reflected on a reflector of the subsurface prior to arriving at the seismic receiver,
apply an algorithm of anti-leakage Fourier transform (ALFT) to transform the source wavefield p F to a wavenumber domain,
apply the ALFT algorithm to the receiver wavefield to transform the receiver wavefield in the wavenumber domain,
determine an imaging condition to the ALFT source and receiver wavefields in the wavenumber domain,
compute a reflection angle θ and an azimuth angle φ of the source wavefield p F and receiver wavefield p B in the wavenumber domain,
calculate the ADCIGs in the wavenumber domain, and
apply an inverse fast Fourier transform (FFT) to determine the ADCIGs in the space domain.