Go to the list of seismic processes.      Go to SIOSEIS introduction.
                          PROCESS PSMIGR
                          ------- ------

Document date:  28 July 2000

                        Phase Shift Migration

Reference: Stoffa et al., Split-step Fourier Migration,
                          Geophysics,55,p.410-421,1990.

This process performs a "split-step" Fourier migration on
stacked seismic data.  The method is designed to provide a
fast f-k approach to migration in laterally varying velocity
media.  The approach is straightforward.  The data are migrated
in small depth increments of dz.  For each dz a loop over a
frequency range fmin to fmax is performed in which the data are:

1) transformed from f-x to f-k;
2) phase shifted by exp(i*dz*Kz) [Kz=csqrt(w^2*u0^2-Kx^2)]
   using a reference slowness u0;
3) transformed from f-k to f-x;
4) phase shifted by exp(i*w*(u(x,z)-u0))

Imaging is then done at depth Zj by a sum over frequencies
from fmin to fmax.

Notes:

The complex square root calculation of Kz results in an exponentially
damped response to inhomogeneous interface waves, a difficulty in
other implementations of this type of phase-shift-plus-correction
method.

Extreme lateral velocity variations, such as seen at some
continental margins, may be handled by breaking the migration
into several overlapping panels which are later spliced together.

Big step, ZSKIP:
   You may migrate in one big step through an upper constant
   slowness region such as the water column.  No imaging is done
   for this region.  If there is a deep water delay in your data,
   you may want to remove delay to negate temporal wrap-around.
   This is most easily done through start-end-time pair in process
   DISKIN (i.e. set 0.0 8.0)

Migration bandwidth, FMIN, FMAX:
   Runtime increases in direct proportion to migration bandwidth
   and the number of frequencies migrated.  The number of frequencies
   is initially computed from the FFT.  The larger the FFT, the larger
   the number of frequencies (a 2048 point fft has 2049 frequencies). 
   It is worthwhile doing a couple of small tests to see what
   frequencies are actually useful and adjusting the parameters
   FMIN and FMAX accordingly.

Depth step, DELTAZ:
   The migration depth step dz need not be tiny, but in general small
   enough to provide unaliased sampling of the smallest vertical
   wavenumber of interest.

Data tapers, BPAD, EPAD:
   The user is given the option of padding either side of the data
   panel being migrated. The padded region is filled with copies of
   the end traces which are tapered down within the pad region. This
   is highly recommended as it reduces Gibbs phenomena and Nyquist
   noise.  It should be kept in mind, however, that a power of 2 FFT
   is used and that the length of the FFT is determined as the next
   power of 2 above (data panel length)+BPAD+EPAD.  So if you are
   migrating 1900 traces and you pad with 100 on each side, you will
   be using a 4096 FFT for each depth step. If you had padded with 50
   on each you would be using a 2048 FFT, much faster.

Migration taper, MTAP:
   In addition to the tapered padding of the first and last trace, an
   exponential taper can be applied to the ends of the spatial window
   at each migration step.  This inhibits "wraparound" of migration
   smiles from the data panel ends sides.  An exponental taper on the
   order of MTAP=25 traces is recommended.  The example below should
   clarify the issue.  It is important to note that the FFT length is
   not based on the length of the exponential taper, only on the
   (data panel length)+BPAD+EPAD, so that if care is not taken data
   may be affected by the exponental taper.

   MTAP=5
   |<5>|                            |<5>|
   eeeee1111111111111111111111111111eeeee  <-*1
   tttttttdddddddddddddddddttttttt0000000
   |< 7 >||               ||< 7 >|      |
   |BPAD=7|               |EPAD=7       |
   |      |               |             |
   |      ||             |
   |                                    |
   |<--- 2^n based on 17+(2*7) = 32 --->| <-*2


   *1) this vector is applied to each frequency at each depth
       step.  In f-k jargon this type of taper is called a sponge.
   *2) if EPAD were given as 9 in this case, then the code would
       use an n=6 (2^n=64) would be used for the FFT.


Present limitations include 8192 traces, 5000 depth steps & 800
frequencies.  If the number of depths (nz) times the number of traces 
(nx) exceeds 500,000 use sioseis.BIG, which allows nz*nx = 3,500,000.
Padding should be incluiding when estimating nx.


EXAMPLE MIGRATION:
  Shallow Sediments near ODP HOLE 504B
  diskin
      ipath 504B.dmo.stack.224.624 set 0.0 8.0 end
  end 
  psmigr
     deltax 12.5 deltaz 10 vskip 1500 zskip 3000 ez 10000
     bpad 50 epad 50 mtap 25 twinlen 0.30
     ref 0 nvsmth 3 vpath v.scratch
     fmin 5 fmax 20 path ps.scratch sgypath vmodel.segy
     fno  224 lno  224 vdp 1500.0 3480.0 1850.0 3490.0 1850.0 3721.3 end
     fno  424 lno  424 vdp 1500.0 3442.5 1850.0 3452.5 1850.0 3785.5 end
     fno  524 lno  524 vdp 1500.0 3397.5 1850.0 3407.5 1850.0 3694.3 end
     fno  624 lno  624 vdp 1500.0 3390.0 1850.0 3400.0 1850.0 3686.8 end 
 end
 
