Skip to content

Commit 6bfbdca

Browse files
authored
Merge pull request #16 from ChevronETC/mbader/slantstack-bugfix
fix bug and simplify slantstack op, and update documentation
2 parents 206b1b5 + 1e32e51 commit 6bfbdca

File tree

4 files changed

+88
-106
lines changed

4 files changed

+88
-106
lines changed

Project.toml

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,6 @@
11
name = "JetPackTransforms"
22
uuid = "e0630bc5-6e1f-509c-a3e4-0aea24bd2b1c"
3-
version = "0.2.0"
3+
version = "0.2.1"
44

55
[deps]
66
FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"

docs/src/slantstack.md

Lines changed: 16 additions & 46 deletions
Original file line numberDiff line numberDiff line change
@@ -1,76 +1,46 @@
11
# JopSlantStack
22

3-
The slant stack is implemented in one of two modes (time or depth). In depth mode, the implemenation is the wave-number domain with a change of variables from the vertical wave-number ``k_z`` and the horizontal offset wave-number ``k_h`` to the ray parameter ``p_h``. In the time mode, the implementation is in the frequency/wave-number domain with a change of variables from frequency ``\omega`` and ``k_h`` to ``p_h``. The change of variables is defined via a mapping from ``k_z`` and ``p`` to ``k_h`` that is derived from dispersion relations for the downward propagating wave
3+
The slant stack is implemented in one of two modes (time or depth). In depth mode, the implementation is the wave-number domain with a change of variables from the vertical wavenumber ``k_z`` and the horizontal subsurface half-offset wavenumber ``k_h`` to the angle of incidence (half of opening angle) ``\theta`` measured with respect to the normal to a planar reflector. An application of this mode would be the conversion of subsurface offset gathers (SSOG) into angle gathers (SSAG). In time mode, the implementation is in the frequency/wavenumber domain with a change of variables from frequency ``\omega`` and surface offset wavenumber ``k_h`` to the ray parameter ``p``. An application of this mode would be the tau-p transform of shot gathers. The change of variables is defined via a mapping from ``k_z`` and ``\theta`` (or ``p``) to ``k_h``.
44

5-
```math
6-
\frac{\omega^2}{c^2} = k_{gx}^2 + k_{gz}^2
7-
```
8-
9-
and the reflected wave
5+
We start by noting the definition of the ray parameter ``p`` with takeoff angle ``\gamma`` (with respect to the vertical)
106

117
```math
12-
\frac{\omega^2}{c^2} = k_{sx}^2 + k_{sz}^2.
8+
p = \frac{\sin\gamma}{c}.
139
```
1410

15-
Next we note the definition of the ray parameter ``p`` with incidence angle ``\theta_s`` and reflected angle ``\theta_g``, and assume that ``\theta=\theta_s=\theta_g``:
11+
Note that the ray parameter is preserved in the system for each wave. The relation between the takoff angle and the plane wave can be derived from the dispersion relation and is given by
1612

1713
```math
18-
p = \frac{1}{c}(\sin\theta_g + \sin\theta_s) = \frac{2}{c}(\sin\theta)
14+
\sin(\gamma) = \frac{k_x}{\omega/c},
1915
```
2016

21-
The relation between the incidence angle and the plane wave is given by:
17+
so that the mapping between ``k_h=k_x`` and ``p`` is
2218

2319
```math
24-
\sin(\theta) = \frac{k_x}{\omega/c}
20+
p = \frac{k_x}{\omega}.
2521
```
2622

27-
So that the mapping between `k_h=2k_x` and `p` is,
28-
23+
In time mode, this provides the mapping between ``k_h`` and ``p`` via the relation
2924
```math
30-
p = 2\frac{k_x}{\omega}
25+
d(\omega,p) = \delta(k_h - p\omega)d(\omega,k_h).
3126
```
3227

