Next: About this document ...
Up: Fomel: Velocity continuation
Previous: INTEGRAL VELOCITY CONTINUATION AND
The differential equation (55), with the first-order
derivative term neglected, has the mathematical form almost similar to
the 15-degree wave extrapolation equation. Its finite-difference
implementation, first described by Claerbout
1986 and Li 1986, is
analogous to that of the 15-degree equation
Claerbout (1976), except for the variant coefficients. We
can write the implicit unconditionally stable finite-difference scheme
for the velocity continuation equation in the form
| ![\begin{displaymath}
({\bf I} + a_{j+1}^{i+1}\,{\bf T})\,{\bf P}_{j+1}^{i+1} -
(...
...}^{i} +
({\bf I} + a_{j}^{i}\,{\bf T})\,{\bf P}_{j}^{i} = 0\;,\end{displaymath}](img139.gif) |
(95) |
where index i corresponds to the time dimension, index j
corresponds to the velocity dimension,
is a vector along the
midpoint direction,
is the identity matrix,
represents the second-derivative operator in midpoint,
and the coefficient a has the expression
| ![\begin{displaymath}
a_{j}^{i} = v_j\,t_i\,{{\Delta v\,\Delta t}\over {(\Delta x)^2}}\;.\end{displaymath}](img143.gif) |
(96) |
In the two-dimensional case, equation (95) reduces to a
tridiagonal system of linear equations, which can be easily inverted.
The direction of stable propagation is either forward in velocity and
backward in time or backward in velocity and forward in time as shown
in Figure 12.
vlcfds
Figure 12 Finite-difference scheme for the
velocity continuation equation. A stable propagation is either forward
in velocity and backward in time (left plot) or backward in velocity
and forward in time (right plot).
|
| ![vlcfds](../Gif/vlcfds.gif) |
The following simple ratfor program implements the finite-difference
velocity continuation. It is slightly modified from Jon Claerbout's
original version and has a more straightforward loop structure than
Zhiming Li's Caso15 program.
module velconx {
use diffxx2
use tridiag
real, dimension (nx), allocatable :: r, ld, rd, rhs
real, dimension (:,:), pointer, private :: pt
real :: v1, v2, t0, dt
real, private :: dv, a0
real, parameter, private :: b = 0.122996
integer :: nt, nx, nv, inv
#% _init (inv, v1,v2,t0, nt,nx,nv, dt,dx)
real, intent (in) :: dx
allocate (pt(nt,nx))
dv = (v1-v2)/nv; a0=(4.*dx*dx)/(dt*dv)
call rtris_init (nx)
#% _lop (p1 (nt,nx), p2 (nt,nx))
integer :: iv1, iv2, ivs, iv
integer :: it1, it2, its, it
real :: v, t, a, b1, b2, offd, diag
real, dimension (:), pointer :: l
if (adj) {
iv1=1; iv2=nv; ivs= 1
it1=nt; it2= 2; its=-1; pt = p2
} else {
iv1=nv; iv2= 1; ivs=-1
it1= 2; it2=nt; its= 1; pt = p1
}
do iv=iv1,iv2,ivs { v=sqrt(v2+(iv-0.5)*dv); rhs = 0.
do it=it1,it2,its { t=t0+dt*(it-1)
select case (inv) {
case (0)
t=sqrt(t)*v; a=a0/t ! Pseudounitary
case(1)
t=t*v; a=a0/v ! Claerbout's
case default
a=a0/(t*v); t=v ! True-amplitude
}
b1=b*a+t; b2=b*a-t; offd=a-b2; diag=a-2*b2
l => pt(it,:); call diffxx(l,ld)
rhs = rhs + a*l + b1*ld
call tris(offd,b2,diag,b2,offd,rhs,r); call diffxx(r,rd)
rhs = a*r + b1*rd - a*l - b2*ld ; l = r
}}
if (adj)
p1 = p1 + pt
else
p2 = p2 + pt
#% _close
deallocate (pt)
call rtris_close ()
}
The parameter adj controls the direction of propagation. It is
equal to zero for backward propagation, which corresponds to the
modeling (demigration) operator. The parameter inv controls the
amplitude behavior by introducing time-dependent divisors to equation
(95). For inv=1, the program implements
Claerbout's velocity continuation equation. For inv=2, it
implements the modified equation (55). The value of
inv=0 corresponds to the intermediate case. It leads to the
pseudounitary velocity continuation, for which the reverse
continuation is the exact adjoint operator
Fomel (1996). We can easily test different types of
amplitude behavior with the dot-product test and its modifications.
The parameter b is required for the ``1/6 trick'' introduced by
Claerbout 1985 to increase the accuracy
of the second-derivative operator
in (95). The
second-order difference in subroutine diffxx implies simple
zero-slope boundary conditions on the midpoint coordinate. The call to
rtris solves the triadiagonal system.
In order to test the performance of the finite-difference velocity
continuation method, I use a simple synthetic model from Basic
Earth Imaging Claerbout (1995). The ``reflectivity'' model is
shown in Figure 13. It contains several features challenging
the migration performance: dipping beds, unconformity, syncline,
anticline, and fault.
vlcmod
Figure 13 Synthetic model for testing
finite-difference migration by velocity continuation.
|
| ![vlcmod](../Gif/vlcmod.gif) |
The following series of figures compares invertability of different
migration methods. In all cases, constant-velocity modeling
(demigration) by the adjoint operator was followed by migration with
the correct velocity v=1500 m/sec. Figures 14,
15, 16, 17 show the results of modeling
and migration with fast (nearest-neighbor) Kirchhoff, antialiased
Kirchhoff Fomel and Biondi (1995), phase-shift (Gazdag), and Stolt methods
respectively. These figures should be compared with Figure
18, showing the analogous result of the finite-difference
velocity continuation. The comparison reveals a remarkable
invertability of velocity continuation, which reconstructs accurately
the main features and frequency content of the model. Since the
forward operators were different for different migrations, this
comparison did not test the migration properties themselves. For such
a test, I compare the results of the phase-shift and
velocity-continuation migrations after Stolt modeling. The result of
velocity continuation, shown in Figure 19 is at least as
accurate as that of the phase-shift method.
vlckir
Figure 14 Result of modeling and
migration with the fast Kirchhoff method. Left plot shows the
reconstructed model. Right plot compares the average amplitude
spectrum of the true model with that of the reconstructed image.
vlckaa
Figure 15 Result of modeling and
migration with the antialiased Kirchhoff method. Left plot shows the
reconstructed model. Right plot compares the average amplitude
spectrum of the true model with that of the reconstructed image.
vlcpha
Figure 16 Result of modeling and
migration with the phase-shift method. Left plot shows the
reconstructed model. Right plot compares the average amplitude
spectrum of the true model with that of the reconstructed image.
vlcsto
Figure 17 Result of modeling and
migration with Stolt method. Left plot shows the reconstructed
model. Right plot compares the average amplitude spectrum of the true
model with that of the reconstructed image.
vlcvel
Figure 18 Result of modeling and
migration with the finite-difference velocity continuation. Left plot
shows the reconstructed model. Right plot compares the average
amplitude spectrum of the true model with that of the reconstructed
image.
vlcspv
Figure 19 Top plot:
modeling with Stolt method, migration with the phase-shift
method. Bottom plot: modeling with Stolt method, migration with the
finite-difference velocity continuation. Left plots show the
reconstructed models. Right plots compare the average amplitude
spectrum of the true model with that of the reconstructed image.
We can conclude that finite-difference velocity continuation is an
attractive migration method. It possesses remarkable invertability
properties, which might be useful in some applications. According to
Li 1986, the computational speed of this method
compares favorably with that of Stolt migration. The advantage is
apparent for cascaded migration or migration with multiple velocity
models. In these cases, the cost of Stolt migration increases in
direct proportion with the number of velocity models, while the cost
of velocity continuation stays the same.
Next: About this document ...
Up: Fomel: Velocity continuation
Previous: INTEGRAL VELOCITY CONTINUATION AND
Stanford Exploration Project
9/12/2000