PARAMETER DICTIONARY
--------- ----------
EZ     - End Depth, in meters.  The number of samples output will be 
         EDEPTH  / DELTAZ + 1.  The first output sample is ALWAYS 0.
       - Required.  e.g.  ez 6000
DELTAX - The distance between traces, in meters.
         Required.
DELTAZ - The output sample interval in meters per sample.  DELTAZ may
         not exceed 32 since it is carried in the SEGY header in
         millimeters and as a 16 bit integer (thus 32767 millimeters
         is the maximum).  This is analogous to the sample interval
         in time which the SEG-Y format carries as nanoseconds
         (milliseconds / 1000).
         Required.
TWINLEN - Time WINdow LENgth of temporal taper at end of trace. This
          taper is exponential in nature and is given in seconds. The
          time taper should reduce energy migrating from truncation
          of time trace.
          Preset = 0.25
MTAP   - Migration TAPer.  The exponential spatial taper discussed above.
         preset = 25
FMIN   - Minimum frequency of interest.
         Preset = 0     e.g. fmin 5
FMAX   - Maximum frequency of interest.  Run time may be reduced
         significantly by using FMIN/FMAX.  The Nyquist frequency is
         2/(sample interval) or (sample rate)/2, thus 4 mil data has
         has FMAX preset to 125/2 or 62.5.
       - preset = nyquist / 2  or   (4/si)    e.g. fmax 50
VDP    - The interval velocity to use in migration, given as a list of
         Velocity Depth Pairs.  The velocity should be in units of
         meters per second and the depth should be in units of
         meters.
         Required.    velocity range 350 to 32000
NVSMTH - This parameter specifies the number of Velocity SMooTHing
         operations desired before constructing velocity file.
         The NVSMTH parameter is useful for smoothing across velocity
         contrasts which otherwise can cause distortion in depth
         migration.  The running average of NVSMTH velocities is
         calculated after the VDP velocity function is expanded
         into an uniformly sampled function.  NVSMTH applies in
         depth, not spatially.
         Preset = 3
FNO    - The first shot/rp number the parameter list applies to.
         Preset = the first shot/rp received.    e.g.   fno 101
LNO    - The last shot/rp number the parameter list applies to.
         Preset = the last shot/rp received.    e.g.   lno 101
BPAD   - The number of zero amplitude traces to insert prior to the
         first trace.
         Preset = 1   range 1 to 500      e.g. bpad 10
EPAD   - The number of zero amplitude traces to append after the last
         trace.
         Preset = 1   range 1 to 500      e.g. epad 10
ZSKIP  - The first depth to compute.  The first output sample is
         ALWAYS 0.  This saves computer time!  The first velocity of
         VDP of the first FNO/LNO list is used as the velocity for
         time to depth conversion unless parameter VSKIP is given.
         Preset = 0.
VSKIP  - Velocity used in extrapolation for parameter ZSKIP, usually
         the water velocity in marine work. If using ZSKIP, then VSKIP
         is required.
REF    - Reference slowness, u0, from a given depth.  The minimum
         slowness may provide a better reference than the average
         slowness for imaging features such as a rough basement 
         surface beneath low velocity sediments.
       =0, minimum slowness.
       =1, average slowness.
       =2, maximum slowness.
         Preset = 1
PATH   - The pathname (filename) of a scratch file PSMIGR should use
         for the intermediate transposed data.  The purpose of this
         parameter is to allow the user to specify the exact disk
         partition to use in case the "current" partition does not 
         have enough space.
         preset = a scratch file in the current directory
         e.g.    path /user/scratch/moreroom
VPATH  - The pathname (filename) of a file PSMIGR should use for the
         transposed velocity slices.  If the file exists, PSMIGR will
         not calculate a new velocity model; VDP will be ignored.  If
         the file does not exist, PSMIGR will write the velocity
         slices to VPATH so that it may be used in other runs of PSMIGR
         without having to recalculate the velocities.  The velocity
         slices may be quite large and VPATH allows the user to specify
         the exact disk partition to use in case the "current" partition
         does not have enough space.
         preset = a scratch file in the current directory
         e.g.    path /user/scratch/vmoreroom
SGYPATH - The pathname (filename) of an additional SEGY compatible
         velocity file to be output for external purposes. Includes
         the smoothing operators.
         preset = none
         e.g.    sgypath /user/scratch/vmoreroom.segy

 
Copyright (C) 1995 by:
  Woods Hole Oceanographic Institution and
  The Regents of The University of California, 1995
Written by Dan Lizarralde, WHOI, adapted by Paul Henkart, SIO
Modifications by Graham Kent & Dan Lizarralde, WHOI 
ALL RIGHTS RESERVED.
Go to the list of seismic processes.      Go to SIOSEIS introduction.