33-
In time mode, this provides the mapping between ``k_h`` and ``p`` via the relation:
34-
```math
35-
d(\omega,p) = \delta(k_h - p\omega)d(\omega,k_h)
36-
```s
37-
38-
The offset wave-number is ``k_h=k_{sx}+k_{gx}`` and we make the assumption that ``k_{sx}=k_{gx}`` so that ``k_{sx}=k_{gx}=k_x``, ``k_h=2k_x`` and ``k_{gz}=k_{sz}=k_z``. Now, we can find the mapping from ``p_h`` and ``k_z`` to ``k_h`` starting from the above dispersion relation we find:
28+
In depth mode, we start from the extended imaging condition between source and receiver plane waves
3929

4030
```math
41-
\begin{aligned}
42-
k_h^2 &= \frac{\omega^2}{c^2} - k_z^2 \\
43-
k_h &= \sqrt{\frac{\omega^2}{c^2} - k_z^2} \tag{1}
44-
\end{aligned}
31+
I(x,z,h) = e^{i(k_{sx}(x+h)+k_{sz}z)} e^{-i(k_{rx}(x-h)+k_{rz}z)}=e^{i((k_{sx}-k_{rx})x + (k_{sz}-k_{rz})z + (k_{sx}+k_{rx})h)}.
4532
```
4633

47-
dividing both sides by ``\omega/c``,
34+
The image vertical wavenumber and subsurface half-offset wavenumber at a reflection point are given by ``(k_z,k_h) = (k_{sz}-k_{rz}, k_{sx}+k_{rx})``. Assuming a horizontal reflector, we have ``k_{sx} = k_{rx}``, ``k_{sz} = -k_{rz}``, and ``\gamma_s = 2\pi - \gamma_r = \theta - \pi``. Given that ``\frac{k_{sx}}{k_{sz}} = -\tan(\gamma_s)`` we deduce
4835

4936
```math
50-
\frac{ck_h}{\omega} = \sqrt{1 - \left(\frac{ck_z}{\omega}\right)^2}
37+
k_h = 2k_{sx} = -2k_{sz}\tan(\gamma_s) = -k_z \tan(\theta).
5138
```
5239

53-
note that ``p=k_h/\omega`` so that,
40+
With some trigonometric manipulations, it can be shown that this relation still holds for non-horizontal reflections.
5441

55-
```math
56-
\begin{aligned}
57-
cp &= \sqrt{1 - \left(\frac{ck_z}{\omega}\right)^2} \\
58-
(cp)^2 &= 1 - \left(\frac{ck_z}{\omega}\right)^2 \\
59-
\left(\frac{c}{\omega}\right)^2 &= \frac{1}{k_z^2}\left(1 - (cp)^2\right) \tag{2}
60-
\end{aligned}
61-
```
62-
63-
Substituting equation 2 into equation 1 gives,
64-
65-
```math
66-
\begin{aligned}
67-
k_h^2 &= k_z^2\left(1 - (cp)^2\right)^{-1} - k_z^2 \\
68-
k_h &= k_z\left[\left(1 - (cp)^2\right)^{-1/2} - 1\right] \tag{3}
69-
\end{aligned}
70-
```
42+
Thus, the mapping between ``k_h`` and ``\theta`` is given by
7143

72-
In depth mode, we use equation 3 to map between ``k_h`` and ``p`` via the relation:
7344
```math
74-
d(k_z,p) = \delta(k_h - (k_z\left[\left(1 - (cp)^2\right)^{-1/2} - 1\right]))d(k_z,k_h)
45+
d(k_z,\theta) = \delta(k_h + k_z \tan(\theta))d(k_z,k_h).
7546
```
76-

src/jop_slantstack.jl

Lines changed: 41 additions & 57 deletions
Original file line numberDiff line numberDiff line change
@@ -1,32 +1,33 @@
11
"""
2-
A = JopSlantStack(dom[; dz=10.0, dh=10.0, h0=-1000.0, ...])
2+
A = JopSlantStack(dom[; dz=1.0, dh=1.0, h0=0.0, ...])
33
4-
where `A` is the 2D slant-stack operator mapping for `z-h` to `tau-θ` (depth mode) or `t-h` to `tau-p` (time mode).
4+
where `A` is the 2D slant-stack operator mapping for `z-h` to `z-θ` (depth mode) or `t-h` to `tau-p` (time mode).
55
The domain of the operator is `nz` x `nh` with precision T, `dz` is the depth spacing (or time interval),
6-
`dh` is the half-offset spacing, and `h0` is the origin of the half-offset axis. The additional named optional arguments
6+
`dh` is the offset spacing, and `h0` is the origin of the offset axis. The additional named optional arguments
77
along with their default values are,
88
99
* `mode="depth` - choose between "depth" and "time" to specify if the input domain is `z-h` or `t-h`.
10-
* `theta=-60:1.0:60` - range of opening angles used when `mode="depth"`. The ray parameter is: p=sin(theta/2)/c
11-
* `p=range(-dt/dx,dt/dx,128)` - ray parameter sampling used when `mode="time"`
10+
* `theta=-60:1.0:60` - range of incidence angles used when `mode="depth"`.
11+
* `p=range(-dz/dh,dz/dh,128)` - ray parameter sampling used when `mode="time"`
1212
* `padz=0.0,padh=0.0` - fractional padding in depth and offset to apply before applying the Fourier transform
13-
* `taperz=(0,0)` - beginning and end taper in the z-direction (or t-direction) before transforming from `z-h` to `kz-kh`
14-
* `taperh=(0,0)` - beginning and end taper in the h-direction before transforming from `z-h` to `kz-kh`
15-
* `taperkz=(0,0)` - beginning and end taper in the kz-direction (or frequency) before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
16-
* `taperkh=(0,0)` - beginning and end taper in the kh-direction before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
13+
* `taperz=(0,0)` - beginning and end taper (fractional) in the z-direction (or t-direction) before transforming from `z-h` to `kz-kh`
14+
* `taperh=(0,0)` - beginning and end taper (fractional) in the h-direction before transforming from `z-h` to `kz-kh`
15+
* `taperkz=(0,0)` - beginning and end taper (fractional) in the kz-direction (or frequency) before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
16+
* `taperkh=(0,0)` - beginning and end taper (fractional) in the kh-direction before transforming from `kz-kh` to `z-h` or `f-kh` to `t-h`
1717
1818
# Notes
1919
2020
* It should be possible to extend this operator to 3D
2121
* If the mode is "time", then `padz` is the padding for the time dimension, `dz` is the time sampling interval, `taperz` is the taper for the time dimension and `tapkerkz` is the taper for frequency.
22+
* For mode="depth", typically theta needs to cover both positive and negative angles and must be < 90 degrees.
2223
"""
2324
function JopSlantStack(
2425
dom::JetAbstractSpace{T};
2526
theta = collect(-60.0:1.0:60.0),
2627
p = nothing,
27-
dz = -1.0,
28-
dh = -1.0,
29-
h0 = NaN,
28+
dz = 1.0,
29+
dh = 1.0,
30+
h0 = 0.0,
3031
padz = 0.0,
3132
padh = 0.0,
3233
taperz = (0,0),
@@ -37,7 +38,7 @@ function JopSlantStack(
3738
mode ("depth", "time") || error("expected mode to be either 'depth' or 'time', got mode=$(mode)")
3839
dz < 0.0 && error("expected dz>0.0, got dz=$(dz)")
3940
dh < 0.0 && error("expected dh>0.0, got dh=$(dh)")
40-
isnan(h0) && error("expected finite h0, got h0=$(h0)")
41+
mode == "depth" && any(abs.(theta) .>= 90.0) && error("for mode='depth', all theta values must be < 90 degrees in absolute value")
4142

4243
nz,nh = size(dom)
4344

@@ -55,54 +56,42 @@ function JopSlantStack(
5556
nhfft = nextprod([2,3,5,7], round(Int, nh*(1+padh)))
5657
kh = 2 * fftfreq(nhfft, pi/dh)
5758

58-
# c*p - used for mode=="depth"
59-
cp = @. sin(deg2rad(theta)/2) # factor of 1/2 is to make theta the opening angle -- i.e. go from (theta_s, theta_g)->theta
59+
# tan(theta) - used for mode=="depth"
60+
tant = @. tan(deg2rad(theta))
6061

6162
# p - used for mode=="time"
6263
if p === nothing
6364
p = collect(range(-dz/dh, dz/dh; length=128))
6465
end
6566

6667
# conversions
67-
kz,kh,cp = map(x->convert(Array{T,1}, x), (kz,kh,cp))
68+
kz,kh,tant,p = map(x->convert(Array{T,1}, x), (kz,kh,tant,p))
6869
nzfft,nhfft = map(x->convert(Int64, x), (nzfft,nhfft))
6970
h0 = T(h0)
7071

7172
# tapers
7273
TX = JopTaper(dom, (1,2), (taperz[1],taperh[1]), (taperz[2],taperh[2]))
73-
TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1,mode=="time" ? length(p) : length(cp)), (1,2), (taperkz[1], taperkh[1]), (taperkz[2], taperkh[2]), mode=(:normal,:fftshift))
74+
TK = JopTaper(JetSpace(Complex{eltype(dom)},div(nzfft,2)+1,mode=="time" ? length(p) : length(tant)), (1,2), (taperkz[1], taperkh[1]), (taperkz[2], taperkh[2]), mode=(:normal,:fftshift))
7475

75-
JopLn(dom = dom, rng = JetSpace(T, nz, mode == "time" ? length(p) : length(cp)), df! = JopSlantStack_df!, df′! = JopSlantStack_df′!,
76-
s = (;mode, nzfft, nz_pad_rng, nhfft, kz, kh, cp, p, h0, TX, TK))
76+
JopLn(dom = dom, rng = JetSpace(T, nz, mode == "time" ? length(p) : length(tant)), df! = JopSlantStack_df!, df′! = JopSlantStack_df′!,
77+
s = (;mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK))
7778
end
7879
export JopSlantStack
7980

80-
function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, cp, p, h0, TX, TK, kwargs...) where {T}
81-
nh, np, dh = size(m,2), mode == "time" ? length(p) : length(cp), abs(kh[2]-kh[1])
81+
function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
82+
nh, np, dh = size(m,2), mode == "time" ? length(p) : length(tant), abs(kh[2]-kh[1])
8283

8384
mpad = zeros(T, nzfft, nhfft)
8485
mpad[nz_pad_rng,1:nh] = TX*m
8586

8687
M = rfft(mpad)
8788

88-
if mode == "depth"
89-
# adjoint of conjugate symmetry in kh
90-
for ikz = 1:div(nzfft,2)+1
91-
for ikh = 2:div(nhfft,2)
92-
M[ikz,ikh] += conj(M[ikz,nhfft-ikh+2])
93-
end
94-
for ikh = div(nhfft,2)+2:nhfft
95-
M[ikz,ikh] = 0
96-
end
97-
end
98-
end
99-
10089
compute_kh = mode == "depth" ? slantstack_compute_kh_from_kz : slantstack_compute_kh_from_frequency
101-
is_out_of_bounds = mode == "depth" ? (ikh_m1, ikh_p1) -> (ikh_m1 < 1 || ikh_p1 > div(nhfft,2)+1) : (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
90+
is_out_of_bounds = (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
10291

10392
D = zeros(eltype(M), size(M,1), np)
10493
for ikz = 1:div(nzfft,2)+1, ip = 1:np
105-
ikh_m1, ikh_p1, _kh = compute_kh(ikz, ip, cp, p, kz, kh, nhfft)
94+
ikh_m1, ikh_p1, _kh = compute_kh(ikz, ip, tant, p, kz, kh, nhfft)
10695

10796
is_out_of_bounds(ikh_m1, ikh_p1) && continue
10897

@@ -111,11 +100,11 @@ function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, n
111100
continue
112101
end
113102

114-
d_p1 = M[ikz,ikh_p1]*exp(-im*kh[ikh_p1]*h0)
115-
a_p1 = abs(kh[ikh_p1] - _kh)/dh
103+
a_m1 = (kh[ikh_p1] - _kh)/dh
104+
a_p1 = 1 - a_m1
116105

106+
d_p1 = M[ikz,ikh_p1]*exp(-im*kh[ikh_p1]*h0)
117107
d_m1 = M[ikz,ikh_m1]*exp(-im*kh[ikh_m1]*h0)
118-
a_m1 = abs(_kh - kh[ikh_m1])/dh
119108

120109
D[ikz,ip] = a_m1*d_m1 + a_p1*d_p1
121110
end
@@ -124,20 +113,20 @@ function JopSlantStack_df!(d::AbstractArray{T,2}, m::AbstractArray{T,2}; mode, n
124113
d .= _d[nz_pad_rng,1:np] ./ nzfft
125114
end
126115

127-
function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, cp, p, h0, TX, TK, kwargs...) where {T}
128-
nh, np, dh = size(m,2), mode == "time" ? length(p) : length(cp), abs(kh[2]-kh[1])
116+
function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode, nzfft, nz_pad_rng, nhfft, kz, kh, tant, p, h0, TX, TK, kwargs...) where {T}
117+
nh, np, dh = size(m,2), mode == "time" ? length(p) : length(tant), abs(kh[2]-kh[1])
129118

130119
dpad = zeros(T, nzfft, np)
131120
dpad[nz_pad_rng,:] = d
132121

133122
D = TK * (rfft(dpad, 1) ./ nzfft)
134123

135124
compute_kh = mode == "depth" ? slantstack_compute_kh_from_kz : slantstack_compute_kh_from_frequency
136-
is_out_of_bounds = mode == "depth" ? (ikh_m1, ikh_p1) -> (ikh_m1 < 1 || ikh_p1 > div(nhfft,2)+1) : (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
125+
is_out_of_bounds = (ikh_m1, ikh_p1)->(ikh_m1 < 1 || ikh_p1 > nhfft)
137126

138127
M = zeros(Complex{T}, div(nzfft,2)+1, nhfft)
139128
for ikz = 1:div(nzfft,2)+1, ip = 1:np
140-
ikh_m1, ikh_p1, _kh = compute_kh(ikz, ip, cp, p, kz, kh, nhfft)
129+
ikh_m1, ikh_p1, _kh = compute_kh(ikz, ip, tant, p, kz, kh, nhfft)
141130

142131
is_out_of_bounds(ikh_m1, ikh_p1) && continue
143132

@@ -146,38 +135,33 @@ function JopSlantStack_df′!(m::AbstractArray{T,2}, d::AbstractArray{T,2}; mode
146135
continue
147136
end
148137

138+
a_m1 = (kh[ikh_p1] - _kh)/dh
139+
a_p1 = 1 - a_m1
140+
149141
m_p1 = D[ikz,ip]*exp(im*kh[ikh_p1]*h0)
150-
a_p1 = (kh[ikh_p1] - _kh)/dh
151142
M[ikz,ikh_p1] += a_p1*m_p1
152143

153144
m_m1 = D[ikz,ip]*exp(im*kh[ikh_m1]*h0)
154-
a_m1 = (_kh - kh[ikh_m1])/dh
155145
M[ikz,ikh_m1] += a_m1*m_m1
156146
end
157147

158-
if mode == "depth"
159-
# conjugate symmetry in kh
160-
for ikz = 1:div(nzfft,2)+1, ikh = 2:div(nhfft,2)
161-
M[ikz,nhfft-ikh+2] = conj(M[ikz,ikh])
162-
end
163-
end
164-
165148
m .= TX * (brfft(M, nzfft)[nz_pad_rng,1:nh])
166149
end
167150

168-
@inline function slantstack_compute_kh_from_kz(ikz::Int64, ip::Int64, cp, p, kz, kh, nhfft)
169-
_kh = 1/(1-cp[ip]^2) - 1
170-
_kh < 0 && return -1,-1,1.0
171-
_kh = 2*kz[ikz]*sqrt(_kh)
151+
@inline function slantstack_compute_kh_from_kz(ikz::Int64, ip::Int64, tant, p, kz, kh, nhfft)
152+
_kh = -kz[ikz]*tant[ip]
172153

173154
ikh_m1 = floor(Int64, _kh/kh[2]) + 1
174155
ikh_p1 = ceil(Int64, _kh/kh[2]) + 1
175156

157+
ikh_m1 = ikh_m1 < 1 ? nhfft + ikh_m1 : ikh_m1
158+
ikh_p1 = ikh_p1 < 1 ? nhfft + ikh_p1 : ikh_p1
159+
176160
ikh_m1, ikh_p1, _kh
177161
end
178162

179-
@inline function slantstack_compute_kh_from_frequency(iω, ip, cp, p, ω, kh, nhfft)
180-
_kh = p[ip]*ω[iω] / 2
163+
@inline function slantstack_compute_kh_from_frequency(iω::Int64, ip::Int64, tant, p, ω, kh, nhfft)
164+
_kh = p[ip]*ω[iω]
181165

182166
ikh_m1 = floor(Int64, _kh/kh[2]) + 1
183167
ikh_p1 = ceil(Int64, _kh/kh[2]) + 1

test/jop_slantstack.jl

Lines changed: 30 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
using JetPackTransforms, Jets, Test
1+
using JetPackTransforms, Jets, Test, LinearAlgebra
22

33
# depth mode
44
@testset "JetSlantStack, dot product test" for T in (Float32, Float64)
@@ -22,7 +22,7 @@ end
2222
d = A*m
2323
v,i = findmax(d)
2424
@test i[1] == 32
25-
@test i[2] == findfirst(x->x0, state(A).cp)
25+
@test i[2] == findfirst(x->x0, state(A).tant)
2626
end
2727

2828
# time mode
@@ -48,3 +48,31 @@ end
4848
v,i = findmax(d)
4949
@test i[1] == 32
5050
end
51+
52+
# depth-time parity test
53+
@testset "JetSlantStack, parity" begin
54+
theta = collect(-45:0.5:45)
55+
# the ray parameters that would give the same results as theta is
56+
# p = -tan(θ)
57+
p = @. -tan(deg2rad(theta))
58+
59+
Ad = JopSlantStack(JetSpace(Float64, 128, 129); mode="depth", theta=theta, p=p, dz=0.01, dh=0.01, h0=-0.64, padz=1, padh=1, taperkz = (0.5,0.5), taperkh = (0.5,0.5))
60+
At = JopSlantStack(JetSpace(Float64, 128, 129); mode="time" , theta=theta, p=p, dz=0.01, dh=0.01, h0=-0.64, padz=1, padh=1, taperkz = (0.5,0.5), taperkh = (0.5,0.5))
61+
m = zeros(domain(Ad))
62+
for i = 1:129
63+
m[div(i,4)+20,i] = 1.0
64+
end
65+
m[16,64] = 1
66+
m[64,96] = 1
67+
dd = Ad*m
68+
dt = At*m
69+
70+
md = Ad'*dd
71+
mt = At'*dt
72+
73+
err_d = maximum(abs.(dd .- dt)) / maximum(abs.(dd))
74+
err_m = maximum(abs.(md .- mt)) / maximum(abs.(md))
75+
@show err_d, err_m
76+
@test err_d < 1e-7
77+
@test err_m < 1e-7
78+
end

0 commit comments

Comments
 (0)