NMES Margrave
NMES Margrave
NMES Margrave
of
Exploration Seismology
with algorithms in MATLAB
Gary F. Margrave
The most important thing to know about this draft is that it is unfinished. This means that it
must be expected to have rough passages, incomplete sections, missing references, and references
to nonexistent material. Despite these shortcomings, I hope this document will be welcomed for
its description of the CREWES MATLAB software and its discussion of velocity, raytracing and
migration algorithms. I only ask that the reader be tolerant of the rough state of these notes.
Suggestions for the improvement of the present material or for the inclusion of other subjects are
most welcome.
Exploration seismology is a complex technology that blends advanced physics, mathematics
and computation. Often, the computational aspect is neglected in teaching because, traditionally,
seismic processing software is part of an expensive and complex system. In contrast, this book is
structured around a set of computer algorithms that run effectively on a small personal computer.
These algorithms are written in MATLAB code and so the reader should find access to a MATLAB
installation. This is not the roadblock that it might seem because MATLAB is rapidly gaining
popularity and the student version costs little more than a typical textbook.
The algorithms are grouped into a small number of toolboxes that effectively extend the function-
ality of MATLAB. This allows the student to experiment with the algorithms as part of the process
of study. For those who only wish to gain a conceptual overview of the subject, this may not be an
advantage and they should probably seek a more appropriate book. On the other hand, those who
wish to master the subject and perhaps extend it through the development of new methods, are my
intended audience. The study of these algorithms, including running them on actual data, greatly
enriches learning.
The writing of this book has been on my mind for many years though it has only become physical
relatively recently. The material is the accumulation of many years of experience in both the hydro-
carbon exploration industry and the academic world. The subject matter encompasses the breadth
of exploration seismology but, in detail, reflects my personal biases. In this preliminary edition,
the subject matter includes raytracing, elementary migration, some aspects of wave-equation mod-
elling, and velocity manipulation. Eventually, it will also include theoretical seismograms, wavelets,
amplitude adjustment, deconvolution, filtering (1D and 2D), statics adjustment, normal moveout
removal, stacking and more. Most of the codes for these purposes already exists and have been used
in research and teaching at the University of Calgary since 1995.
During the past year, Larry Lines and Sven Treitel have been a constant source of encouragement
and assistance. Pat Daley’s guidance with the raytracing algorithms has been most helpful. Pat,
Larry, Sven, Dave Henley, and Peter Manning have graciously served as editors and test readers.
Many students and staff with CREWES have stress-tested the MATLAB codes. Henry Bland’s
technical computer support has been essential. Rob Stewart and Don Lawton have been there with
moral support on many occasions.
I thank you for taking the time to examine this document and I hope that you find it rewarding.
ii
Contents
Preface ii
1 Introduction 1
1.1 Scope and Prerequisites . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.1.1 Why MATLAB? . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.1.2 Legal matters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.2 MATLAB conventions used in this book . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.3 Dynamic range and seismic data display . . . . . . . . . . . . . . . . . . . . . . . . . 6
1.3.1 Single trace plotting and dynamic range . . . . . . . . . . . . . . . . . . . . . 6
1.3.2 Multichannel seismic display . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
1.3.3 The plotimage picking facility . . . . . . . . . . . . . . . . . . . . . . . . . . 16
1.3.4 Drawing on top of seismic data . . . . . . . . . . . . . . . . . . . . . . . . . . 17
1.4 Programming tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.1 Scripts . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
1.4.2 Functions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
The structure of a function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
1.4.3 Coping with errors and the MATLAB debugger . . . . . . . . . . . . . . . . . 21
1.5 Programming for efficiency . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.5.1 Vector addressing . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1.5.2 Vector programming . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
1.5.3 The COLON symbol in MATLAB . . . . . . . . . . . . . . . . . . . . . . . . 26
1.5.4 Special values: NaN, Inf, and eps . . . . . . . . . . . . . . . . . . . . . . . . . 27
1.6 Chapter summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28
2 Velocity 29
2.1 Instantaneous velocity: vins or just v . . . . . . . . . . . . . . . . . . . . . . . . . . . 30
2.2 Vertical traveltime: τ . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31
2.3 vins as a function of vertical traveltime: vins (τ ) . . . . . . . . . . . . . . . . . . . . . 31
2.4 Average velocity: vave . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
2.5 Mean velocity: vmean . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
2.6 RMS velocity: vrms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 34
2.7 Interval velocity: vint . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
2.8 MATLAB velocity tools . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
2.9 Apparent velocity: vx , vy , vz . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
2.10 Snell’s Law . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
2.11 Raytracing in a v(z) medium . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
iii
iv CONTENTS
3 Wave propagation 65
3.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
3.2 The wave equation derived from physics . . . . . . . . . . . . . . . . . . . . . . . . . 66
3.2.1 A vibrating string . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
3.2.2 An inhomogeneous fluid . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
3.3 Finite difference modelling with the acoustic wave equation . . . . . . . . . . . . . . 71
Introduction
1 Migration refers to the fundamental step in creating an earth image from scattered data.
1
2 CHAPTER 1. INTRODUCTION
with many MATLAB applications, it helps greatly if the reader has had some experience with
linear algebra. Hopefully, the concepts of matrix, row vector, column vector, and systems of linear
equations will be familiar.
Traditional languages like C and Fortran originated in an era when computers were room-sized
behemoths and resources were scarce. As a result, these languages are oriented towards simplifying
the task of the computer at the expense of the human programmer. Their cryptic syntax leads to
efficiencies in memory allocation and computation speed that were essential at the time. However,
times have changed and computers are relatively plentiful, powerful, and cheap. It now makes sense
to shift more of the burden to the computer to free the human to work at a higher level. Spending
an extra $100 to buy more RAM may be more sensible than developing a complex data handling
scheme to fit a problem into less space. In this sense, MATLAB is a higher-level language that frees
the programmer from technical details to allow time to concentrate on the real problem.
Of course, there are always people who see these choices differently. Those in disagreement with
the reasons cited here for MATLAB can perhaps take some comfort in the fact that MATLAB syntax
is fairly similar to C or Fortran and translation is not difficult. Also, The MathWorks markets a
MATLAB ”compiler” that emits C code that can be run through a C compiler.
0.1
0.05
−0.05
−0.1
0 0.02 0.04 0.06 0.08 0.1 0.12 0.14 0.16 0.18 0.2
time (sec)
1
linear scale
0.5
0
0 50 100 150 200 250
Hz
0
−20
dd down
−40
−60
0 50 100 150 200 250
Hz
3 Amp=abs(Wavem);
4 dbAmp=20*log10(Amp/max(Amp));
5 subplot(3,1,1);plot(t,wavem);xlabel(’time (sec)’)
6 subplot(3,1,2);plot(f,abs(Amp));xlabel(’Hz’);ylabel(’linear scale’)
7 subplot(3,1,3);plot(f,dbAmp);xlabel(’Hz’);ylabel(’db down’)
End Code
The actual MATLAB code is displayed in an upright typewriter font while introductory remarks
are emphasized like this. The Code Snippet does not employ typographic variations to indicate
which functions are contained in the NumMethToolbox as is done in the text proper.
It has proven impractical to discuss all of the input parameters for all of the programs shown
in Code Snippets. Those germane to the current topic are discussed, but the remainder are left for
the reader to explore using MATLAB ’s interactive help facility. For example, in Code Snippet 1.1
wavemin creates a minimum phase wavelet sampled at .002 seconds, with a dominant frequency of
20 Hz, and a length of .2 seconds. Then fftrl computes the Fourier spectrum of the wavelet and
abs constructs the amplitude spectrum from the complex Fourier spectrum. Finally, the amplitude
spectrum is converted to decibels (db) using the formula
Amp(f )
Ampdecibels (f ) = 20 log . (1.1)
max(Amp(f ))
To find out more about any of these functions and their inputs and outputs, type, for example, “help
fftrl” at the MATLAB prompt.
1.2. MATLAB CONVENTIONS USED IN THIS BOOK 5
This example illustrates several additional conventions. Seismic traces2 are discrete time series
but the textbook convention of assuming a sample interval of unity in arbitrary units is not appro-
priate for practical problems. Instead, two vectors will be used to describe a trace, one to give its
amplitude and the other to give the temporal coordinates for the first. Thus wavemin returns two
vectors (that they are vectors is not obvious from the Code Snippet) with wavem being the wavelet
amplitudes and t being the time coordinate vector for wavem. Thus the upper part of Figure 1.1 is
created by simply cross plotting these two vectors: plot(t,wavem). Similarly, the Fourier spectrum,
Wavem, has a frequency coordinate vector f. Temporal values must always be specified in seconds
and will be returned in seconds and frequency values must always be in Hz. (One Hz is one cycle per
second). Milliseconds or radians/second specifications will never occur in code though both cyclical
frequency f and angular frequency ω may appear in formulae (ω = 2πf ).
Seismic traces will always be represented by column vectors whether in the time or Fourier
domains. This allows an easy transition to 2-D trace gathers, such as source records and stacked
sections, where each trace occupies one column of a matrix. This column vector preference for signals
can lead to a class of simple MATLAB errors for the unwary user. For example, suppose a seismic
trace s is to be windowed to emphasize its behavior in one zone and de-emphasize it elsewhere.
The simplest way to do this is to create a window vector win that is the same length as s and use
MATLAB’s .* (dot-star) operator to multiply each sample on the trace by the corresponding sample
of the window. The temptation is to write a code something like this:
Code Snippet 1.2.2. This code creates a synthetic seismogram using the convolutional model but
then generates an error while trying to apply a (triangular) window.
1 [r,t]=reflec(1,.002,.2); %make a synthetic reflectivity
2 [w,tw]=wavemin(.002,20,.2); %make a wavelet
3 s=convm(r,w); % make a synthetic seismic trace
4 n2=round(length(s)/2);
5 win=[linspace(0,1,n2) linspace(1,0,length(s)-n2)]; % a triangular window
6 swin=s.*win; % apply the window to the trace
End Code
MATLAB’s response to this code is the error message:
??? Error using ==> .* Matrix dimensions must agree.
On line 6 ==> swin=s.*win;
The error occurs because the linspace function generates row vectors and so win is of size 1x501
while s is 501x1. The .* operator requires both operands to be have exactly the same geometry.
The simplest fix is to write swin=s.*win(:); which exploits the MATLAB feature that a(:) is
reorganized into a column vector (see page 26 for a discussion) regardless of the actual size of a .
As mentioned previously, two dimensional seismic data gathers will be stored in ordinary matri-
ces. Each column is a single trace, and so each row is a time slice. A complete specification of such
a gather requires both a time coordinate vector, t, and a space coordinate vector, x.
Rarely will the entire code from a function such as wavemin or reflec be presented. This is
because the code listings of these functions can span many pages and contain much material that
is outside the scope of this book. For example, there are often many lines of code that check input
parameters and assign defaults. These tasks have nothing to do with numerical algorithms and so
will not be presented or discussed. Of course, the reader is always free to examine the entire codes
at leisure.
2 A seismic trace is the recording of a single component geophone or hydrophone. For a multicomponent geophone,
3 Previous systems used a 16 bit word and variable gain. A four-bit gain word, an 11 bit mantissa, and a sign bit
of Calgary.
1.3. DYNAMIC RANGE AND SEISMIC DATA DISPLAY 7
1000 m offset
0.02
0.01
−0.01
−0.02
0 0.5 1 1.5 2 2.5 3
seconds
10 m offset
1
0.5
−0.5
−1
0 0.5 1 1.5 2 2.5 3
seconds
Figure 1.2: (Top) A real seismic trace recorded about 1000 m from a dynamite
shot. (Bottom) A similar trace recorded only 10 m from the same shot. See Code
Snippet 1.3.1
The close proximity of tracenear to a large explosion produces a very strong first arrival while
later information (at 3 seconds) has decayed by ∼72 decibels. (To gain familiarity with decibel
scales, it is useful to note that 6 db corresponds to a factor of 2. Thus 72 db represents about 72/6
∼ 12 doublings or a factor of 212 = 4096.) Alternatively, tracefar shows peak amplitudes that are
40db (a factor of 26.7 ∼ 100) weaker than tracenear.
Code Snippet 1.3.1. This code loads near and far offset test traces, computes the Hilbert envelopes
of the traces (with a decibel scale), and produces Figures 1.2 and 1.3.
1 clear; load testtrace.mat
2 subplot(2,1,1);plot(t,tracefar);
3 title(’1000 m offset’);xlabel(’seconds’)
4 subplot(2,1,2);plot(t,tracenear);
5 title(’10 m offset’);xlabel(’seconds’)
6 envfar = abs(hilbert(tracefar)); %compute Hilbert envelope
7 envnear = abs(hilbert(tracenear)); %compute Hilbert envelope
8 envdbfar=todb(envfar,max(envnear)); %decibel conversion
9 envdbnear=todb(envnear); %decibel conversion
10 figure
11 plot(t,[envdbfar envdbnear],’b’);xlabel(’seconds’);ylabel(’decibels’);
12 grid;axis([0 3 -140 0])
End Code
The first break time is the best estimate of the arrival time of the first seismic energy. For
tracefar this is about .380 seconds while for tracenear it is about .02 seconds. On each trace,
energy before this time cannot have originated from the source detonation and is usually taken as
8 CHAPTER 1. INTRODUCTION
−20
−40
−60
decibels
−80
−100
−120
−140
0 0.5 1 1.5 2 2.5 3
seconds
Figure 1.3: The envelopes of the two traces of Figure 1.2 plotted on a decibel scale.
The far offset trace is about 40 db weaker than the near offset and the total dynamic
range is about 120 db. See Code Snippet 1.3.1
an indication of ambient noise conditions. That is, it is due to seismic noise caused by wind, traffic,
and other effects outside the seismic experiment. Only for tracefar is the first arrival late enough
to allow a reasonable sampling of the ambient noise conditions. In this case, the average background
noise level is about 120 to 130 db below the peak signals on tracenear. This is very near the
expected instrument performance. It is interesting to note that the largest peaks on tracenear
appear to have square tops, indicating clipping, at an amplitude level of 1.0. This occurs because
the recording system gain settings were set to just clip the strongest arrivals and therefore distribute
the available dynamic range over an amplitude band just beneath these strong arrivals.
Dynamic range is an issue in seismic display as well as recording. It is apparent from Figure 1.2
that either signal fades below visual thresholds at about 1.6 seconds. Checking with the envelopes
on Figure 1.3, this suggests that the dynamic range of this display is about 40-50 db. This limitation
is controlled by two factors: the total width allotted for the trace plot and the minimum perceivable
trace deflection. In Figure 1.2 this width is about 1 inch and the minimum discernible wiggle is
about .01 inches. Thus the dynamic range is about 10−2 ∼ 40db, in agreement with the earlier
visual assessment.
It is very important to realize that seismic displays have limited dynamic range. This means
that what-you-see is not always what-you’ve-got. For example if a particular seismic display being
interpreted for exploration has a dynamic range of say 20 db, then any spectral components (i.e.
frequencies in the Fourier spectrum) that are more than 20 db down will not affect the display.
If these weak spectral components are signal rather than noise, then the display does not allow
the optimal use of the data. This is an especially important concern for wiggle trace displays of
multichannel data where each trace gets about 1/10 of an inch of display space.
1.3. DYNAMIC RANGE AND SEISMIC DATA DISPLAY 9
0.04
0.03
0.02
0.01
−0.01
−0.02
0.2 0.3 0.4 0.5 0.6 0.7 0.8
Figure 1.4: A portion of the seismic trace in Figure 1.2 is plotted in WTVA format
(top) and WT format (bottom). See Code Snippet 1.3.2
More popular than the wiggle-trace display is the wiggle-trace, variable-area (WTVA) display.
Function wtva (Code Snippet 1.3.2) was be used to create Figure 1.4 where the two display types
are contrasted using tracefar. The WTVA display fills-in the peaks of the seismic trace (or troughs
if polarity is reversed) with solid color. Doing this requires determining the zero crossings of the
trace which can be expensive if precision is necessary. Function wtva just picks the sample closest
to each zero crossing. For more precision, the final argument of wtva is a resampling factor which
causes the trace to be resampled and then plotted. Also, wtva works like MATLAB ’s low-level
line in that it does not clear the figure before plotting. This allows wtva to be called repeatedly in
a loop to plot a seismic section. The return values of wtva are MATLAB graphics handles for the
“wiggle” and the “variable area” that can be used to further manipulate their graphic properties.
(For more information consult your MATLAB reference.) The horizontal line at zero amplitude on
the WTVA plot is not usually present in such displays and seems to be a MATLAB graphic artifact.
Code Snippet 1.3.2. The same trace is plotted with wtva and plot. Figure 1.4 is the result.
1 clear;load testtrace.mat
2 plot(t,tracefar)
3 [h,hva]=wtva(tracefar+.02,t,’k’,.02,1,-1,1);
4 axis([.2 .8 -.02 .04])
End Code
Clipping was mentioned previously in conjunction with recording but also plays a role in display.
Clipping refers to the process of setting all values on a trace that are greater than a clip level equal
to that level. This can have the effect of moving the dynamic range of a plot display into the lower
amplitude levels. Figure 1.5 results from Code Snippet 1.3.3 and shows the effect of plotting the
10 CHAPTER 1. INTRODUCTION
0
6db 12db 18db 24db 30db 36db 42db 48db 54db 60db
0.5
1
seconds
1.5
2.5
3
−0.05 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4
Figure 1.5: The seismic trace in Figure 1.2 is plotted repeatedly with different clip
levels. The clip levels are annotated on each trace. See Code Snippet 1.3.3
same trace (tracefar at progressively higher clip levels. Function clip produces the clipped trace
that is subsequently rescaled so that the clip level has the same numerical value as the maximum
absolute value of the original trace. The effect of clipping is to make the weaker amplitudes more
evident at the price of completely distorting the stronger amplitudes. Clipping does not increase the
dynamic range, it just shifts the available range to a different amplitude band.
Code Snippet 1.3.3. This code makes Figure 1.5. The test trace is plotted repeatedly with pro-
gressively higher clip levels.
1 clear;load testtrace.mat
2 amax=max(abs(tracefar));
3 for k=1:10
4 trace_clipped=clip(tracefar,amax*(.5)ˆ(k-1))/((.5)ˆ(k-1));
5 wtva(trace_clipped+(k-1)*2*amax,t,’k’,(k-1)*2*amax,1,1,1);
6 text((k-1)*2*amax,.1,[int2str(k*6) ’db’]);
7 end
8 flipy;ylabel(’seconds’);
End Code
Exercise 1.3.1. Use MATLAB to load the test traces shown in Figure 1.2 and display them. By
appropriately zooming your plot, estimate the first break times as accurately as you can. What is the
approximate velocity of the material between the shot and the geophone? (You may find the function
simplezoom useful. After creating your graph, type simplezoom at the MATLAB prompt and then
use the left mouse mutton to draw a zoom box. A double-click will un-zoom.)
Exercise 1.3.2. Use MATLAB to load the test traces shown in Figure 1.2 and compute the envelopes
as shown in Code Snippet 1.3.1. For either trace, plot both the wiggle trace and its envelope and
show that the trace is contained within the bounds defined by ±envelope.
1.3. DYNAMIC RANGE AND SEISMIC DATA DISPLAY 11
0 0.1
0.1
0.15
0.2
0.2
0.3
0.25
0.4
seconds
seconds
0.5 0.3
0.6
0.35
0.7
0.4
0.8
0.45
0.9
1 0.5
0 100 200 300 400 500 600 700 800 900 1000 0 50 100 150 200 250 300 350 400 450 500
Figure 1.6: A synthetic seismic section plotted Figure 1.7: A zoom (enlargement of a portion)
in WTVA mode with clipping. See Code Snippet of Figure 1.6. Note the clipping indicated by
1.3.4 square-bottomed troughs
Exercise 1.3.3. What is the dynamic range of a seismic wiggle trace display plotted at 10 traces/inch?
What about 30 traces/inch?
Code Snippet 1.3.4. Here we create a simple synthetic seismic section and plot it as a WTVA
plot with clipping. Figure 1.6 is the result .
1 global NOSIG;NOSIG=1;
2 [r,t]=reflec(1,.002,.2);%make reflectivity
3 nt=length(t);
4 [w,tw]=wavemin(.002,20,.2);%make wavelet
5 s=convm(r,w);%make convolutional seismogram
6 ntr=100;%number of traces
7 seis=zeros(length(s),ntr);%preallocate seismic matrix
8 shift=round(20*(sin([1:ntr]*2*pi/ntr)+1))+1; %a time shift for each trace
9 %load the seismic matrix
10 for k=1:ntr
11 seis(1:nt-shift(k)+1,k)=s(shift(k):nt);
12 CHAPTER 1. INTRODUCTION
0
0.1
0.1
0.15
0.2
0.2
0.3
0.4 0.25
0.5 0.3
0.6
0.35
0.7
0.4
0.8
0.45
0.9
1 0.5
0 100 200 300 400 500 600 700 800 900 0 50 100 150 200 250 300 350 400 450 500
Figure 1.8: A synthetic seismic section plotted Figure 1.9: A zoom (enlargement of a portion)
as an image using plotimage . Compare with of Figure 1.8. Compare with Figure 1.7
Figure 1.6
12 end
13 x=(0:99)*10; %make an x coordinate vector
14 plotseis(seis,t,x,1,5,1,1,’k’);ylabel(’seconds’)
End Code
Code Snippet 1.3.4 illustrates the use of global variables to control plotting behavior. The global
NOSIG controls the appearance of a signature at the bottom of a figure created by plotseis . If
NOSIG is set to zero, then plotseis will annotate the date, time, and user’s name in the bottom
right corner of the plot. This is useful in classroom settings when many people are sending nearly
identical displays to a shared printer. The user’s name is defined as the value of another global
variable, NAME_ (the capitalization and the underscore are important). Display control through
global variables is also used in other utilities in the NumMethToolbox (see page 14). An easy way
to ensure that these variables are always set as you prefer is to include their definitions in your
startup.m file (see your MATLAB reference for further information).
The popularity of the WTVA display has resulted partly because it allows an assessment of the
seismic waveform (if unclipped) as it varies along an event of interest. However, because its dynamic
range varies inversely with trace spacing it is less suitable for small displays of densely spaced data
such as on a computer screen. It also falls short for display of multichannel Fourier spectra where
the wiggle shape is not usually desired. For these purposes, image displays are more suitable. In
this technique, the display area is divided equally into small rectangles (pixels) and these rectangles
are assigned a color (or gray-level) according to the amplitude of the samples within them. On
the computer screen, if the number of samples is greater than the number of available pixels, then
each pixel represents the average of many samples. Conversely, if the number of pixels exceeds the
number of samples, then a single sample can determine the color of many pixels, leading to a blocky
appearance.
Figures 1.8 and 1.9 display the same data as Figures 1.6 and 1.7 but used the function plotimage
in place of plotseis in Code Snippet 1.3.4. The syntax to invoke plotimage is plotimage(seis,t,x).
It may also be called like plotimage(seis) in which case it creates x and t coordinate vectors as
simply column and row numbers. Unlike plotseis , plotimage adorns its figure window with con-
trols to allow the user to interactively change the polarity, brightness, color map, and data scaling
1.3. DYNAMIC RANGE AND SEISMIC DATA DISPLAY 13
white white
1 1
0.8 0.8
0.6 0.6
graylevel
graylevel
0.4 0.4
0.2 0.2 40 80
20 60 100
black
0 0
black
0 10 20 30 40 50 60 70 0 10 20 30 40 50 60 70
level number level number
Figure 1.10: A 50% gray curve from seisclrs Figure 1.11: Various gray level curves from
(center) and its brightened(top) and darkened seisclrs for different gray pct values.
(bottom) versions.
scheme (though these controls are shown in these figures, in most cases in this book they will be
suppressed). The latter item refers to the system by which amplitudes are mapped to gray levels
(or colors).
Code Snippet 1.3.5. This example illustrates the behavior of the seisclrs color map. Figures
1.10 and 1.11 are created.
1 figure;
2 global NOBRIGHTEN
3
4 NOBRIGHTEN=1;
5 s=seisclrs(64,50); %make a 50% linear gray ramp
6 sb=brighten(s,.5); %brighten it
7 sd=brighten(s,-.5); %darken it
8 plot(1:64,[s(:,1) sb(:,1) sd(:,1)],’k’)
9 axis([0 70 -.1 1.1])
10 text(1,1.02,’white’);text(50,.02,’black’)
11 xlabel(’level number’);ylabel(’gray_level’);
12 figure; NOBRIGHTEN=0;
13 for k=1:5
14 pct=max([100-(k-1)*20,1]);
15 s=seisclrs(64,pct);
16 line(1:64,s(:,1),’color’,’k’);
17 if(rem(k,2))tt=.1;else;tt=.2;end
18 xt=near(s(:,1),tt);
19 text(xt(1),tt,int2str(pct))
20 end
21 axis([0 70 -.1 1.1])
22 text(1,1.02,’white’);text(50,-.02,’black’)
23 xlabel(’level number’);ylabel(’gray_level’);
End Code
By default, plotimage uses a gray level scheme, defined by seisclrs , which is quite nonlinear.
Normally, 64 separate gray levels are defined and a linear scheme would ramp linearly from 1 (white)
14 CHAPTER 1. INTRODUCTION
meters meters
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000
0 0
0.5 0.5
1 1
seconds
seconds
1.5 1.5
2 2
2.5 2.5
3 3
Figure 1.12: A raw shot record displayed using Figure 1.13: The same raw shot record as Figure
plotimage and maximum scaling 1.12 but displayed using mean scaling with a clip
level (parameter CLIP) of 4.
at level 1 to 0 (black) at level 64. However, this was found to produce overly dark images with little
event standout. Instead, seisclrs assigns a certain percentage, called gray pct, of the 64 levels
to the transition from white to black and splits the remaining levels between solid black and solid
white. For example, by default this percentage is 50 which means that levels 1 through 16 are solid
white, 17 through 48 transition from white to black, and 49 through 64 are solid black. The central
curve in Figure 1.10 illustrates this. As a final step, seisclrs brightens the curve using brighten
to produce the upper curve in Figure 1.10. Figure 1.11 shows the gray scales that result from varying
the value of gray pct.
Given a gray level or color scheme, function plotimage has two different schemes for determining
how seismic amplitudes are mapped to the color bins. The simplest method, called maximum scaling
determines the maximum absolute value in the seismic matrix, assigns this the color black, and the
negative of this number is white. All other values map linearly into bins in between. This works
well for synthetics, or well-balanced real data, but often disappoints for noisy or raw data because of
the data’s very large dynamic range. The alternative, called mean scaling, measures both the mean,
s̄ and standard deviation σs of the data and assigns the mean to the center of the gray scale. The
ends of the scale are assigned to the values s̄ ± CLIPσs where CLIP is a user-chosen constant. Data
values outside this range are assigned to the first or last gray-level bin thereby clipping them. Thus,
mean scaling centers the gray scale on the data mean and the extremes of the scale are assigned
to a fixed number of standard deviations from the mean. In both of these schemes, neutral gray
corresponds to zero (if the seisclrs color map is used).
Figure 1.12 displays a raw shot record using the maximum scaling method with plotimage .
(The shot record is contained in the file smallshot.mat, and the traces used in Figure 1.2 are the first
and last trace of this record.) The strong energy near the shot sets the scaling levels and is so much
stronger than anything else that most other samples fall into neutral gray in the middle of the scale.
On the other hand, figure 1.13 uses mean scaling (and very strong clipping) to show much more of
the data character.
In addition to the user interface controls on its window, the behavior of plotimage can be con-
trolled through global variables. Most variables defined in MATLAB are local and have a scope that
1.3. DYNAMIC RANGE AND SEISMIC DATA DISPLAY 15
0.1
0.1
0.2
0.2
0.3
0.3
0.4
0.4
seconds
seconds
0.5
0.5
0.6 0.6
0.7 0.7
0.8 0.8
0.9 0.9
1 1
Figure 1.14: A synthetic created to have a Figure 1.15: The same synthetic as Figure 1.14
6db decrease with each event displayed by displayed with plotseis .
plotimage . See Code Snippet 1.3.7.
is limited to the setting in which they are defined. For example, a variable called x defined in the
base workspace (i.e. at the MATLAB prompt ) can only be referenced by a command issued in the
base workspace. Thus, a function like plotimage can have its own variable x that can be changed
arbitrarily without affecting the x in the base workspace. In addition, the variables defined in a
function are transient, meaning that they are erased after the function ends. Global variables are
an exception to these rules. Once declared (either in the base workspace or a function) they remain
in existence until explicitly cleared. If, in the base workspace, x is declared to be global, then the
x in plotimage is still independent unless plotimage also declares x global. Then, both the base
workspace and plotimage address the same memory locations for x and changes in one affect the
other.
Code Snippet 1.3.6 illustrates the assignment of global variables for plotimage as done in the
author’s startup.m file. These variables only determine the state of a plotimage window at the
time it is launched. The user can always manipulate the image as desired using the user interface
controls. Any such manipulations will not affect the global variables so that the next plotimage
window will launch in exactly the same way.
Code Snippet 1.3.6. This is a portion of the author’s startup.m file that illustrates the setting of
global variables to control plotimage.
End Code
16 CHAPTER 1. INTRODUCTION
Code Snippet 1.3.7. This code creates a synthetic with a 6db decrease for each event and displays
it without clipping with both plotimage and plotseis. See Figures 1.14 and 1.15.
1 % put an event every 100 ms. Each event decreases in amplitude by 6 db.
2 r=zeros(501,1);dt=.002;ntr=100;
3 t=(0:length(r)-1)*dt;
4 amp=1;
5 for tk=.1:.1:1;
6 k=round(tk/dt)+1;
7 r(k)=(-1)ˆ(round(tk/.1))*amp;
8 amp=amp/2;
9 end
10 w=ricker(.002,60,.1);
11 s=convz(r,w);
12 seis=s*ones(1,ntr);
13 x=10*(0:ntr-1);
14 global SCALE_OPT GRAY_PCT
15 SCALE_OPT=2;GRAY_PCT=100;;
16 plotimage(seis,t,x);ylabel(’seconds’)
17 plotseis(seis,t,x,1,1,1,1,’k’);ylabel(’seconds’)
End Code
The dynamic range of image and WTVA displays is generally different. When both are displayed
without clipping using a moderate trace spacing (Figures 1.14 and 1.15), the dynamic range of
plotimage (≈30db) tends to be slightly greater than that of plotseis (≈20db). However, the
behavior of plotimage is nearly independent of the trace density while plotseis becomes useless
at high densities.
Exercise 1.3.4. Load the file smallshot.mat into your base workspace using the command load
and display the shot record with plotimage. Study the online help for plotimage and define the six
global variables that control its behavior. Practice changing their values and observe the results. In
particular, examine some of the color maps: hsv, hot, cool, gray, bone, copper, pink, white, flag, jet,
winter, spring, summer, and autumn. (Be sure to specify the name of the color map inside single
quotes.)
Exercise 1.3.5. Recreate Figures 1.14 and 1.15 using Code Snippet 1.3.7 as a guide. Experiment
with different number of traces (ntr) and different program settings to see how dynamic range is
affected.
0 100 200 300 400 500 600 700 800 900 1000 0 100 200 300 400 500 600 700 800 900 1000
0 0
0.1 0.1
0.2 0.2
0.3 0.3
0.4 0.4
0.5 0.5
0.6 0.6
0.7 0.7
0.8 0.8
0.9 0.9
1 1
Figure 1.16: This figure is displayed by Code Figure 1.17: After the user picks the first breaks
Snippet 1.3.8 before allowing the user to select in Figure 1.16, the picks are drawn on top of the
points. seismic data. See Code Snippet 1.3.8.
for the first window will be lost. If desired, this can be avoided by copying the pick list into another
variable prior to activating the second picking process.
A pick is made by clicking the left mouse button down at the first point of the pick, holding the
button down while dragging to the second point, and releasing the button. The pick will be drawn
on the screen in a temporary color and then drawn in a permanent color when the mouse is released.
The permanent color is controlled by the global variable PICKCOLOR. If this variable has not been
set, picks are drawn in red. A single mouse click, made without dragging the mouse and within the
bounds of the axes, will delete the most recent pick.
This picking facility is rudimentary in the sense that the picks are established without any
numerical calculation involving the data. Furthermore, there is no facility to name a set of picks
or to retain their values. Nevertheless, a number of very useful calculations, most notably raytrace
migration (see sections 4.2.2 and 4.3.3) are enabled by this. It is also quite feasible to implement a
more complete picking facility on this base.
coordinates of the nodes of the line to be drawn. If no z coordinates are given, they are assumed
zero which causes them to be drawn in the same z plane as the seismic (MATLAB graphics are
always 3D even when they appear to be only 2D). Plotting the line at z = 0 usually works but if
more mouse clicks are done on the figure, such as for zooming, then it can happen that the line may
get re-drawn first with the seismic on top of it. This causes the line to “disappear”. Instead, it is
more robust to give the line a vector of z coordinates of unity so that it is guaranteed to always be
in front of the seismic.
In addition to (x, y, z) coordinates, line accepts an arbitrary length list of (attribute, property)
pairs that define various features of the line. In this case, its color is set to red, its line thickness is
set to twice normal, and markers are added to each point. There are many other possible attributes
that can be set in this way. To see a list of them, issue the command get(h) where h is the handle
of a line and MATLAB will display a complete property list. Once the line has been drawn, you can
alter any of its properties with the set command. For example, set(h,’color’,’c’) changes the
line color to cyan and set(h,’ydata’,get(h,’ydata’)+.1) shifts the line down .1 seconds.
Code Snippet 1.3.8. This example displays the small sample shot record using plotimage and
then uses ginput to allow the user to enter points. These points are plotted on top of the seismic
data using line. The results are in Figures 1.16 and 1.17.
1 load smallshot
2 global SCALE_OPT CLIP
3 SCALE_OPT=1;CLIP=1;
4 plotimage(seis,t,x)
5 axis([0 1000 0 1])
6 [xpick,tpick]=ginput;
7 h=line(xpick,tpick,ones(size(xpick)),’color’,’r’,...
8 ’linewidth’,2,’marker’,’*’);
End Code
1.4.1 Scripts
MATLAB is designed to provide a highly interactive environment that allows very rapid testing
of ideas. However, unless some care is taken, it is very easy to create results that are nearly
irreproducible. This is especially true if a long sequence of complex commands has been entered
directly at the MATLAB prompt. A much better approach is to type the commands into a text file
and execute them as a script. This has the virtue that a permanent record is maintained so that
the results can be reproduced at any time by simply re-executing the script.
A script is the simplest possible MATLAB programming structure. It is nothing more than a
sequence of syntactically correct commands that have been typed into a file. The file name must
end in .m and it must appear in the MATLAB search path. When you type the file name (without
the .m) at the MATLAB prompt, MATLAB searches its path for the first .m file with that name
and executes the commands contained therein. If you are in doubt about whether your script is the
1.4. PROGRAMMING TOOLS 19
first so-named file in the path, use the which command. Typing which foo causes MATLAB to
display the complete path to the file that it will execute when foo is typed.
A good practice is to maintain a folder in MATLAB’s search path that contains the scripts
associated with each logically distinct project. These script directories can be managed by creating
a further set of control scripts that appropriately manipulate MATLAB’s search path. For example,
suppose that most of your MATLAB tools are contained in the directory
C:\Matlab\toolbox\local\mytools
and that you have project scripts stored in
C:\My Documents\project1\scripts
and
C:\My Documents\project2\scripts.
Then, a simple management scheme is to include something like
1 global MYPATH
2 if(isempty(MYPATH))
3 p=path;
4 path([p ’;C:\MatlabR11\toolbox\local\mytools’]);
5 MYPATH=path;
6 else
7 path(MYPATH);
8 end
in your startup.m file and then create a script called project1.m that contains
1 p=path;
2 path([p ’;C:\My Documents\project1\scripts’])
and a similar script for project2. When you launch MATLAB your startup.m establishes your base
search path and assigns it to a global variable called MYPATH. Whenever you want to work on
project1, you simply type project1 and your base path is altered to include the project1 scripts.
Typing startup again restores the base path and typing project2 sets the path for that project.
Of course you could just add all of your projects to your base path but then you cannot have scripts
with the same name in each project.
Exercise 1.4.1. Write a script that plots wiggle traces on top of an image plot. Test your script
with the synthetic section of Code Snippet 1.3.4. (Hint: plotseis allows its plot to be directed to an
existing set of axes instead of creating a new one. The “current” axes can be referred to with gca.)
1.4.2 Functions
The MATLAB function provides a more advanced programming construct than the script. The
latter is simply a list of MATLAB commands that execute in the base workspace. Functions execute
in their own independent workspace that communicates with the base workspace through input and
output variables.
[output_variable_list]=function_name(input_variable_list).
Rectangular brackets [ . . . ] enclose the list of output variables, while regular parentheses enclose the
input variables. Either list can be of any length, and entries are separated by commas. For clip ,
there is only one output variable and the [ . . . ] are not required. When the function is invoked,
either at the command line or in a program, a choice may be made to supply only the first n input
variables or to accept only the first m output variables. For example, given a function definition like
[x1,x2]=foo(y1,y2), it make be legally invoked with any of
[a,b]=foo(c,d);
[a,b]=foo(c);
a=foo(c,d);
a=foo(c);
The variable names in the function declaration have no relation to those actually used when the
function is invoked. For example, when a=foo(c) is issued the variable c becomes the variable (y1)
within foo’s workspace and similarly, foo’s x1 is returned to become the variable a in the calling
workspace. Though it is possible to call foo without specifying y2 there is no simple way to call it
and specify y2 while omitting y1. Therefore, it is advisable to structure the list of input variables
such that the most important ones appear first. If an input variable is not supplied, then it must be
assigned a default value by the function.
Code Snippet 1.4.1. A simple MATLAB function to clip a signal.
1 function trout=clip(trin,amp);
2 % CLIP performs amplitude clipping on a seismic trace
3 %
4 % trout=clip(trin,amp)
5 %
6 % CLIP adjusts only those samples on trin which are greater
7 % in absolute value than ’amp’. These are set equal to amp but
8 % with the sign of the original sample.
9 %
10 % trin= input trace
11 % amp= clipping amplitude
12 % trout= output trace
13 %
14 % by G.F. Margrave, May 1991
15 %
16 if(nargin˜=2)
17 error(’incorrect number of input variables’);
18 end
19 % find the samples to be clipped
20 indices=find(abs(trin)>amp);
21 % clip them
22 trout=trin;
23 trout(indices)=sign(trin(indices))*amp;
End Code
In the example of the function clip there are two input variables, both mandatory, and one out-
put variable. The next few lines after the function definition are all comments (i.e. non-executable)
as is indicated by the % sign that precedes each line. The first contiguous block of comments follow-
ing the function definition constitutes the online help. Typing help clip at the MATLAB prompt
will cause these comment lines to be echoed to your MATLAB window.
The documentation for clip follows a simple but effective style. The first line gives the function
name and a simple synopsis. Typing help directory name, where directory name is the name of a
1.4. PROGRAMMING TOOLS 21
directory containing many .m files, causes the first line of each file to be echoed giving a summary of
the directory. Next appears one or more function prototypes that give examples of correct function
calls. After viewing a function’s help file, a prototype can be copied to the command line and edited
for the task at hand. The next block of comments gives a longer description of the tasks performed
by the function. Then comes a description of each input and output parameter and their defaults,
if any. Finally, it is good form to put your name and date at the bottom. If your code is to be
used by anyone other than you, this is valuable because it establishes who owns the code and who
is responsible for fixes and upgrades.
Following the online help is the body of the code. Lines 16-18 illustrate the use of the auto-
matically defined variable nargin that is equal to the number of input variables supplied at calling.
clip requires both input variables so it aborts if nargin is not two by calling the function error.
The assignment of default values for input variables can also be done using nargin. For example,
suppose the amp parameter were to be allowed to default to 1.0. This can be accomplished quite
simply with
if(nargin<2) amp=1; end
Lines 20, 22, and 23 actually accomplish the clipping. They illustrate the use of vector addressing
and will be discussed more thoroughly in Section 1.5.1.
Code Snippet 1.4.2. Assume that a matrix seis already exists and that nzs and nzt are already
defined as the sizes of zero pads in samples and traces respectively. Then a padded seismic matrix
is computed as:
1 %pad with zero samples
2 [nsamp,ntr]=size(seis);
3 seis=[seis;zeros(nzs,ntr)];
4 %pad with zero traces
5 seis=[seis zeros(nsamp+nzs,nzt)];
End Code
The use of [ . . . ] on lines 3 and 5 is key here. Unless they are being used to group the return
variables from a function, [ . . . ] generally indicate that a new matrix is being formed from the
concatenation of existing matrices. On line 3, seis is concatenated with a matrix of zeros that
has the same number of traces as seis but nzs samples. The semi-colon between the two matrices
is crucial here as it is the row separator while a space or comma is the column separator. If the
elements in [ ] on line 3 were separated by a space instead of a semi-colon, MATLAB would emit
the error message:
??? All matrices on a row in the bracketed expression must have the same
number of rows.
22 CHAPTER 1. INTRODUCTION
This occurs because the use of the column separator tells MATLAB to put seis and zeros(nxs,ntr)
side-by-side in a new matrix; but, this is only possible if the two items have the same number of
rows. Similarly, if a row separator is used on line 5, MATLAB will complain about “All matrices
on a row in the bracketed expression must have the same number of columns”.
Another common error is an assignment statement in which the matrices on either side of the
equals sign do not evaluate to the same size matrix. This is a fundamental requirement and calls
attention to the basic MATLAB syntax: MatrixA = MatrixB. That is, the entities on both sides of
the equals sign must evaluate to matrices of exactly the same size. The only exception to this is
that the assignment of a constant to a matrix is allowed.
One rather insidious error deserves special mention. MATLAB allows variables to be “declared”
by simply using them in context in a valid expression. Though convenient, this can cause problems.
For example, if a variable called plot is defined, then it will mask the command plot. All further
commands to plot data (e.g. plot(x,y)) will be interpreted as indexing operations into the matrix
plot. More subtle,√if the letter i is used as the index of a loop, then it can no longer serve its
predefined task as −1 in complex arithmetic. Therefore some caution is called for in choosing
variable names, and especially, the common Fortran practice of choosing i as a loop index is to be
discouraged.
The most subtle errors are those which never generate an overt error message but still cause
incorrect results. These logical errors can be very difficult to eliminate. Once a program executes
successfully, it must still be verified that it has produced the expected results. Usually this means
running it on several test cases whose expected behavior is well known, or comparing it with other
codes whose functionality overlaps with the new program.
The MATLAB debugger is a very helpful facility for resolving run-time errors, and it is simple
to use. There are just a few essential debugger commands such as: dbstop, dbstep, dbcont,
dbquit, dbup, dbdown. A debugging session is typically initiated by issuing the dbstop command
to tell MATLAB to pause execution at a certain line number, called a breakpoint. For example
dbstop at 200 in plotimage will cause execution to pause at the executable line nearest line 200
in plotimage . At this point, you may choose to issue another debug command, such as dbstep
which steps execution a line at a time, or you may issue any other valid MATLAB command. This
means that the entire power of MATLAB is available to help discover the problem. Especially
when dealing with large datasets, it is often very helpful to issue plotting commands to graph the
intermediate variables.
As an example of the debugging facility, consider the problem of extracting a slice from a matrix.
That is, given a 2D matrix, and a trajectory through it, extract a sub-matrix consisting of those
samples within a certain half-width of the trajectory. The trajectory is defined as a vector of row
numbers, one-per-column, that cannot double back on itself. A first attempt at creating a function
for this task might be like that in Code Snippet 1.4.3.
Code Snippet 1.4.3. Here is some code for slicing through a matrix along a trajectory. Beware,
this code generates an error.
1 function s=slicem(a,traj,hwid)
2
3 [m,n]=size(a);
4 for k=1:n %loop over columns
5 i1=max(1,traj(k)-hwid); %start of slice
6 i2=min(m,traj(k)+hwid); %end of slice
7 ind=(i1:i2)-traj(k)+hwid; %output indices
8 s(ind,k) = a(i1:i2,k); %extract the slice
9 end
End Code
1.4. PROGRAMMING TOOLS 23
a=((5:-1:1)’)*(ones(1,10))
a =
5 5 5 5 5 5 5 5 5 5
4 4 4 4 4 4 4 4 4 4
3 3 3 3 3 3 3 3 3 3
2 2 2 2 2 2 2 2 2 2
1 1 1 1 1 1 1 1 1 1
traj=[1:5 5:-1:1]
traj =
1 2 3 4 5 5 4 3 2 1
The samples of a that lie on the trajectory are 5 4 3 2 1 1 2 3 4 5. Executing slicem with
the command s=slicem(a,traj,1) generates the following error:
This error message suggests that there is a problem with indexing on line 8; however, this could
be occurring either in indexing into a or s. To investigate, issue the command “dbstop at 8 in
slicem” and re-run the program. Then MATLAB stops at line 8 and prints
8 s(ind,k) = a(i1:i2,k);
K
Now, suspecting that the index vector ind is at fault we list it to see that it contains the values
1 2 and so does the vector i1:i2. These are legal indices. Noting that line 8 is in a loop, it seems
possible that the error occurs on some further iteration of the loop (we are at k=1). So, execution
is resumed with the command dbcont. At k=2, we discover that ind contains the values 0 1 2
while i1:i2 is 1 2 3. Thus it is apparent that the vector ind is generating an illegal index of 0. A
moments reflection reveals that line 7 should be coded as ind=(i1:i2)-traj(k)+hwid+1; which
includes an extra +1. Therefore, we exit from the debugger with dbquit, make the change, and
re-run slicem to get the correct result:
s =
0 5 4 3 2 2 3 4 5 0
5 4 3 2 1 1 2 3 4 5
4 3 2 1 0 0 1 2 3 4
The samples on the trajectory appear on the central row while the other rows contain neighboring
values unless such values exceed the bounds of the matrix a, in which case a zero is returned. A more
polished version of this code is found as slicemat in the NumMethToolbox . Note that a simpler
debugger procedure (rather than stopping at line 8) would be to issue the command “dbstop if
error” but the process illustrated here shows more debugging features.
24 CHAPTER 1. INTRODUCTION
Code Snippet 1.5.1. A Fortran programmer new to MATLAB would probably code clip in this
way.
1 function trout=fortclip(trin,amp)
2
3 for k=1:length(trin)
4 if(abs(trin(k)>amp))
5 trin(k)=sign(trin(k))*amp;
6 end
7 end
8
9 trout=trin;
End Code
Code Snippet 1.5.2. This script compares the execution time of clip and fortclip
1 [r,t]=reflec(1,.002,.2);
2
3 tic
4 for k=1:100
5 r2=clip(t,.05);
6 end
7 toc
8
9 tic
10 for k=1:100
11 r2=fortclip(t,.05);
12 end
13 toc
End Code
In contrast to the vector coding style just described Code Snippet 1.5.1 shows how clip might be
coded by someone stuck-in-the-rut of scalar addressing (a C or Fortran programmer). This code is
logically equivalent to clip but executes much slower. This is easily demonstrated using MATLAB’s
built in timing facility. Code Snippet 1.5.2 is a simple script that uses the tic and toc commands
for this purpose. tic sets an internal timer and toc writes out the elapsed time since the previous
tic. The loops are executed 100 times to allow minor fluctuations to average out. On the second
execution of this script, MATLAB responds with
1.5. PROGRAMMING FOR EFFICIENCY 25
elapsed_time = 0.0500
elapsed_time = 1.6000
which shows that clip is 30 times faster than fortclip . (On the first execution, clip is only
about 15 times faster because MATLAB spends some time doing internal compiling. Run time tests
should always be done a number of times to allow for effects like this.)
A simple blunder that can slow down fortclip even more is to write it like
Code Snippet 1.5.3. An even slower version of fortclip
1 function trout=fortclip(trin,amp)
2
3 for k=1:length(trin)
4 if(abs(trin(k)>amp))
5 trout(k)=sign(trin(k))*amp;
6 end
7 end
End Code
This version of fortclip , still produces correct results, but is almost 50 times slower than clip .
The reason is that the output trace, trout, is being addressed sample-by-sample but it has not been
pre-allocated. This forces MATLAB to resize the vector each time through the loop. Such resizing
is slow and may require MATLAB to make repeated requests to the operating system to grow its
memory.
Exercise 1.5.1. Create a version of fortclip and verify the run time comparisons quoted here.
Your computer may give different results. Show that the version of fortclip in Code Snippet 1.5.3
can be made to run as fast as that in Code Snippet 1.5.1 by using the zeros function to pre-allocate
trout.
End Code
This works correctly but is needlessly slow. An experienced MATLAB programmer would know that
the operation of matrix-times-diagonal matrix results in a new matrix whose columns are scaled by
the corresponding element of the diagonal matrix. You can check this out for yourself by a simple
MATLAB exercise:
ans =
1 2 3
1 2 3
1 2 3
1 2 3
Thus the MATLAB programmer would write Code Snippet 1.5.4 with the very simple single line
in Code Snippet 1.5.5.
Code Snippet 1.5.5. A MATLAB programmer would write the example of Code Snippet 1.5.4 in
this way:
1 seis=seis*diag(scales);
End Code
Tests indicate that Code Snippet 1.5.5 is about ten times as fast as Code Snippet 1.5.4.
Vector programming is also facilitated by the fact that MATLAB’s mathematical functions are
automatically set up to operate on arrays. This includes functions like sin, cos, tan, atan,
atan2, exp, log, log10, sqrt and others. For example, if t is a time coordinate vector then
sin(2*pi*10*t) is a vector of the same size as t whose entries are a 10 Hz. sine wave evaluated
at the times in t. These functions all perform element-by-element operations on matrices of any
size. Code written in this way is actually more visually similar to the corresponding mathematical
formulae than scalar code with its many explicit loops.
Some of MATLAB’s vectorized functions operate on matrices but return matrices of one di-
mension smaller. The functions sum, cumsum, mean, main, max, std all behave in this manner.
For example, if trace is a vector and seis is a matrix, then mean(trace) results in a scalar that
is the mean value of the vector while mean(seis) results in a row vector containing the mean
value of each column of seis. The mean value of the entire matrix, seis, can be calculated with
mean(mean(seis)).
seis(:,10) refers to the 10th trace. The colon indicates that all rows (samples) are desired.
seis(520:2:720,:) grabs very other sample between sample number 520 and number 720 on all
traces.
Of course, the indices need not be entered as explicit numbers but could be a variable instead as in
the case of the return from find in Code Snippet 1.4.1.
The last item needs more explanation. Regardless of the dimensions of a matrix, it can always
be referred to as a single column vector using a single index. For example, the two dimensional
matrix, seis, can be indexed like seis(irow,icol) and as seis(ipos). This is possible because
computer memory is actually a one-dimensional space and MATLAB stores a matrix in this space
in column order. This means that the actual storage pattern of seis is [seis(1,1) seis(2,1)
...seis(nrows,1) seis(1,2) seis(2,2) ... ...seis(nrows,ncols)]. Thus the last sample of
column one and the first sample of column two are contiguous in memory. This means that the first
sample of column two can be addressed either by seis(1,2) or by seis(nrows+1). In accordance
with this formulation, MATLAB allows any matrix to be entered into a formula as a equivalent
column vector using the notation seis(:).
with isnan. Also, the debug command dbstop if naninf is helpful to stop execution at the point
that NaN or Inf occurs in any function.
MATLAB uses double precision floating point arithmetic for its calculations. This is actually
a higher standard of precision than has traditionally been used in seismic data processing. Most
processing systems have been written in Fortran or C and use single precision arithmetic. MATLAB
supplies the internal constant eps whose value is the smallest positive floating point number available
on your machine. In a loose sense, eps represents the machine “noise level” and is about 2.2204x10−16
on the author’s computer. (How many decibels down from unity is eps?) As a general rule, tests for
the equality of two floating point variables should be avoided. Instead of a test like if(x==y) consider
something like if(abs(x-y)<10*eps). eps can also be useful in computing certain limits that
involve a division by zero. For example, computing a sinc function with x=0:.1:pi;y=sin(x)./x;
generates a zero divide and y(1) is NaN. However, x=0:.1:pi;y=sin(x+eps)./(x+eps); avoids
the zero divide and results in y(1) of 1.
Velocity
Exploration seismology is literally overflowing with velocities. Just to name a few, there are interval
velocity, instantaneous velocity, apparent velocity, rms velocity, average velocity, mean velocity,
stacking velocity, horizontal velocity, vertical velocity, phase velocity, group velocity, P-wave velocity,
S-wave velocity, migration velocity, weathering velocity, and almost certainly others. This chapter is
meant to bring some order to this chaos of velocities. For starters, there is a fundamental distinction
between physical velocities and velocity measures. The former refers to velocities that are actually
the speed at which some physical wave propagates. Examples are instantaneous velocity, P- and
S-wave velocities, phase and group velocity. On the other hand, velocity measures are typically
quantities derived from data analysis that have the physical dimensions of velocity but are related to
physical velocities in some indirect (and possibly unknown) fashion. Examples of velocity measures
include average, mean, and rms velocities, interval velocity, stacking velocity, apparent velocity, and
migration velocity. In contrast to physical velocities, it cannot generally be expected that a physical
wave actually propagates at the speed of one of these velocity measures.
Using the measured data for the analysis of velocity and the application to the data of corrections
that depend upon velocity are fundamental to seismic processing. Velocity, together with the geom-
etry of the seismic acquisition and the shapes of the reflectors, causes the characteristic traveltime
shapes (such as hyperbolae) in the data. Spatial variations in velocity cause distortions from the
simple, canonical shapes, and these distortions must be accounted for in processing. Sometimes ve-
locity information can be obtained from supporting data such as well logs, core measurements, and
geologic maps but this information is inherently sparse with respect to the coverage of the seismic
data itself. Ultimately, the only way to obtain the spatially dense velocity information required in
processing is from the data itself. This is known as velocity analysis and is the source of some of
the velocity measures such as stacking and migration velocities. The interpretation of the velocity
measures such that they can be converted into measurements of physical velocities is called velocity
inversion and is an ongoing problem. A careful discussion of the definitions of these measures is
essential to gain understanding of the complexities of velocity inversion.
For simplicity, consider the case of P-wave propagation in a heterogeneous earth. (S-waves can
be treated similarly.) A velocity measure will be defined by either relating it mathematically to the
actual P-wave speed or by describing how it can be derived from seismic data.
29
30 CHAPTER 2. VELOCITY
500
1000
depth(m)
1500
2000
2500
3000
1500 2000 2500 3000 3500 4000
velocity(m/s)
The universal velocity function is the name given to this function when v0 = 1800 m/sec and
c = .6 sec-1 . Figure 2.1 shows this function for a depth range of 0 to 3000 m. The name originates
in the days of analog computing (the 1950s and 1960s) when changing a velocity function was a
considerable labor. For example, normal moveout was removed by a “computer” (who was a human
being) who “traced the trace” (hence the name trace for a seismic recording) with a stylus that
was connected to a pen by a lever arm driven by an eccentric cam. The pen would redraw the
trace with normal moveout removed. However, the shape of the cam was determined by the velocity
function and changing functions often meant going to the machine shop to cut a new cam. Thus the
possibility of a universal function was quite attractive. And, it actually works quite well in many
sedimentary basins.
2.2. VERTICAL TRAVELTIME: τ 31
The total one-way traveltime from the surface to depth, z, is simply the sum of many such small
contributions. In the limit as dz → 0 this becomes the integral
z
dz̃
τ (z) = . (2.3)
0 vins (z̃)
The z dependence appears as the upper limit of the integral while z̃ is just a dummy variable of
integration. Defined in this way, τ (z) is a function that increases monotonically with z. As such it
is guaranteed to have an inverse. That is, given τ (z) it is always possible to find z(τ ) and vice-versa.
Function τ (z) is called a time-depth curve and can be used to find depth given vertical traveltime
or the reverse.
Continuing with the universal velocity function of Equation 2.1, the vertical traveltime of equation
(2.3) becomes
z
dz̃ 1 v0 +cz dξ 1 cz
τ (z) = = = ln 1 + . (2.4)
0 v(z̃) c v0 ξ c v0
Thus a linear variation of velocity with depth yields a logarithmic time-depth curve. It is a simple
matter to invert Equation 2.4 for z(τ ) to get
v0 cτ
z(τ ) = [e − 1] . (2.5)
c
For the case of the universal velocity function, the curves τ (z) and z(τ ) are shown in Figures 2.2
and 2.3 respectively. The graph of z(τ ) is just the transpose of the graph of τ (z).
0 0
0.2 500
0.4 1000
depth (m)
time (s)
0.6 1500
0.8 2000
1 2500
1.2 3000
0 500 1000 1500 2000 2500 3000 0 0.2 0.4 0.6 0.8 1 1.2
depth (m) time (s)
Figure 2.2: The time-depth curve (τ (z)) of equa- Figure 2.3: The depth-time curve (z(τ )) of equa-
tion (2.4) for the case of v0 = 1800m/sec and tion (2.4) for the case of v0 = 1800 m/sec and
c = .6/sec c = .6 sec-1
Comparing equations (2.3) and (2.7) shows that knowledge of vins (z) allows computation of τ (z)
and knowing vins (τ ) enables computation of z(τ ). In practice, vins (z) might be directly obtained
from a sonic well log while vins (τ ) can be estimated from the analysis of stacking velocities.
For the universal velocity function of equation (2.1) the vertical traveltime results in the loga-
rithmic expression given in equation (2.4) and z(τ ) is given by equation (2.5). Thus vins (τ ) becomes
v
0
vins (τ ) = vins (z(τ )) = v0 + c (ecτ − 1) (2.8)
c
which reduces to
vins (τ ) = v0 ecτ . (2.9)
Thus, when vins is linear with depth it is actually exponential with vertical traveltime. For case of
the universal velocity function, vins (τ ) is plotted in Figure 2.4.
0.2
0.4
time (sec)
0.6
0.8
1.2
2000 2500 3000 3500
velocity (m/s)
Figure 2.4: For the universal velocity function shown in Figure 2.1, vins is an exponential
function of vertical traveltime.
while vave (τ ) is τ
z(τ ) 1
vave (τ ) = = vins (τ̃ )dτ̃ . (2.11)
τ τ 0
Equation (2.11) shows that the average velocity is a mathematical average only with respect to
vertical traveltime. (Recall from calculus, that the average (or mean) value of a function f (x) over
x
the interval from 0 → x is x−1 0 f (x )dx .) When considered as a function of depth, the average
velocity is not a mathematical average.
For the running example of linear variation of velocity with depth, vave (z) becomes
cz
vave (z) = (2.12)
cz
ln 1 + v0
0
0
0.2
500
1000 0.4
time (sec)
depth (m)
1500 0.6
2000 0.8
v v
ins ins
2500 1 v
v ave
ave
3000 1.2
2000 2500 3000 3500 2000 2500 3000 3500
velocity (m/s) velocity (m/s)
Figure 2.5: vave and vins are shown versus depth Figure 2.6: vave and vins are shown versus time
for the case of the universal velocity functions for the case of the universal velocity functions
(equation (2.1)) with v0 = 1800 m/s and c = (equation (2.1)) with v0 = 1800 m/s and c =
.6 sec-1 .6 sec-1
which shows that the mean velocity increases at half the rate of vins (z). Of course, vmean (z) can be
re-expressed as a function of τ by substituting z(τ ), this results in
v0
vmean (τ ) = [1 + ecτ ] . (2.16)
2
Equation (2.15) is compared with vave (z) and vins (z) in Figure 2.7. vmean (z) and vins (z) are
both linear but vave (z) is not. Also, vmean (z) is always greater than vave (z) which will be proven
true in the next section.
Two important relationships exist between the three velocity measures: vave , vmean , and vrms .
First, they are linked by the relation
2
vrms = vave vmean (2.18)
2.7. INTERVAL VELOCITY: VIN T 35
500
1000
depth (m)
1500
2000 v
ins
vmean
2500 v
ave
3000
2000 2500 3000 3500
velocity (m/s)
Figure 2.7: vmean and vave are compared with vins for the case of the linear increase of
vins with z.
which means that only two of the three velocity measures can be considered
independent. The proof
2
is straight forward vrms = τ −1 vins
2
dτ = τz z1 vins [vins dτ ] = τz [ z1 vins dz] = vave vmean .
The second relationship follows from a mathematical inequality known as Schwartz’s Inequality.
This result is quite general and confirms the intuition that the square root of the average of the
squares of a set of numbers is always greater than the direct average
of the numbers
themselves.
That is, let {xk } be an arbitrary set of n real numbers, then n−1 k x2k ≥ n−1 k xk where the
equality only occurs if all xk are identical.
For the continuing example of vins linear with depth, equation (2.9) showed that vint was expo-
nential. Thus vrms becomes
1 τ 2 2cτ̃ v2
2
vrms (τ ) = v0 e dτ̃ = 0 e2cτ − 1 . (2.19)
τ 0 2cτ
Figure 2.8 compares vrms vave and vins for the universal velocity function. Though vrms is
always greater than vave the difference is not greater than 2%.
Exercise 2.6.1. Use MATLAB to numerically demonstrate Schwartz’s inequality. Use rand to
generate a vector of random numbers with zero mean. Compute both the mean and the rms average
of these numbers. Repeat the process for many different vectors with steadily increasing length. Plot
a graph of these two averages versus length. Repeat your work using random numbers that are all
positive with a mean value somewhere near that expected for seismic velocities. Show that if the
velocity vector has identical entries then the rms average and the mean are equal but if one value
differs from the rest, then the rms average exceeds the mean.
0.2
0.4
time (sec)
0.6
0.8
v
rms
vins
1 v
ave
1.2
2000 2500 3000 3500
velocity (m/s)
Figure 2.8: vrms and vave are compared with vins for the case of the linear increase of vins
with z.
For example, the average and rms velocities across an interval defined by τ1 and τ2 are simply
τ2
1
vave (τ2 , τ1 ) = vins (τ̃ )dτ̃ (2.20)
τ2 − τ1 τ1
and τ2
2 1 2
vrms (τ2 , τ1 ) = vins (τ̃ )dτ̃ . (2.21)
τ2 − τ1 τ1
These quantities are functions of both the upper and lower bounds of the integrals. In the limit, as
the interval shrinks so small that vins (τ )can be considered constant, both interval velocities approach
the instantaneous velocity.
It follows directly from the definition of average velocity (equation (2.10)) that the average
velocity across a depth interval is just the ratio of the depth interval to the vertical traveltime across
the interval. That is
z2 − z1 ∆z
vave (τ2 , τ1 ) = = . (2.22)
τ2 − τ1 ∆τ
Thus, if the range from 0 to z is divided into n finite intervals defined by z0 , z1 , z2 , . . . , zn−1 , zn
(where z0 = 0 and zn = z) then
n n
z(τ ) k=1 [zk − zk−1 ] 1
vave (τ ) = = = [τk − τk−1 ]vave (τk , τk−1 ) (2.23)
τ τ τ
k=1
(2.24) is the average velocity of the kth finite interval and is thus the time average of vins across the
interval. Of course, if vins (τ ) is constant across each interval then vk is an instantaneous velocity
and this distinction vanishes.
Equation (2.24) can be used in a number of ways. Most obviously, it shows how to combine a
set of local average velocities to obtain the macro average velocity across a larger interval. Also, it
can be used to estimate a local average velocity given two macro-average velocities. Suppose the
average velocities vave1 and vave2 from z = 0 to depths z1 and z2 are known. Then an expression for
the average velocity, from z1 → z2 or equivalently from τ1 → τ2 , follows from equation (2.24) and is
1
vave (τ2 , τ1 ) = [τ2 vave2 − τ1 vave1 ] . (2.25)
τ2 − τ1
The two macro averages are each weighted by their time intervals and then the shallower average is
subtracted from the deeper. This difference is then divided by the time interval of the local average.
If vave1 , τ1 , vave2 , andτ2 are all measured without error, then this process (i.e. equation 2.25)
works perfectly. However, in a real case, errors in the measurements can lead to wildly incorrect
estimates of vave (τ2 , τ1 ). With noisy data, there is no guarantee that the term [τ2 vave2 − τ1 vave1 ]
will always be positive leading to the possibility of negative velocity estimates. Also, the division by
τ2 − τ1 can be unstable if the two times are very close together.
The discussion so far has been only for average velocities across an interval. The entire derivation
above can be repeated for rms velocities with similar results but a few more subtleties arise in
interpretation. Rather than repeating the derivation just given with ‘rms’ in place of ‘ave’, it is
c b c
instructive to consider an alternative approach. It is a property of integrals that a = a + b where
a < b < c. Given τ0 < τ1 < τ2 and applying this rule to equation (2.17) results in
τ1 τ2
2 1 2 2
vrms (τ2 , τ0 ) = vins (τ̃ )dτ̃ + vins (τ̃ )dτ̃ . (2.26)
τ2 − τ0 τ0 τ1
Recognizing τ the integrals in [ . . . ] as rms interval velocities squared multiplied by their interval
2 2
times (i.e. τ01 vins (τ̃ )dτ̃ = [τ1 − τ0 ]vrms (τ1 , τ0 ) and similarly for the other integral) leads to
2 1 2 2
vrms (τ2 , τ0 ) = (τ1 − τ0 )vrms (τ1 , τ0 ) + (τ2 − τ1 )vrms (τ2 , τ1 ) . (2.27)
τ2 − τ0
For the case of n subdivisions between τ2 and τ0 this generalizes to
n
2 1
vrms (τ2 , τ0 ) = vk2 ∆τk (2.28)
τ2 − τ0
k=1
n
where vk = vrms (τk , τk−1 ) and, as before, ∆τk = τk − τk−1 and τ = k=1 ∆τk . This is the rms
equivalent to equation (2.24) and all of the comments made previously about vave apply here for
vrms . In particular, equation (2.28) should be thought of a combining interval rms velocities into
a macro rms velocity. Only in the case when vins does not vary significantly across an interval can
the vk in equation (2.28) be considered to be instantaneous velocities.
Equation (2.27) is the addition rule for rms velocities. To combine rms velocities, the squared ve-
locities are added and each must be weighted by its time interval. Equation (2.27) can be rearranged
38 CHAPTER 2. VELOCITY
to give an expression for estimating a local rms velocity from two macro velocities
2 1 2 2
vrms (τ2 , τ1 ) = (τ2 − τ0 )vrms (τ2 , τ0 ) − (τ1 − τ0 )vrms (τ1 , τ0 ) (2.29)
τ2 − τ1
or, with simplified notation
2 1 2 2
vrms (τ2 , τ1 ) = τ2 vrms2 − τ1 vrms1 . (2.30)
τ2 − τ1
In this expression τ0 has been set to 0 and vrms2 = vrms (τ2 , τ0 ) and similarly for vrms1 . Equation
(2.30) is often called the Dix equation for interval rms velocities because it was C.H. Dix (Dix, 1955)
who first recognized the connection between stacking velocities and rms velocities and showed how
to calculate an interval velocity from two stacking velocity measurements.
The application of equation (2.30) in practice requires the measurement of stacking velocities and
vertical traveltimes to two closely spaced reflectors. Under the assumption that stacking velocities
are well approximated by rms velocities (Dix (1955) or Taner and Koehler (1969)), the rms velocity
of the interval is then estimated with equation (2.30). However, as with averagevelocities, errors in
2 2
measurement lead to problems with the estimation of vrms (τ2 , τ1 ). If the term τ2 vrms2 − τ1 vrms1
becomes negative then imaginary interval velocity estimates result. Thus one essential condition for
the use of this technique is that
2 2 τ1
vrms2 > vrms1 . (2.31)
τ2
Since τ1 /τ2 < 1, this is a constraint upon how fast vrms estimates can decrease with increasing
time. There is no mathematical basis to constrain the rate at which vrms can increase; however, it
is reasonable to formulate a constraint on physical grounds. Since P-wave seismic velocities are not
expected to exceed some vmax (say 7000 m/s) a constraint would be
2 2 τ1 2 τ2 − τ1
vrms2 < vrms1 + vmax . (2.32)
τ2 τ2
Exercise 2.7.1. Show that the right-hand-side of inequality (2.32) is always greater than vrms1
provided that vmax > vrms1 so that this is a constraint on the rate of increase of vrms .
Exercise 2.7.2. Suppose that the times and depths to two reflectors are known, say τ1 , τ2 , z1 , and z2 ,
and that the rms velocities to the reflectors are also known, say vrms1 and vrms2 . Consider the interval
velocities defined by
z2 − z 1 2
vrms2 2
τ2 − vrms1 τ1
vint = and ṽint = .
τ2 − τ1 τ2 − τ1
Under what condition(s) will these interval velocity estimates be similar? If the interval between the
reflectors is highly heterogeneous, which estimate will be larger. Why?
0
0 0
0.1
0.1
0.2
500 0.2
0.3
0.3
seconds
seconds
0.4 b
meters
1000 0.4
0.5
0.5 a
0.6
1500 0.6 0.7
0.7 0.8
c d
2000 0.8 0.9
1000 2000 3000 4000 1000 2000 3000 4000 1000 1500 2000 2500 3000 3500 4000
meters/sec meters/sec meters/sec
Figure 2.9: (Left) A five legged vins (z) with lin- Figure 2.10: (a) vins (τ ) from Figure 2.9. (b)
ear variation between the knees. (Right) The finely sampled vrms (τ ). (c) de-sampled vrms (τ ).
curve on the left has been converted to vins (τ ). (d) ten legged interval rms approximation to
The variation between the knees is no longer lin- vins (τ ).
ear.
from one form to another is often necessary in data processing and software tools are required for
this purpose. Two alternate approaches are (1) fit the numerical velocity functions with an nth
order polynomial and perform the necessary integrations using polynomial integration rules, or (2)
implement a numerical integration scheme for the conversion formulae. The second approach is
taken here.
Five functions are available to convert velocity functions from one form to another. There is also
a utility plotting function to draw piecewise constant lines.
No utility is provided to deal with instantaneous velocity since this can be treated simply by
generating a finely sampled interval velocity. Like seismic traces, velocity functions are specified
by prescribing two vectors, one for the velocities and the other for the depth or time that they
apply. Interval velocities are inherently piecewise constant and the convention used here is that the
kth velocity, vk prescribed at depth zk , persists as a constant until depth zk+1 . Therefore, the last
velocity in a velocity vector applies for all z greater than the last entry in the depth vector.
As a first example of the use of these functions, suppose that vins (z) is a piecewise linear function
with five “knees” prescribed by (v, z) pairs as (1200 m/s, 0 m), (2700 m/s, 500 m), (3500 m/s, 1400
m), (3000 m/s, 1600 m), (4000 m/s, 2000 m). Code Snippet 2.8.1 shows how to compute vins (τ ) for
40 CHAPTER 2. VELOCITY
this vins (z) and displays its results in Figure 2.9. Line 2 uses piecewise linear interpolation (pwlint )
to resample the five legged function to a fine interval so that the piecewise constant interval velocity
approximates the desired piecewise linear instantaneous velocity. Line 5 computes τ (z) by calling
vint2t and thus defining vins (τ ). (Note that the vint2t returns one-way time so these values must
be doubled for two-way time.)
Code Snippet 2.8.1. This code defines a five legged, piecewise linear vins (z) and then computes
vins (τ ). It makes Figure 2.9.
1 v=[1200 2700 3500 3000 4000];z=[0 500 1400 1600 2000];
2 z2=0:10:2000;v2=pwlint(z,v,z2);
3 subplot(1,2,1);plot(v2,z2,v,z,’r*’);flipy
4 xlabel(’meters/sec’);ylabel(’meters’);
5 t2=vint2t(v2,z2);
6 t=interp1(z2,t2,z);
7 subplot(1,2,2);plot(v2,t2,v,t,’r*’);flipy;
8 xlabel(’meters/sec’);ylabel(’seconds’);
End Code
The computation of τ (z) requires the numerical computation of the integral in equation (2.3).
Given a densely sampled function, the simplest way to do a numerical integration is with the functions
sum and cumsum. The difference between these is that the first does a definite integral while the
second does an indefinite integral. For example, the command sum([1:5].ˆ2) just adds the squares
of the integers from 1 to 5 to get 55. This is a discrete approximation to the definite integral of the
function y = x2 from .5 to 5.5 for which the analytic answer is (5.53 −.53 )/3 ≈ 55.4167. Alternatively
cumsum([1:5].ˆ2) has the result [1, 5, 14, 30, 55] which gives the value of the integral of y = x2 for a
lower
5.5 2limit of .5 and five successive upper limits of 1.5, 2.5, 3.5, 4.5, xand2 5.5 . Thus sum approximates
.5
x dx and is just a single number while cumsum approximates .5 x̃ dx̃ and results in a vector the
same length as the input. Thus, the traveltime integration of equation (2.3) is done within vint2t
with the single command t1(2:nz)=cumsum(dz./vint(1:nz-1)). The vector t1 is τ (z) and has
been initialized to zero, dz is a vector of the depth intervals, and vint is the interval velocity vector.
(For more details, browse the source code file.)
Often velocity functions from well logs do not begin at z = 0. Therefore the computation of τ (z)
from z = 0 requires that some initial velocity be used for the depth range above the log. This can
be input to vint2t , in the form of a constant time shift τ0 , as the fourth argument. τ0 defaults
to τ0 = z(1)/v(1), that is the first depth divided by the first velocity. This can be estimated if the
depth and time to some marker formation top are already known. Then vint2t can be run once
with the initial velocity set to zero, and the time to the marker can be observed to be off by some
∆τ in one-way time.
Now consider the computation of vrms (τ ) and the construction of a ten layer approximation to
vins (τ ) using interval rms velocities. This is done in Code Snippet 2.8.2 and the result is shown in
Figure 2.10. Line 2 creates the dense vrms (τ ) function using vint2vrms shown as curve ‘b’ in the
figure. Line 3 defines a blocky ten legged time vector and line 4 creates the coarse vrms (τ ) (curve
‘c’) by calling vint2vrms with a third argument. Essentially, this is just a de-sampled version of
the finely sampled vrms (τ ). (It was not necessary to create the finely sampled vrms (τ ) as this was
done to demonstrate the software.) Finally line 6 creates the interval rms approximation to vins (τ )
by sending the coarse vrms (τ ) to vrms2vint .
Code Snippet 2.8.2. This carries on after Code Snippet 2.8.1 to compute: vrms from the sur-
face for every point in vins (τ ) a ten point sparse vrms (τ ) and a ten-legged rms-interval velocity
approximation to vins (τ ). Figure 2.10 is the result.
2.9. APPARENT VELOCITY: VX , VY , VZ 41
x x + ∆x Receivers
x x + ∆x Receivers
θ λx
θ
Receivers
tim
e
t
λ
+
w λz
∆
av
t
ef
ro
nt z - ∆z
Wa z - ∆z
tWt
m i
eW
tW θ
z
z
Figure 2.11: An impulsive plane wave ap- Figure 2.12: Similar to Figure 2.11 except that
proaches horizontal and vertical receiver arrays the plane wave is monochromatic.
1 plot(v2,t2,v,t,’r*’);flipy
2 vrms=vint2vrms(v2,t2);
3 tblock=linspace(min(t2),max(t2),10);
4 vrmsblock=vint2vrms(v2,t2,tblock);
5 drawvint(t2,vrms);drawvint(tblock,vrmsblock);
6 vint=vrms2vint(vrmsblock,tblock);
7 drawvint(tblock,vint);
8 xlabel(’meters/sec’);ylabel(’seconds’);
End Code
Exercise 2.8.1. Load the file vp_from_well.mat that contains a vector of P-wave velocities and a
corresponding vector of depths. Create two ten-leg interval-velocity approximations to vins (z) one
using average velocities and one using rms velocities. Compare the accuracy of time to depth con-
version using these two alternate approximations. (For time to depth conversions, compute average
velocities from both interval velocity functions). Is there any significant difference?
Exercise 2.8.2. Create a function called vrms2vave that converts average velocities to rms ve-
locities. Then create the inverse function vave2vrms. Test your codes with the universal velocity
function by comparing with the analytic forms for vrms and vave .
Exercise 2.8.3. Working with the universal velocity function, create vrms (τ ) for the depth range of
0 to 3000 m. Then use MATLAB’s random number generator, rand, to add a 2% random fluctuation
to the vrms (τ ) function. Finally calculate vins (τ ) from the perturbed vrms (τ ) function and determine
the percentage error in the vins (τ ) estimates caused by the 2% error in vrms (τ ). What does this say
about the sensitivity of interval velocity calculations to noise?
immediately beneath the receivers has the acoustic velocity, v0 , then the wavefront arrives at x + ∆x
later than it does at x by the delay
sin θ
∆t = ∆x. (2.33)
v0
Thus it appears to move along the array at the speed
∆x v0
vx ≡ = . (2.34)
∆t sin θ
The quantity vx is called an apparent velocity because the wavefront appears to move at that speed
even though its true speed is v0 . Apparent velocities are one of the four fundamental observables in
seismic data, the others being position, time and amplitude. The apparent velocity is never less than
the real velocity and can range up to infinity. That is v0 ≤ vx ≤ ∞. Since infinities are cumbersome
to deal with it is common to work with the inverse of vx , called the time-dip, ∆t/∆x or horizontal
slowness, sx . Another common convention for horizontal slowness is to call it p which signifies the
ray parameter.
Now, enlarge the thought experiment to include an array of receivers in a vertical borehole and
spaced at ∆z. Reasoning similar to that used before shows that the arrival at z − ∆z is delayed
from that at z by
cos θ
∆t = ∆z. (2.35)
v0
Therefore, the vertical apparent velocity, vz , is
∆z v0
vz ≡ = . (2.36)
∆t cos θ
Using simple trigonometry, it results that
1 1 1
= 2 + 2. (2.37)
v02 vx vz
If this argument is extended to 3D, the result is that there is an apparent velocity in the y direction
as well and that the sum of the inverse squares of the apparent velocities is the inverse square of the
real velocity:
1 1 1 1
2 = 2 + 2 + 2. (2.38)
v0 vx vy vz
We observe apparent velocity, or equivalently time-dip, by simply choosing an event of interest on
a seismic record and measuring the slope, ∆t/∆x. It will be seen shortly that this is a measurement
of the ray parameter of the ray associated with the chosen event. Knowledge of the ray parameter
essentially determines the raypath, provided that the velocities in the subsurface are known. This
allows the ray to be projected down into the earth and to possibly determine it’s reflection point.
This technique, called raytrace migration, will be discussed in section 4.2.2.
Another way to measure apparent velocities is with a multi-dimensional Fourier transform. In
2-D for example, the f-k transform represents a seismic record on a grid of horizontal wavenumber,
kx , and temporal frequency, f . Each point in the (kx ,f ) plane has a complex number associated
with it that gives the amplitude and phase of a fundamental Fourier component: e2πi(kx x−f t) . In
3D these fundamental components are monochromatic (i.e. a single f ) plane waves so in 3D or 2D
they are called Fourier plane waves. These waves are infinite in extent so they have no specific
arrival time. However, it does make sense to ask when a particular wave crest arrives at a particular
2.9. APPARENT VELOCITY: VX , VY , VZ 43
location. Mathematically, a wave crest is identified as a point of constant phase where phase refers
to the entire argument of the complex exponential. If the Fourier transform has the value Aeiφ at
some (kx ,f ), then the Fourier plane wave at that point is Ae2πi(kx x−f t)+iφ . The motion of a point
of constant phase is tracked by equating the phase at (x, t) with that at (x + ∆x, t + ∆t) as in
If equation (2.38) is evaluated for the Fourier plane wave e2πi(kx x+ky y+kz z−f t) the result is
1 k2 ky2 k2
2
= x2 + 2 + z2 . (2.41)
v f f f
Now, for any monochromatic wave, we have λf = v and, using k = λ−1 , then equation (2.41) leads
to
k 2 = kx2 + ky2 + kz2 . (2.42)
This very important result shows that the coordinate wavenumbers, (kx ,ky ,kz ), are the components
of a wavenumber vector, #k. Using a position vector, #r = (x, y, z), the Fourier plane wave can be
written compactly as e2πi(k·r−f t) . Equation (2.42) can be expressed in terms of apparent wavelengths
(e.g. kx = λ−1
x , etc.) as
1 1 1 1
2
= 2 + 2 + 2. (2.43)
λ λx λy λz
which shows that, like apparent velocities, the apparent wavelengths are always greater than the
true wavelength. An important property of the wavenumber vector is that it points in the direction
of wave propagation (for an isotropic, homogeneous medium). An easy way to see this is to take the
gradient of the phase of a Fourier plane wave. Since the wavefront is defined as a surface of constant
phase, then wavefronts can be visualized as contours of the phase function φ̃ = πi(#k · #r − f t). The
gradient of this phase function points in the direction normal to the wavefronts, that is the direction
of a raypath. This gradient is
# φ̃ = ∂ φ̃ x̂ + ∂ φ̃ ŷ + ∂ φ̃ ẑ
∇ (2.44)
∂x ∂y ∂z
where x̂, ŷ, and ẑ are unit vectors in the coordinate directions. Calculating the indicated partial
derivatives leads to
∇# φ̃ = kx x̂ + ky ŷ + kz ẑ = #k. (2.45)
This equation can be achieved more simply by writing ∇ # = r̂ ∂ , where r = (x2 + y 2 + z 2 ) and r̂
∂r
is a unit vector pointing to a particular location, so that
#
# φ̃ = r̂ ∂ k · #r = r̂ #k = #k.
∇ (2.46)
∂r
44 CHAPTER 2. VELOCITY
z z
δ x
ŷyẑz
+
x
x̂
∂
=
φ̃
∇
#
x
medium 1 θ
λ1
θ ϕ ϕ
δ medium 2
λ2
Figure 2.13: Snell’s law results from the physical requirement that the apparent velocity
along an interface be conserved.
So, the inverse apparent velocities are proportional to the components of the wavenumber vector
that points in the direction of wave propagation.
z
θ v
φ
z
φ
v
θ
z
v φ θ
z
v φ
Figure 2.14: In a v(z) medium, all velocity interfaces are horizontal and ray propagation
is especially simple.
transmitted energy. In an acoustic medium, where there is only one mode of wave propagation, this
results in the angle of incidence and the angle of reflection being equal. However, in elastic media,
the incident and reflected waves can be either P-waves or S-waves. For example, in the case of an
incident P-wave reflecting as an S-wave, Snell’s law states
sin θp sin θs
= . (2.48)
vp vs
As a final summary statement, Snell’s law states that wavefront propagation across a velocity inter-
face always conserves the apparent velocity along the interface. Though it is not proven here, this
holds for arbitrary wave types and for non-planar wavefronts and interfaces.
This analysis may be repeated at any other interface with a similar conclusion. Thus the quantity
p = sin θj /vj is conserved throughout the entire wave propagation. p is generally referred to as the
ray parameter since it is a unique constant for any ray and so parameterizes (or names) the possible
rays. The conservation of p and its identification as the horizontal apparent slowness is a direct
consequence of Snell’s law. This analysis generalizes to a continuous variation of v with z such that
sin (θ(z))
p≡ = a constant for any particular ray. (2.51)
v(z)
46 CHAPTER 2. VELOCITY
dx
dz
ds
θ(z)
Figure 2.15: A differential ray element ds, at depth z, and traveling at an angle θ(z) with
respect to the vertical.
General expressions for the traveltime and horizontal distance traveled by any ray in a v(z)
medium can be easily derived. Figure 2.15 shows a differential element of a ray. From the geometry,
it follows that
dx = tan (θ(z)) dz (2.52)
and
ds dz
dt = = . (2.53)
v(z) v(z) cos (θ(z))
Snell’s law is incorporated by replacing the trigonometric functions using pv(z) = sin (θ(z)) so that
pv(z)
dx = dz (2.54)
1 − p2 v 2 (z)
and
dz
dt = . (2.55)
v(z) 1 − p2 v 2 (z)
Expressions for macroscopic raypaths are obtained by simply integrating these results. For a ray
traveling between depths z1 and z2 , the horizontal distance traveled becomes
z2
pv(z)
x(p) = dz (2.56)
z1 1 − p2 v 2 (z)
If the v(z) medium is discretely layered rather than continuously variable, then summation forms
2.11. RAYTRACING IN A V (Z) MEDIUM 47
s Receivers r
ZOOM
∆r θ
θ θ
EmergentWRay ∆s
Figure 2.16: The ray parameter can be measured directly by measuring the delay between
wavefront arrivals at successive receivers.
of equations (2.56) and (2.57) are more appropriate. These are easily seen to be
n
pvk
x(p) = ∆zk (2.58)
k=1
1 − p2 vk2
and
n
∆z
t(p) = k . (2.59)
v 1 − p2 vk2
k=1 k
∆t 1 sin θ0
= = = p. (2.61)
∆r vr v0
where vr is the horizontal apparent velocity of the wave at the measurement location.
Generally, it is expected that p changes with position due to changes in emergence angle and due
to lateral variations in the instantaneous velocity. Horizontal events, i.e. waves traveling vertically
when they reach the receiver array, have a ray parameter of 0. Waves traveling horizontally across
the array have a ray parameter of 1/v0 and on a seismic (x, t) plot have the maximum possible
48 CHAPTER 2. VELOCITY
A B kx
x
t
Figure 2.17: The physical considerations that place limits upon ray parameter manifest as
limits on the maximum time-dip in (x, t) space (A) and segment (kx , f ) space into allowed
and forbidden regions (B).
slope. Even though the wavefronts are vertical, the slope on the time section cannot exceed 1/v0 .
Taking sign into account, the range of possible values for the ray parameter is −v0−1 ≤ p ≥ v0−1 .
This maximum possible steepness in (x, t) space is a fundamental property of a seismic time section.
Assuming a value for v0 in the course of data processing means that any events with slopes steeper
than v0−1 are interpreted as non-physical and must be rejected. Since apparent velocities can also
be expressed in (kx , f ) space as f /kx (equation (2.40)), this phenomenon means that (kx , f ) space
is segmented into allowed and forbidden regions.
Figure 2.17 shows this situation for both (x, t) and (kx , f ) spaces. In (x, t) space, the fan of
allowed time dips forms a bow-tie when plotted at a specific point, though such events may exist
all (x, t) locations. In (kx , f ) space, the zero ray parameter plots vertically down the f axis while
it is horizontal in (x, t). Thus the fan of allowed p values in (kx , f ) forms a symmetric shape about
the f axis. (Only one half of the bow tie is shown here because negative frequencies are generally
redundant for real-valued data.) In (kx , f ) space, p values are found only along radial lines emanating
from the origin, not along lines of the same slope emanating from some other point. The portion
of (kx , f ) space corresponding to |kx /f | < v0−1 (i.e. outside the shaded region in Figure 2.17B) is
“forbidden” to simple rays and is called the evanescent region. (It will be seen later that certain
exponentially decaying wave types can occur here.) The shaded portion of Figure 2.17B is called
the body wave region and is the portion of (kx , f ) space that is available for seismic imaging.
and
1 v0 + cz 1 + 1 − p2 v02 .
t(z, p) = ln (2.63)
c v0 2
1 + 1 − p2 {v0 + cz}
Slotnick shows that equation (2.62) describes an arc of a circle, having radius 1/(pc) and centered
at x0 = 1 − p2 v02 /(pc) and z0 = −v0 /c. Manifestly, equation (2.62) describes a ray as it goes from
zero to some depth z. Since velocity is always increasing in this case, the Snell’s law angles are always
getting larger as a ray goes deeper. Eventually the ray flattens out, when θ(z) = sin−1 (pv(z)) = 90◦ ,
2.11. RAYTRACING IN A V (Z) MEDIUM 49
0 0
5 5
10
10
z kilometers
z kilometers
15
15
20
20
25
30 25
35 30
0 5 10 15 20 25 30 35 0 5 10 15 20 25 30 35
x kilometers x kilometers
Figure 2.18: A selection of raypaths are shown Figure 2.19: A set of wavefronts associated with
for the universal velocity function of Figure 2.1. the raypaths of Figure 2.18. This was created
Code Snippet 2.11.1 created this plot. with Code Snippet 2.11.2.
and turns upward. Therefore, the depth at which a ray bottoms out, called its turning point, is found
from
1 − pv0
zmax = . (2.64)
pc
The complete raypath for a ray that reaches its turning point and then returns to the surface
must have two values of x for each z so the function x(z) is mathematically double valued. However,
z(x) is still single valued so the complete raypath is most easily computed by solving equation (2.62)
for z to get
1 2
z= 1 − {pcx − cos θ0 } − pv0 . (2.65)
pc
Code Snippet 2.11.1 implements equation (2.65) to create the raypath display shown in Figure 2.18.
This code establishes a horizontal distance range and a set of takeoff angles on lines 1-3. Then, in
looping over the desired rays, it calculates a depth for each element of the vector x. However, many
of these rays will not cover the full range of x. Equation (2.65) returns a complex number for an x
distance that cannot be reached by a particular ray. Lines 11-14 search for these points and replace
their depths with NaN. (Recall that NaN’s do not display when plotted.) Also, the code generates
z values that are negative so lines 15-18 set these to NaN. The ray bending in Figure 2.18 is much
more realistic than a simple constant velocity calculation using straight rays. It shows that most
of the energy of a seismic source cannot penetrate to great depths because it is turned around by
refraction. All of the raypaths are obviously circular arcs whose centers move to the right as the ray
parameter decreases.
Code Snippet 2.11.1. This code implements equation (2.65) and makes Figure 2.18.
1 x=1:50:35000;
2 vo=1800;c=.6;nrays=10;
3 thetamin=5;thetamax=80;
4 deltheta=(thetamax-thetamin)/nrays;
5 zraypath=zeros(nrays,length(x));
6 for k=1:nrays
50 CHAPTER 2. VELOCITY
7 theta=thetamin+(k-1)*deltheta;
8 p=sin(pi*theta/180)/vo;
9 cs=cos(pi*theta/180);
10 z = (sqrt( 1/(pˆ2*cˆ2) - (x-cs/(p*c)).ˆ2) -vo/c);
11 ind=find(imag(z)˜=0.0);
12 if(˜isempty(ind))
13 z(ind)=nan*ones(size(ind));
14 end
15 ind=find(real(z)<0.);
16 if(˜isempty(ind))
17 z(ind)=nan*ones(size(ind));
18 end
19 zraypath(k,:) = real(z);
20 end
21 figure;plot(x/1000,zraypath/1000);flipy;
22 xlabel(’x kilometers’);ylabel(’z kilometers’)
End Code
Slotnick also derives expressions for the wavefronts (surfaces of constant traveltime) and shows
them to be circles whose centers are along the z axis at
v0
zw0 = [cosh(cT ) − 1] (2.66)
c
where T is the traveltime defining the wavefront. Each circular wavefront has a radius of
v0
r= sinh(cT ). (2.67)
c
Code Snippet 2.11.2 calculates and plots 10 wavefronts for a set of traveltimes from 0 to 5 seconds.
For each wavefront, the strategy is to calculate the center of the wavefront circle (line z) and its
radius (line 8). Then, for a predefined set of depths (line 1) the horizontal positions are calculated
from the equation of a circle (line 9). As in the case of the raypaths, this results in both complex
and negative values which are found and set to NaN (lines 10-17). The resulting wavefronts clearly
show the effects of increasing velocity with depth, being more closely spaced at shallow depths than
when deeper. Figure 2.20 superimposes the raypaths and wavefronts and uses the command axis
equal to ensure a 1:1 aspect ratio. Clearly the raypaths are normal to the wavefronts.
Code Snippet 2.11.2. This code assumes that Code Snippet 2.11.1 has already been run and
proceeds from there to calculate and plot the wavefronts. Figure 2.19 is the result.
1 zw=0:50:30000;
2 nwaves=10;tmax=5;
3 xwavefront=zeros(nwaves,length(zw));
4 times=linspace(0,tmax,nwaves);
5 zo=zeros(1,nwaves);
6 for k=1:nwaves
7 zo(k)=vo*(cosh(c*times(k))-1)/c;
8 r=vo*sinh(c*times(k))/c;%radius
9 xw=sqrt(r.ˆ2-(zw-zo(k)).ˆ2);
10 ind=find(real(xw)<0.);
11 if(˜isempty(ind))
12 xw(ind)=nan*ones(size(ind));
13 end
14 ind=find(imag(xw)˜=0.0);
15 if(˜isempty(ind))
16 xw(ind)=nan*ones(size(ind));
2.11. RAYTRACING IN A V (Z) MEDIUM 51
10
z kilometers
15
20
25
30
0 5 10 15 20 25 30
x kilometers
Figure 2.20: This figure superimposes the raypaths of Figure 2.18 onto the wavefronts of
Figure 2.19.
17 end
18 xwavefront(k,:) = real(xw);
19 end
20 figure;plot(xwavefront/1000,zw/1000);flipy;
21 xlabel(’x kilometers’);ylabel(’z kilometers’)
End Code
Exercise 2.11.1. Use equations (2.62) and (2.63) to calculate and plot the (x, t) curve for a P-P
reflection from 2000 meters depth. Assume the universal velocity function and a source and receivers
at depth 0. Repeat your calculations for a P-S reflection assuming a constant vp /vs of 2. In each
case, what is the maximum source-receiver offset for which a reflection is expected?
Exercise 2.11.2. Derive equations (2.62) and (2.63) from equations (2.56) and (2.57).
Exercise 2.11.3. For the linear increase of velocity with depth, a source at (x, z) = (0, 0), and
considering transmitted rays only, is there a unique p for each point in the (x, z) plane? Will your
conclusion remain valid for more general v(z)?
z = 0 to z = zr at offset x/2 and then double the results. However, equation (2.56) does not easily
facilitate this since the value of p that will do the job is not known. If velocity were constant, then
simple geometry would predict the takeoff angle of θ0 = arctan(x/(2zr )) and thus the ray parameter
is trivially found. For increasing velocity with depth, the takeoff angle will be smaller than that
predicted with constant velocity while for decreasing velocity with depth the situation is reversed.
In the general case, the takeoff angle and hence the ray parameter and cannot be found analytically.
There are seven functions and one demonstration script provided for v(z) raytracing. These are
rayfan Shoots a fan of rays given their ray parameters for v(z).
The approach taken here is to shoot a fan of rays and to iteratively refine the fan until a
convergence criterion is met. For a general case, this iteration will shoot many fans of rays so the
ray fan algorithm needs to be very efficient. The functions rayfan , rayfan a , and shootray all
performs this task in s similar fashion. However, only the latter is actually used in the two-point
raytracers, traceray pp , traceray ps , and traceray , because it is the most efficient. However,
the efficiency of shootray is achieved at the sacrifice of some flexibility so rayfan and rayfan a
are provided as well.
As shown in Code Snippet 2.11.3, shootray has three inputs: v, z, and p that are (respectively)
column vectors of velocity and depth and a row vector of ray parameters. The vector of rays defined
by p is traced from z(1) to z(end) using the velocities in v. As discussed in section 2.8, the velocity
model is assumed to be piecewise constant, with velocity v(k) beginning at z(k) and persisting as
a constant until z(k+1). Accordingly, the last velocity in v is irrelevant for shootray , and line 3
establishes an index vector to the relevant velocities. Lines 4 through 8 are key to the efficiency
of shootray . On line 4, v(iprop) is an m length column vector and p is an n length row vector.
This product of an m × 1 matrix representing v(z) with an 1 × n matrix representing p is an m × n
matrix of sin θ = pv(z) called sn. The kth column of sn represents the product pv(z) for the kth
ray parameter. Line 6 builds the matrix cs that represents cos θ(z) = 1 − p2 v 2 (z). Lines 7 and
8 build the m × n matrices vprop and thk that represent v(z) and ∆z for each ray. Since these
quantities are the same for each ray, they are represented by matrices with identical columns. Thus,
we have four m × n matrices, with z as the row coordinate and p as the column coordinate. Not all
of these combinations of z and p can correspond to physical rays since the depth of penetration of
a ray is limited. Line 5 detects these non-physical combinations by finding any values of pv(z) that
are greater than one.
Code Snippet 2.11.3. The function shootray shoots a fan of rays as defined by the row vector p
from z(1) to z(end) through the velocity model defined by the column vectors v and z.
2.11. RAYTRACING IN A V (Z) MEDIUM 53
1 function [x,t]=shootray(v,z,p)
2 ... code not displayed ... see shootray.m
3 iprop=1:length(z)-1;
4 sn = v(iprop)*p;
5 [ichk,pchk]=find(sn>1);
6 cs=sqrt(1-sn.*sn);
7 vprop=v(iprop)*ones(1,length(p));
8 thk=abs(diff(z))*ones(1,length(p));
9 if(size(sn,1)>1)
10 x=sum( (thk.*sn)./cs);
11 t=sum(thk./(vprop.*cs));
12 else
13 x=(thk.*sn)./cs;
14 t=thk./(vprop.*cs);
15 end
16 %assign infs
17 if(˜isempty(ichk))
18 x(pchk)=inf*ones(size(pchk));
19 t(pchk)=inf*ones(size(pchk));
20 end
End Code
The ray tracing is completed in lines 9 through 15 and an if statement is used to handle the
single layer case. Using the m × n matrices thk, sn, cs, line 10 computes equation (2.58) by using
the .* operator to form the expression pvk ∆zk / 1 − p2 vk2 as an m × n matrix. Then, as discussed
in section 2.8, the function sum is used to compute the discrete sum that approximates the integral
in equation (2.56). When applied to a matrix, sum produces a row vector that is the sum of each
column of the matrix. This behavior dictated the construction of the m × n matrices with p constant
in each column. The if statement is required for the single layer case also because of the behavior
of sum. If there is only one layer, then the m × n matrices are all 1 × n row vectors. When given
a row vector, sum adds its entries to get a single value instead of “summing” the columns of length
1. The if statement circumvents this behavior by omitting the sum entirely in the single layer case.
The final step, lines 17 through 20 assigns Inf to those non-physical values that were flagged earlier
on line 5.
The function shootray does not check for consistency of the inputs to ensure that v and z are
column vectors and p is a row vector. Also, it lacks flexibility to specify the start and end depths
and always shoots the rays from z(1) to z(end). rayfan traces rays with the same scheme as
shootray but incorporates error checking and allows the start and end depths to specified. This
makes rayfan easier to use but slower. shootray is intended to be used within a larger raytracing
scheme while rayfan is more easily used from the command line.
The two-point ray tracing functions traceray pp , traceray ps , and traceray all work simi-
larly with an iterative scheme involving multiple calls to shootray . The codes are quite intricate
and will not be displayed. traceray pp and traceray ps are constructed to trace P-P and P-S
primary reflections and nothing else. traceray is more general and can trace an arbitrary ray
that has multiple bounces (reflections) and mode conversions. All three codes use similar schemes
of building an equivalent layered model that allows the problem to be addressed by shooting rays
in one direction only. For example, the P-P reflection problem requires a ray to be traced down
through a set of layers to a specified depth and then back up again though nearly the same set of
layers. (It will be exactly the same set if the source and receiver depths are the same.) Instead of
this two-way problem, a one-way problem is set up by determining the layers for the upgoing leg
and placing them beneath the reflector in the order that they are encountered (Figure 2.21). The
layers containing the source, receiver, and reflector are adjusted in thickness so that shootray can
54 CHAPTER 2. VELOCITY
A B
v source v source
v receiver v
v v
v v
v v
v v
v v
v v
v v
v v
v
v
v
v
v
v
v
v
receiver
Figure 2.21: (A) A P-P reflection between a source and receiver at different depths. (B)
The equivalent one-way raypath to the two-way ray shown in (A). The dashed line is the
reflector in both cases.
be used to trace rays from z(1) to z(end). A similar equivalent layering can always be constructed
for a ray that changes mode (between P and S) and that takes any number of extra bounces. In the
former case, the equivalent layering simply must use the correct velocities. In the latter case, the
equivalent layering can contain many repeats of a section. This is a technically simple scheme for
computing any ray in v(z).
The two-point raytracing functions can all trace rays from a single starting point to multiple
receivers at once. That is, they are vectorized to produce simple source gathers. Though the source
depth need not equal the receiver depth, if multiple receivers are specified, they must all be at the
same depth. This means that geometries like that of a VSP must be handled by calling the ray
tracer separately for each receiver depth. All three programs iterate until all receivers have a ray
traced within a capture radius of their position or until a maximum number of iterations is reached.
A capture radius is the radius of an imaginary circle around each receiver which, if a ray falls within
it, is considered “good enough” for that receiver. If the flag optflag is set to 1, then the actual
traveltime and ray parameter returned are linearly interpolated between the captured ray and the
next closest ray. If optflag is zero, then the captured ray is used directly.
Unlike many raytracers, the codes here can reflect, or mode convert, a ray at any depth not just
at a layer boundary. Though this may not be physically correct, it allows the common (and very
practical) practice of separating the reflectivity from the propagation model. The idea is that the
rays are traced through a simplified background medium which is chosen to give sufficiently accurate
traveltimes for the task at hand. Reflections are assumed to come from a level of detail beyond that
required for the traveltime calculations. For example, a linear variation of v(z) may be used with
locally determined constants as a background velocity function. Of course, the user who disagrees
with this approach is free to prescribe reflections at actual velocity layer boundaries.
Code Snippet 2.11.4 exhibits the use of traceray pp and traceray ps to model P-P and P-S
reflections with the universal velocity function as the background medium. The results are shown
in Figures 2.22 and 2.23. Line 1 defines the vp and vs velocity models with a vp /vs of 2. Lines 2
and 3 establish the basic geometry of source, receivers, and reflector. The capture radius is set to
10% of the receiver spacing and the maximum number of iterations is set to 4 which is adequate for
smooth media. The call to traceray pp is on line 7 and the final parameter dflag instructs the
2.11. RAYTRACING IN A V (Z) MEDIUM 55
Vertical gradient simulation, P−P mode zsrc=100 zrec=500 Vertical gradient simulation, P−S mode zsrc=100 zrec=500
0 0
1000 1000
meters
meters
2000 2000
3000 3000
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
meters meters
2 2.8
seconds
seconds
2.1 3
2.2 3.2
2.3 3.4
1000 1500 2000 2500 3000 1000 1500 2000 2500 3000
meters meters
Figure 2.22: A P-P reflection is shown for a Figure 2.23: A P-S reflection is shown for a back-
background medium of vp (z) = 1800 + .6zm/s. ground medium of vs (z) = 900 + .3z. See Code
See Code Snippet 2.11.4. Snippet 2.11.4
program to draw the rays into the current figure window. Lines 9 and 10 draw symbols indicating
the source and receiver and line 12 plots the traveltime in the bottom half of the figure. Lines 15-22
are very similar except that traceray ps is called to create a P-S reflection. An S-P reflection could
be created by reversing the order of the velocity arguments (i.e.traceray ps(vs,zs,vp,zp,...)).
Similarly, an S-S reflection can be modelled using traceray pp and providing it with the S-wave
velocity functions.
Code Snippet 2.11.4. This example raytraces a P-P and a P-S reflection in a linear gradient
medium. It creates Figures 2.22 and 2.23.
VSP Vertical gradient simulation, P−P mode VSP Vertical gradient simulation, P−S mode
0 0
1000 1000
meters
meters
2000 2000
3000 3000
0 500 1000 1500 0 500 1000 1500
meters meters
500 500
depth (meters)
depth (meters)
1000 1000
1500 1500
2000 2000
2500 2500
1.4 1.6 1.8 2 2.2 1.5 2 2.5 3
seconds seconds
Figure 2.24: A P-P reflection is shown for an Figure 2.25: A P-S reflection for an offset VSP
offset VSP. See Code Snippet 2.11.5. is shown. See Code Snippet 2.11.5.
25 subplot(2,1,2);plot(xoff,t);grid;
26 xlabel(’meters’);ylabel(’seconds’);flipy;
End Code
Often very interesting calculations can be done by calling these functions in a loop. For example,
since the receivers are required to be at a constant depth on each call to traceray pp a VSP cannot
be modelled with a single call. However, as Code Snippet 2.11.5 shows, the desired results can be
obtained by looping over receiver depth. Though not as efficient in this mode, the calculation is still
simple and accurate. The parameter pfan is set to -2 in this case which causes the ray tracer to
start with the final ray fan determined the last time it was called. The modelled events are shown
in Figures 2.24 and 2.25.
Code Snippet 2.11.5. This code traces P-P and P-S reflections for an offset VSP by looping over
receiver depth. It creates Figures 2.24 and 2.25.
1 zp=0:10:4000;vp=1800+.6*zp;vs=.5*vp;zs=zp;%velocity model
2 zrec=500:100:2500;zsrc=0;zd=3000;xoff=1500;%geometry
3 caprad=10;itermax=4;%cap radius, and max iter
4 pfan=-2;optflag=1;pflag=1;dflag=2;% default ray fan, and various flags
5 figure;subplot(2,1,1);flipy
6 t=zeros(size(zrec));
7 for kk=1:length(zrec);
8 if(kk==1)dflag=-gcf;else;dflag=2;end
9 [t(kk),p]=traceray_pp(vp,zp,zsrc,zrec(kk),zd,xoff,caprad,pfan,...
10 itermax,optflag,pflag,dflag);
11 end
12 title([’ VSP Vertical gradient simulation, P-P mode ’])
13 line(xoff,zrec,’color’,’b’,’linestyle’,’none’,’marker’,’v’)
14 line(0,zsrc,’color’,’r’,’linestyle’,’none’,’marker’,’*’)
15 grid;xlabel(’meters’);ylabel(’meters’);
16 subplot(2,1,2);plot(t,zrec);
17 xlabel(’seconds’);ylabel(’depth (meters)’);grid;flipy;
18
19 figure;subplot(2,1,1);flipy;
2.11. RAYTRACING IN A V (Z) MEDIUM 57
A P−P−P−P−P−P−P−P−P−P mode in vertical gradient media A P−S−S−S−S−P−P−P−P−S mode in vertical gradient media
0 0
500 500
1000 1000
meters
meters
1500 1500
2000 2000
2500 2500
3000 3000
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
meters meters
3.2 4.65
3.25 4.7
3.3 4.75
time
time
3.35 4.8
3.4 4.85
3.45 4.9
1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000 1000 1200 1400 1600 1800 2000 2200 2400 2600 2800 3000
offset offset
Figure 2.26: A complicated multiple is shown Figure 2.27: A complicated ray is shown that
that remains a P-ray throughout. See Code converts back and forth from P to S as it
Snippet 2.11.6 bounces. See Code Snippet 2.11.6
20 t=zeros(size(zrec));
21 for kk=1:length(zrec);
22 if(kk==1)dflag=-gcf;else;dflag=2;end
23 [t(kk),p]=traceray_ps(vp,zp,vs,zs,zsrc,zrec(kk),zd,xoff,caprad,pfan,...
24 itermax,optflag,pflag,dflag);
25 end
26 title([’ VSP Vertical gradient simulation, P-S mode ’])
27 line(xoff,zrec,’color’,’b’,’linestyle’,’none’,’marker’,’v’)
28 line(0,zsrc,’color’,’r’,’linestyle’,’none’,’marker’,’*’)
29 grid;xlabel(’meters’);ylabel(’meters’);
30 subplot(2,1,2);plot(t,zrec);
31 xlabel(’seconds’);ylabel(’depth (meters)’);grid;flipy;
End Code
As a final example of these raytracing facilities, consider the calculation of a complicated mode
that takes multiple bounces and converts back and forth between P and S modes. This calculation
can be done with traceray and is shown in Code Snippet 2.11.6. The key idea here is to use a
ray code that defines the depths and mode types of the different legs of the raypath. A ray code is
an n × 2 matrix with the first column being a list of depths and the second containing either the
integer 1 (P-wave) or 2 (S-wave). For example, the ray code [0 1;1000 2;800 2;1000 1;100 1]
specifies a P-ray from depth 0 to 1000 m. At 1000 m it converts to an S-ray and travels back up to
800 m and then back down to 1000 m while remaining an S-ray. On the second reflection at 1000
m it converts back to a P-ray and travels up to 100 m where it ends. The final entry in column two
is meaningless. Code Snippet 2.11.6 demonstrates this facility by creating a complicated multiple
bounce P-ray (no mode conversions) on line 8 and a multimode ray on line 17. Generally, it might
be expected that significant rays like this will be symmetric (i.e. the same bounce pattern on both
the up and down legs); however, creating an asymmetric ray is perhaps a better demonstration.
There are many more possible ways to use these raytracing facilities. For more examples, execute
the script raytrace demo and follow the on-screen instructions. The code of this script is similar
to the examples presented here and should be comprehensible.
58 CHAPTER 2. VELOCITY
Code Snippet 2.11.6. Here is a demonstration of the use of traceray to compute multiple-bounce
and multiple-mode raypaths. This code creates Figures 2.26 and 2.27.
1 zp=0:10:4000;vp=1800+.6*zp;vs=.5*vp;zs=zp;
2 xoff=1000:100:3000;
3 caprad=10;itermax=4;%cap radius, and max iter
4 pfan=-1;optflag=1;pflag=1;dflag=2;% default ray fan, and various
5
6 raycode=[0 1;1500 1;1300 1;2000 1;1800 1;3000 1;2000 1;2300 1;1000 1; 1500 1; 0 1];
7 figure;subplot(2,1,1);flipy
8 [t,p]=traceray(vp,zp,vs,zs,raycode,xoff,caprad,pfan,itermax,optflag,pflag,dflag);
9 title(’A P-P-P-P-P-P-P-P-P-P mode in vertical gradient media’);
10 xlabel(’meters’);ylabel(’meters’)
11 line(xoff,zeros(size(xoff)),’color’,’b’,’linestyle’,’none’,’marker’,’v’)
12 line(0,0,’color’,’r’,’linestyle’,’none’,’marker’,’*’);grid
13 subplot(2,1,2);plot(xoff,t);
14 grid;flipy;xlabel(’offset’);ylabel(’time’)
15
16 raycode=[0 1;1500 2;1300 2;2000 2;1800 2;3000 1;2000 1;2300 1;1000 1; 1500 2; 0 1];
17 figure;subplot(2,1,1);flipy
18 [t,p]=traceray(vp,zp,vs,zs,raycode,xoff,caprad,pfan,itermax,optflag,pflag,dflag);
19 title(’A P-S-S-S-S-P-P-P-P-S mode in vertical gradient media’);
20 xlabel(’meters’);ylabel(’meters’)
21 line(xoff,zeros(size(xoff)),’color’,’b’,’linestyle’,’none’,’marker’,’v’)
22 line(0,0,’color’,’r’,’linestyle’,’none’,’marker’,’*’);grid
23 subplot(2,1,2);plot(xoff,t);
24 grid;flipy;xlabel(’offset’);ylabel(’time’)
End Code
and vp /vs of 2. Create a display showing the P-P and P-S raypaths for source-receiver offsets from
10 to 3010 at 100 meter intervals for a reflector at a depth of 2500 m. Let the source and receivers
be a z = 0 Make another plot showing the angle of incidence on the reflector (in degrees) versus
offset for both P-P and P-S reflections.
Exercise 2.11.5. For the velocity model in the previous exercise, what is the largest offset that can
be successfully raytraced for the reflector at z = 2500m for both P-P and P-S reflections. Explain
why the ray tracing fails. How do these conclusions change if the reflector is at 1500 m?
calculations at the interfaces. The first approach is generally simpler because the raytracing can
be reduced to the solution of an ordinary differential equation that implicitly incorporates Snell’s
law. In the second approach, the description of the geologic interfaces must be done with great
care to ensure that they connect without overlap or voids and this can require complex geometric
calculations. Furthermore, it is generally a complex matter to determine the next interface that
a ray will encounter. However, for three dimensional simulations, the gridded model approach can
rapidly consume a great deal of computer memory can rapidly become impractical. This is especially
true if rays are to be traced through a complex structure that lies far beneath a simple one. In this
case, the grid must be made quite small for the complex structure which results in too much detail
for the upper medium. In this book, only the gridded approach in two dimensions will be explored.
1 ∂2ψ
∇2 ψ − =0 (2.68)
v 2 (#x) ∂t2
where A(#x) and T (#x) are unknown functions describing amplitude and traveltime that are expected
to vary with position. In the constant velocity limit, A becomes constant while T still varies rapidly
which leads to the expectation that variation in A will often be negligible compared with variation
in T . Substitution of equation (2.69) into equation (2.68) will require computation of the Laplacian
of the assumed solution. This is done as
# ·∇
∇2 ψ(#x, t) = ∇ # · e2piif (t−T ) ∇A
# A(#x)e2πif (t−T (x)) = ∇ # − 2πiAe2πif (t−T ) ∇T
# (2.70)
Then, using this result and ∂t2 ψ(#x, t) = −4π 2 f 2 Ae2πif (t−T (x)) in equation (2.68) and equating real
and imaginary parts gives the two equations
2 ∇2 A 1
#
∇T − − =0 (2.72)
4π 2 f 2 A v(#x)2
and
A 2 # · ∇T
# = 0.
∇ T − ∇A (2.73)
2
So far, these results are exact and no simplification has been achieved. However, the second term
in equation (2.72) is expected to be negligible when velocity gradients are weak or when frequencies
are high regardless of the velocity gradient. When this term is discarded, the result is the nonlinear,
60 CHAPTER 2. VELOCITY
T
T
dx P
P
Figure 2.28: A differential raypath element d#x extends from a point P1 on a wavefront at
time T1 to a point P2 on a wavefront at time T2 . The ray segment is perpendicular to
both wavefronts and has length |d#x| = ∆s.
Though not exact, for many realistic situations, a solution to the eikonal equation gives accurate
traveltimes through complex media. Equation (2.73) is called the geometrical spreading equation
because its solution can be shown to describe the flow on energy along a raypath.
The differential equation for raypaths is the vector equation that is implied by the eikonal equa-
tion (2.74). For isotropic media, raypaths are normal to the wavefronts and the latter are described
# must be a vector that points in the direction of
as surfaces where T (#x) = constant. Therefore, ∇T
#
the raypath and the eikonal equation shows that ∇T = v −1 . Consider a wavefront at time T1 (#x1 )
and another at a slightly later time T2 (#x2 ) where T2 − T1 is very small. Let P1 be a point on the
surface T1 and P2 be the corresponding point on T2 along the normal from P1 (Figure 2.28). Then
∆s = (T2 (P2 )−T1 (P1 ))/v(P1 ) is an estimate of the perpendicular distance between these wavefronts,
and
# (P1 ) = lim T2 (P2 ) − T1 (P1 ) ŝ = ŝ
∇T (2.75)
P2 →P1 ∆s v(P1 )
where ŝ is a unit vector that is normal to the surface T1 (P1 ) or, in other words, ŝ points along the
raypath. If d#x is the differential vector pointing from P1 to P2 then ŝ may be written
d#x
ŝ = (2.76)
ds
where ds = |d#x| and s is the arclength along the raypath. Combining equations (2.75) and (2.76)
gives the raypath equation
# (#x) = ŝ = 1 d#x .
∇T (2.77)
v(#x) v(#x) ds
The traveltime T can be eliminated from equation (2.77) by calculating its derivative with respect
to arclength, which is
#
d ∇T d 1 d#x
= . (2.78)
ds ds v ds
2.12. RAYTRACING FOR INHOMOGENEOUS MEDIA 61
#
d ∇T # dT # 1
=∇ =∇ . (2.79)
ds ds v(#x)
The last step in this equation is justified using equation (2.77) and the identity d/ds = ŝ · ∇. #
#
Intuitively, this step is valid because ∇T points along the raypath and therefore dT /ds, that is the
# . Thus
scalar derivative with respect to arclength along the raypath, gives the full magnitude of ∇T
the ray equation is recast as
1 # d 1 d#x
2
∇v(#x) = − (2.80)
v (#x) ds v(#x) ds
which is a second order ordinary differential equation for the raypath vector #x. This result can be
further recast into a system of first-order equations by modifying the right-hand-side using
d dt d 1 d
= = (2.81)
ds ds dt v(#x) dt
1 d#x 1 d#x
p# = = 2 . (2.82)
v(#x) ds v (#x) dt
These last two equations allow equation (2.80) to be recast as the first-order system
d#x
= v 2 (#x)#
p (2.83)
dt
and
d#
p # x)
∇v(#
=− . (2.84)
dt v(#x)
Verification that these two equations (2.83) and (2.84) are equivalent to equation (2.80) can be done
by solving equation (2.83) for p# and substituting it into equation (2.84).
Example 2.12.1. As an illustration of the validity of equations (2.83) and (2.84), it is instructive to
show that they reduce to the v(z) case considered previously. If v(#x) = v(z), then in two dimensions,
let #x = [x, z] and p# = [px , pz ] so that equation (2.84) becomes
1 dx sin θ
px = = = constant (2.87)
v(z) ds v(z)
62 CHAPTER 2. VELOCITY
500
1000
1500
meters
2000
2500
3000
3500
4000
0 2000 4000 6000 8000 10000 12000
meters
Figure 2.29: A ray is shown traced through v(x, z) = 1800 + .6z + .4x using ode45 as
shown in Code Snippet 2.12.1.
where sin θ = dx/ds has been used. Of course, this is Snell’s law for the v(z) medium. Of the
remaining two equations in the system (2.85) and (2.86), the first is redundant because p2x + p2z =
v −2 (z) and the second gives
dz dz dz
dt = = = . (2.88)
v 2 (z)pz v(z) cos θ v(z) 1 − v 2 (z)p2x
When integrated, this will result in equation (2.57). In that equation, p is the same as px here.
d#r
= #a. (2.90)
dt
Then, given this prescription for d#r/dt and initial values for #r, the fourth-order Runge-Kutta (RK4)
method is invoked to integrate equation (2.90) for #r(t).
The MATLAB implementation requires a velocity matrix giving v(x, z) on a grid with ∆x =
2.12. RAYTRACING FOR INHOMOGENEOUS MEDIA 63
∆z ≡ ∆g and uses the convention that (x = 0, z = 0) is in the upper left corner. Prior to
raytracing, the function rayvelmod is invoked with arguments being the velocity matrix and the
grid spacing ∆g. This function creates a number of global variables that will be used repeatedly
in the raytracing. These include matrices of v 2 , ∂x ln v, and ∂z ln v that are needed in equation
(2.89). This precomputation speeds the raytracing and is especially beneficial if a great many rays
are to be traced; however, the cost is three times the memory overhead of the simple velocity matrix.
rayvelmod need only be called once at the beginning of the raytracing unless the velocity model is
changed.
The function drayvec implements the computation of d#r/dt according to equation (2.89). This
function is designed with the interface required by MATLAB’s built-in ODE solver ode45. This
latter function implements an RK4 scheme by calling a user-designed function like drayvec that
computes the time derivative of the solution vector. The general interface required by ode45 for
such a function is that it must have two input arguments that are the current time and the current
value of the solution vector #r. Thus, even though equation (2.90) does not require the value of t to
compute d#r/dt, drayvec requires t as input but does not use it.
Code Snippet 2.12.1. This code builds a velocity matrix representing v(x, z) = 1800 + .6z + .4x
and then uses ode45 to trace a ray. It creates Figure 2.29.
1 dg=10; %grid spacing
2 tstep=0:.004:3; %time step vector
3 x=0:dg:10000;z=x’; %x and z coordinates
4 v=1800+.6*(z*ones(1,length(x)))+.4*(ones(length(z),1)*x);%velocity
5 rayvelmod(v,dg);clear v;%initialize the velocity model
6 theta=pi*45/180;%takeoff angle
7 r0=[0,0,sin(theta)/1800,cos(theta)/1800]’;%initial value of r
8 [t,r]=ode45(’drayvec’,tstep,r0);%solve for raypath
9 plot(r(:,1),r(:,2));flipy%plot
End Code
Code Snippet 2.12.1 illustrates the tracing of a single ray in a medium with both vertical and
horizontal velocity gradients. The time-step vector is established on line 2 and a .004 second timestep
is used. Smaller timesteps will improve accuracy but will also lengthen the computation. For a
particular velocity model, some testing may be required to determine the optimal time step for the
problem at hand. The velocity matrix is built on line 4 and then passed to rayvelmod on line 5. It
is then cleared to save space because rayvelmod has established the required velocity information
in global matrices. The takeoff angle and initial values for #r are calculated on lines 7 and 8. As
p = v −2
mentioned in exercise 2.12.1, the components of p# are not independent of one another since p# ·#
and this is illustrated in the calculation of the initial values. Finally, ode45 is invoked to integrate
equation (2.90) and thus compute #r for the vector of times established on line 2. The calculated ray
is shown in Figure 2.29.
In the computation of d#r/dt it is generally true that the ray coordinates at a particular time will
not coincide with grid points in the model. Thus the estimation of the velocity terms in equations
(2.89) requires some form of interpolation. By default, drayvec uses nearest neighbor interpolation
because this is the fastest method. For more accuracy, bilinear interpolation is also available and its
use is controlled through a global variable that is explained in the drayvec help file.
A problem with the use of ode45, that is apparent upon close inspection of Figure 2.29, is that
the ray has been traced beyond the bounds of the velocity model. To avoid indexing into the velocity
array beyond its bounds, drayvec has been constructed to use the first (or last) column to represent
a simple v(z) medium for rays than go beyond the beginning (or end) of the model. A similar strategy
is used to cope with rays that go beyond the minimum or maximum depths. Though this allows
64 CHAPTER 2. VELOCITY
the solution to proceed, a better approach would be to detect when a ray has reached the boundary
of the model and stop the raytracing. For this purpose, an RK4 solver has been built as described
in Press et al. (1992) and is incorporated in the functions shootrayvxz and shootrayvxz g . In
these functions, the ray is tested as each time step to determine if it is within the model bounds
and raytracing is halted before the maximum time if the ray leaves the model. The syntax to
trace a ray with either program is essentially the same as with Code Snippet 2.12.1 except that the
command [t,r]=shootrayvxz(tstep,r0) replaces [t,r]=ode45(’drayvec’,tstep,r0) on line 8.
The functions shootrayvxz and shootrayvxz g differ in that the latter calls drayvec to compute
d#r/dt while the former does this computation directly without the function call. The result is that
shootrayvxz is much more efficient but does not offer bilinear interpolation as does shootrayvxz g .
If nearest-neighbor interpolation is satisfactory, then shootrayvxz should be preferred because it is
much faster.
Another useful raytracer is shootraytosurf . This function works similarly to shootrayvxz but
the ray is traced until it reaches z = 0 or a maximum time limit is exceeded. This is useful in normal
incidence raytrace modelling that will be discussed in chapter 4 in connection with normal incidence
raytrace migration.
Chapter 3
Wave propagation
3.1 Introduction
Seismic wave propagation in the upper layers of the earth’s crust is a complex physical process.
For example, a wave emanating from an explosive charge in a shallow (5 to 20 m) hole generally
undergoes an immediate and rapid spectral decay before it leaves the near surface1 . Theory predicts
that this spectral decay is associated with minimum phase effects that produce the characteristic
source waveform (Futterman, 1962). In addition to the primary downgoing pulse, there is always
upgoing energy from the source that reflects off of the earth’s free surface to produce a so-called
source ghost. In fact, there will be an entire series of similar ghosts that are produced when either
the primary or a particular ghost reflects off the base of the near surface and travels up again to
reflect off the free surface. Thus, it must always be expected that the downgoing pulse that serves
to illuminate the subsurface consists of a complex reverberatory series.
The subsurface is commonly idealized as a stacked sequence of nearly homogeneous layers, called
formations, whose geometry may be simple or complex. The contacts between the formations are
called horizons and are assumed to be in welded contact. (Of course, slip does occur along faults in
the subsurface but this is negligible over the time-span of a seismic experiment.) In the case of a
sedimentary basin, the formations are assumed to be horizontal and hence their material description
depends only on the vertical coordinate. This book adopts the common convention of a right-
handed Cartesian coordinate system with the z coordinate oriented in the vertical and increasing
downward. Thus it is common to describe the sedimentary basin environment as a v(z) setting by
which it is meant that the propagation speed of seismic waves depends only upon the formation
depth. (However, even in this setting, it is rarely a good model to assume that the base of the near
surface is a horizontal boundary.)
More complex settings are found near mountain belts, along passive continental margins, and
near the boundaries of sedimentary basins. Here, formations can be greatly distorted from the
horizontal due to tectonic compressional or extensional forces. Common jargon for such settings is
to say that they are v(x, z) environments which means that seismic wave speed can depend arbitrarily
upon position.
As a wave propagates into the subsurface, a number of important physical effects occur. When
1 The phrase near surface refers to a deliberately vague region just below the earth’s surface. It is generally
considered to be a region of high attenuation where static delays occur. Near surface effects are expected to be
surface consistent.
65
66 CHAPTER 3. WAVE PROPAGATION
a wave encounters a horizon (with different material parameters on either side) both reflected and
transmitted waves result. Regardless of the type of incident wave, reflections and transmissions are
generally found of each possible type of waves called compressional or P , shear or S, and perhaps
interface waves. Even if a source can be devised to produce only one type of wave, the propagating
wavefield rapidly evolves to contain all possible types going in all possible directions. For idealized
elastic media and plane waves 2 the equations governing reflection and transmission at a plane
interface are known exactly. These algebraically complex equations are given by Aki and Richards
(1980) and many approximate forms are also known.
Waves also lose energy as they propagate. Even in a perfectly elastic medium, a spherical
wavefront suffers a 1/r amplitude loss that is a geometric consequence of spherical divergence (or
spreading). In real media, there is always some additional loss due to anelasticity and typically
quantified by the dimensionless quality factor Q. These anelastic effects are often modelled by the
constant Q theory (Kjartansson, 1980) which refers to a Q that is independent of frequency but may
vary with position. All anelastic loss theories predict that attenuation is necessarily accompanied
by dispersion or they violate causality.
These introductory remarks only hint at the complexity of the true seismic wave propagation
process in the earth’s crust. They are mentioned here to give the reader some idea of the drastic
simplifications that will be seen in the theoretical discussions in this chapter. Producing realistic
numerical models of seismic wave propagation is a blend of both science and intuition. The science
must be pushed to limits dictated by the available computation power but then approximations must
be made. The science can be dealt with logically; but, knowing which physical effects are small and
can be neglected or developing an approximate expression for a mathematical operator, is often a
very subjective process.
This chapter will examine computational methods of wave propagation based mostly on the
scalar wave equation. Some initial effort will be expended in justifying the mathematical form of the
equation and understanding its range of validity. Then raytracing methods for simple media will be
developed and used to illustrate fundamental effects. Next methods for the formal solution of the
wave equation will be developed. The Fourier method will be shown to emerge from the classical
technique of separation of variables. Then Green’s theorem will be used to develop solutions which
proceed by integrating over the boundary (or initial) values. Finally, finite difference methods will
be explored to develop very general solution methods. (Note to the reader: this paragraph refers to
the chapter as I hope it to be. It is far from that now. Sorry!)
2 Plane waves are a mathematical idealization because they are necessarily infinite in extent. A point source actually
emits spherical wavefronts which may be approximately planar if the radius of the wavefront is large. A spherical
wavefront may be mathematically decomposed into plane waves for analysis.
3.2. THE WAVE EQUATION DERIVED FROM PHYSICS 67
6
6
ψ(x)
F
equilibriumWposition
x x+dx
Figure 3.1: Tension, T acts on a vibrating string whose displacement is ψ(x) to produce a restoring
force F
known to fall into three general classes: hyperbolic, parabolic, and elliptic. A discussion of the origin
of this classification and the meaning of these terms is beyond the scope of this book and the reader
is invited to consult a more general reference. There are many excellent works and some that are
particular favorites of this author are Morse and Feshbach (1953), Zauderer (1989), Ames (1992),
and Durran (1999). These books are very practical in their orientation and cannot be said to present
the subject of pde’s with full mathematical rigor. For this purpose, the sophisticated reader may
wish to consult Gustafson (1987) or, for the very brave, Taylor (1996).
In this section, some physical systems that lead to wave equations are examined. Only systems
that are relevant to the seismic exploration problem will be discussed.
where F (x) is the force per unit length and the minus sign in the square brackets is needed because
the tension acts at x and x + dx in opposite directions.
Now, recall that the definition of the second derivative of ψ(x) is
∂ 2 ψ(x) 1 ∂ψ(x + dx) ∂ψ(x)
= lim − . (3.3)
∂x2 dx→0 dx ∂x ∂x
∂ 2 ψ(x)
F (x) = T . (3.4)
∂x2
This important result says that the restoring force (per unit length) depends on the local curvature
of the string. When the second derivative is positive, the curvature is concave (<) and when it
is negative the curvature is convex (=). Since a positive force is directed upward in Figure 3.1
and a negative force is downward, it is easy to see that this force does act to restore the string to
equilibrium.
If the string’s weight is negligible, then Newton’s second law (force = mass × acceleration) gives
∂ 2 ψ(x, t) ∂ 2 ψ(x, t)
T 2
=µ (3.5)
∂x ∂t2
where µ is the mass per unit length and the dependence of ψ on both space and time is explicitly
acknowledged. It is customary to rearrange equation (3.5) to give the standard form for the one-
dimensional scalar wave equation as
∂ 2 ψ(x, t) 1 ∂ 2 ψ(x, t)
2
− 2 =0 (3.6)
∂x v ∂t2
where v = Tµ is the wave velocity 3 This analysis has shown that the scalar wave equation arises
for the vibrating string as a direct consequence of Newton’s second law. Wave equations invariably
arise in this context, that is when a displacement is initiated in a continuum. Also typical is the
assumption of small displacements. Generally, large displacements lead to nonlinear equations.
The waves modelled here are known as transverse waves because the particle displacement is in a
directional orthogonal to the direction of wave propagation. In section 3.2.2 longitudinal waves, that
have particle oscillation in the direction of wave propagation, are discussed.
If there are external forces being applied to the string as specified by the force density 4 function
3 It is quite common to use velocity for this scalar quantity even though velocity is classically a vector quantity in
∂φ(x, t) ∂φ(x, t)
−a =0 (3.8)
∂x ∂t
where a is a constant.
Exercise 3.2.1. Show that the left side of equation (3.6) can be factored into two operators of the
form of the left side of equation (3.8). What values does the constant a take in each expression.
Show that ψ(x, t) = ψ1 (x + vt) + ψ2 (x − vt) is a solution to (3.6) by showing that the two factors
annihilate ψ1 or ψ2 . (ψ1 (x+vt) means an arbitrary one-dimensional function of x+vt and similarly
for ψ2 (x − vt)). Can you explain the physical meaning of this result?
The quantity θ is the volume strain or dilatation and is a measure of local fractional volume change.
This relation says that the pressure fluctuations are directly proportional to induced volume changes.
# · #ν . The minus sign in equation (3.9) is needed
θ is related to particle velocity, #ν , though ∂t θ = ∇
because an increase in pressure causes a decrease in volume and vice-versa. This is an example of a
constitutive relation that defines how stress relates to strain.
We assume that the fluid is at rest except for those motions induced by the pressure disturbance
P (#x, t). As with the vibrating string, the motion of the fluid is defined by Newton’s second law
4 In this context this is just a fancy way of saying force per unit length.
5 It is unfortunate that inhomogeneous is used in at least two very different contexts in this theory. Here the word
refers to physical parameters that vary with position. In the previous section the word referred to a pde with a source
term and that usage has nothing to do with this.
70 CHAPTER 3. WAVE PROPAGATION
pressure gradient point in the direction opposite to the gradient. That is, the gradient points from
low to high pressure but the fluid will move from high to low pressure.
The interpretation of the ddtν requires some care. Generally in fluid dynamics, such time deriva-
tives consist of two terms as dT ∂T #
dt = ∂t + (# ν · ∇)T where T denotes any scalar fluid property such as
temperature. The first term describes local temporal changes (e.g. local heating) while the second
describes changes due to new material being transported to the analysis location by fluid motion.
For this reason, the second term is known as the convective term. The assumption that the fluid is
at rest except for the motions caused by the passing pressure wave allows the convective term to be
neglected. If the term is included, a nonlinear wave equation will result. Thus, equation (3.10) is
rewritten approximately as
−∇P# (#x, t) = ρ(#x) ∂#ν (#x, t) . (3.11)
∂t
A wave equation for pressure can be derived by taking the divergence of (3.11),
#
# · ∇P
−∇ # x) · ∂#ν (#x, t) + ρ(#x) ∂ ∇ · #ν (#x, t) ,
# (#x, t) = ∇ρ(# (3.12)
∂t ∂t
and substituting equation (3.11) into the first term on the right side and ∂t of equation (3.9) into
the second. This gives
# ln ρ(#x) · ∇P
# (#x, t) − ρ(#x) ∂ 2 P (#x, t)
∇2 P (#x, t) − ∇ =0 (3.13)
K(#x) ∂t2
where ∇2 = ∇ # ·∇
# and 1 ∇ρ
# = ∇ln# ρ have been used. If the second term were absent in equation
ρ
(3.13) then we would have a classic scalar wave equation for pressure. This occurs exactly if ρ(#x) is
constant. In many other cases when ρ(#x) varies slowly the term will be negligible. Neglecting the
# ρ term then gives
∇ln
1 ∂ 2 P (#x, t)
∇2 P (#x, t) − =0 (3.14)
v2 ∂t2
where v = K ρ . Equation (3.14) is a scalar wave equation that approximately models the propa-
gation of pressure waves through an inhomogeneous fluid. Though a number of assumptions have
been made in deriving this result, it still has quite general application. Also, it is possible to use
a solution of (3.14) to develop an approximate solution to the more complex Equation (3.13). Let
# ln ρ · ∇P
P̂ (#x, t) be a solution to equation (3.14), then it can be used to approximately express the ∇ #
term in (3.13) to give
1 ∂ 2 P (#x, t) # ln ρ(#x) · ∇
# P̂ (#x, t).
∇2 P (#x, t) − =∇ (3.15)
v 2 (#x) ∂t2
The term on the right does not depend upon the unknown P (#x, t) but rather plays the role of a
source term in scalar wave equation. Thus the effect of strong density inhomogeneity introduces an
approximate source term whose strength is proportional to the logarithmic gradient of density. This
process can be repeated indefinitely with the solution from iteration n-1 being used to compute the
effective source term for iteration n.
Exercise 3.2.2. Show that The volume dilatation, θ also satisfies a scalar wave equation.
Exercise 3.2.3. Consider the 1-D inhomogeneous fluid. Write down the wave equation that ap-
3.3. FINITE DIFFERENCE MODELLING WITH THE ACOUSTIC WAVE EQUATION 71
This section describes a simple facility for finite difference modelling using the acoustic wave equa-
tion. This toolkit was originally developed by Carrie Youzwishen who was employed by this author
for that purpose. Subsequently it has been evolved by this author and Dr. Zhengsheng Yao. The
package provides a basic facility, using second order finite difference approximations, for modelling
with the variable velocity acoustic wave equation in two dimensions. The facility provides tools
for model building, source specification, time-stepping a wavefield, seismogram creation, and mak-
ing wavefield movies. Sources and receivers can be positioned anywhere within the 2D medium
(including on the boundaries) and absorbing boundary conditions are implemented.
Consider the variable-velocity scalar wave equation in two dimensions
where ∆t is the time discretization interval. This expression can be solved for the wavefield at t + ∆t
to give
ψ(x, z, t + ∆t) = 2 + ∆t2 v 2 (x, z)∇2 ψ(x, z, t) − ψ(x, z, t − ∆t). (3.18)
This is an expression for time stepping the wavefield. It shows that estimation of the wavefield at
t + ∆t requires knowledge of the two earlier wavefields at t and t − ∆t. Each of these wavefields is
called a snapshot and, in a computer simulation, they are all two dimensional matrices.
Equation (3.18) shows that ψ(x, z, t + ∆t) is estimated from the superposition of three terms:
2ψ(x, z, t), ∆t2 v 2 (x, z)∇2 ψ(x, z, t), and −ψ(x, z, t − ∆t). Of these, the second requires the compu-
tation of the Laplacian (∇2 ) on the current wavefield which is an expensive operation. MATLAB
supplies the function del2 that computes .25∇2 of a matrix using centered second-order finite oper-
ators that are modified near the boundary. Experimentation showed that del2 was not optimal for
use with absorbing boundary conditions so two alternate Laplacians, del2 5pt and del2 9pt , were
created. These functions both pad the entire boundary of the input matrix with extra rows and
columns of zeros and this works better with absorbing boundaries. del2 5pt implements ∇2 using
second-order finite-difference operators while del2 9pt uses fourth-order. Thus del 5pt computes
the approximation
∇2 ψ(x, z, t) =
1
[−ψ(x + 2∆x, z, t) + 16ψ(x + ∆x, z, t) − 30ψ(x, z, t) + 16ψ(x − ∆x, z, t) − ψ(x − 2∆x, z, t)] +
12∆x2
1
[−ψ(x + 2∆z, z, t) + 16ψ(x + ∆z, z, t) − 30ψ(x, z, t) + 16ψ(x − ∆z, z, t) − ψ(x − 2∆z, z, t)]
12∆z 2
(3.20)
The core function in the finite difference toolbox is afd snap (the ‘afd ’ prefix stands for acous-
tic finite difference). Function afd snap requires two input wavefields representing ψ(x, z, t) and
ψ(x, z, t − ∆t) and computes ψ(x, z, t + ∆t) according to equation (3.18). This function is intended
to be used in a computation loop that time-increments a wavefield any number of steps. Two initial
wavefields must be constructed to start the simulation and sources may be prescribed by placing
appropriate impulses in these two wavefields. Receivers are simulated by extracting samples at each
receiver location from each ψ(x, z, t) as it is computed. These extracted samples may be accumulated
in vectors representing the recorded traces.
The use of equation (3.18) in a time-stepping simulation is known to be unstable in certain
circumstances (Mitchell and Griffiths, 1980). Instability means that the amplitudes of the wavefield
grow without bound as it is stepped through time. The key to this behavior is the amplitude of
the ∇2 ψ(x, z, t) term in equation (3.18). Using the five-point Laplacian of equation (3.19) (with
∆z = ∆x) in equation (3.18) leads to
∆t2 v 2
ψ(x, z, t + ∆t) = 2ψ(x, z, t) − ψ(x, z, t − ∆t) + [δxx ψ(x, z, t) + δzz ψ(x, z, t)] (3.21)
∆x2
where δxx ψ(x, z, t) = ψ(x + ∆x, z, t) − 2ψ(x, z, t) + ψ(x − ∆x, z, t) and similarly for δzz . In this
expression, all of the ψ terms can be considered to have similar magnitude. Thus, the factor
∆t2 v 2 ∆x−2 is a possible amplification factor if it becomes too large. Lines et al. (1999) show that
the condition for stability is
v∆t 2
≤√ (3.22)
∆x a
where the constant a is the sum of the absolute values of the weights for the various wavefield terms
in the finite-difference approximation for ∇2 . For the Laplacian of equation (3.19), a = 8 while
for equation (3.20), a = 128/12 = 32/3. Also, since v is a function of (x, z) it suffices to use the
maximum velocity in the model. Thus the stability conditions are
vmax ∆t √1
2
second-order Laplacian,
≤ (3.23)
∆x 3 fourth-order Laplacian.
8
Thus, the time and space sample rates cannot be chosen independently. Generally, finite-
difference operators need many more samples than the Nyquist criterion of two per wavelength.
Technically, this is because the operators cause an artificial dispersion called grid dispersion. Grid
dispersion preferentially affects the shorter wavelengths so oversampling reduces the dispersion. A
good rule of thumb is about five to ten samples per wavelength. Typically, in the creation of a
model, a desired temporal frequency range is known. Then, the minimum wavelength is given
by λmin = vmin /fmax and the spatial sample rate can be chosen to achieve a desired number of
3.3. FINITE DIFFERENCE MODELLING WITH THE ACOUSTIC WAVE EQUATION 73
Code Snippet 3.3.1 illustrates the use of these finite difference facilities to model a single shot
record. First, a three-layer velocity model is defined on lines 1 through 11 by calling afd vmodel
twice. The velocity model is a matrix representing vins in (x, z) space. The function afd vmodel is a
convenience function that fills in a polygonal area in a matrix with a constant value. The polygonal
area is defined by the (x, z) coordinates of its nodes and the coordinate origin is assumed to be the
upper left corner of the model. In all of these finite difference codes, the vertical and horizontal
sample sizes must be identical. In this example, the velocity model has three layers with horizontal
interfaces and layer velocities of 2000 m/s, 2800 m/s, and 3200 m/s.
The model is built by first defining (x, z) coordinates (lines 2-3) and filling the velocity matrix
with 3200 m/s (line 5). Subsequently, two polygons are defined to represent the two upper layers.
On line 6, z1 and z2 are the depths to the bottom of the first and second layers. Line 8 defines the
(x, z) coordinates of the four vertices of a rectangle representing the first layer. A property of the
74 CHAPTER 3. WAVE PROPAGATION
meters
0 200 400 600 800 1000 1200
0
200
400
meters
600
800
1000
1200
algorithm used by afd vmodel is that points that lie exactly on the boundary of the polygon are
considered outside the polygon and so do not acquire the new velocity. The variable dx2, defined on
line 7 as half of the grid spacing, is used on line 8 to define a rectangle that is half a grid spacing above
depths 0 and z2 and half a grid spacing outside the minimum and maximum x coordinates. This
ensures that the rectangle extends to the limits of the velocity matrix. Line 9 fills this rectangle with
the velocity of 2000 m/s and then lines 10 and 11 repeat this process for the next layer. The resulting
velocity model is shown in Figure 3.2. This plot was made using plotimage(vmodel-3000,z,x).
A constant is subtracted from the velocity model so that the resulting matrix has both positive
and negative values as expected by plotimage . The raypaths shown in this figure correspond to
traveltimes shown in Figures 3.3 and 3.4.
Code Snippet 3.3.1 creates two seismograms, the first (line 21) uses a second-order Laplacian
(equation 3.19) and the second (line 23) uses a fourth-order Laplacian (equation 3.20). The prepa-
ration for the seismograms defines the time step (line 13), temporal sample rate (line 14), maximum
time (line 14), the receiver locations (lines 15-16) and the source strength and geometry (line 19).
The time step is generally chosen to be small enough to satisfy the stability requirement (equation
3.23) and the temporal sample rate is usually much more coarse. afd shotrec internally calculates
the seismogram at the sample rate of the time step and then resamples it to the desired temporal
sample rate. This is sensible because it is a well-known property of finite-difference methods that
the higher frequencies are physically inaccurate. In this case, the Nyquist frequency for a sample
rate of .001 seconds is 1000 Hz while for .004 seconds it is 125 Hz. Considering that the antialias
filter for resampling to .004 seconds will attenuate frequencies above about half of Nyquist, this ex-
ample anticipates using only frequencies below about 60 Hz and that is less than 10% of the original
Nyquist.
The variables snap1 and snap2 created on lines 17-19 represent the the wavefield ψ(x, z, t = −∆t)
and ψ(x, z, t = 0) respectively. They are created to be the same size as the velocity model and snap1
is a matrix of zeros. The source is specified by placing appropriate nonzero samples in snap2. In
this case, a point source in the center of the x axis at z = 0 is simulated.
At this stage, all of the input parameters to afd shotrec have been described except for the final
three. These specify the filter (or wavelet) to be convolved with the seismogram, the phase of the
filter, and the order of the Laplacian. In this case, Ormsby filter specifications are given as [5 10 30
40] and this means that the filter pass band begins at 5 Hz., reaches full ampltude at 10 Hz., begins to
3.3. FINITE DIFFERENCE MODELLING WITH THE ACOUSTIC WAVE EQUATION 75
0 0
0.2 0.2
0.4 0.4
seconds
seconds
0.6 0.6
0.8 0.8
1 1
0 200 400 600 800 1000 1200 0 200 400 600 800 1000 1200
meters meters
Figure 3.3: The second-order seismogram cre- Figure 3.4: The fourth-order seismogram cre-
ated on line 21 of Code Snippet 3.3.1. Superim- ated on line 23 of Code Snippet 3.3.1. See Figure
posed on the right-hand side are raytraced trav- 3.3 for a description of the superimposed travel-
eltimes for (from top to bottom) the first pri- times.
mary reflection, the first-order multiple in the
top layer, and the second primary reflection.
The corresponding raypaths are shown in Fig-
ure 3.2.
ramp down at 30 Hz., and rejects frequencies above 40 Hz. The penultimate parameter specifies the
phase of the filter that is, in this case, zero (a value of 1 gives minimum phase). The last parameter
is either 1, indicating a second-order Laplacian, or 2, meaning a fourth-order approximation.
The seismograms are shown in Figures 3.3 and 3.4 respectively. Both of these figures are plotted
with a highly clipped display, otherwise, the direct wave near the source would be the only visible
arrival. Though similar at first glance, these two seismograms are significantly different. The three
hyperbolic reflections on each figure are (from the top) the primary reflection off the bottom of the
first layer, the first-order multiple reflecting between the bottom of the first layer and the surface,
and the primary reflection off the bottom of the second layer. The traveltimes for these events
are superimposed on the right-hand sides of the figures and their raypaths are shown in Figure
3.2. The reflections for these events lag significantly behind the raytraced traveltimes but more
so on the second-order solution. This time lag occurs because the finite difference approximations
to the derivative in the wave equation have a frequency-dependent performance. Essentially, for
wavelengths that are large compared to the computation grid, the approximate derivatives are
acceptable but, as the wavelength approaches the grid spacing, the approximation becomes very poor.
The result is a phenomenon known as grid dispersion, meaning that the actual propagation speed of
waves on the grid is not simply vins (x, z) but is a complex function of wavelength (or, equivalently,
frequency). Grid dispersion makes the reflections appear to lag behind the corresponding raytrace
traveltimes. It also makes the apparent wavelet appear more elongated (dispersed). Comparison of
the first reflection at far offsets on both seismograms shows that the second-order result has a more
dispersed waveform. Also, the fourth-order result has a smaller lag behind the raytrace traveltimes
than the second-order for all reflections.
Careful inspection of both Figures 3.3 and 3.4 shows a set of events that originate where each
76 CHAPTER 3. WAVE PROPAGATION
reflection meets the sides of the seismogram and then cross in the center of the figures. These are
typical artifacts known as edge effects. They arise because a computer simulation of wave propagation
must always operate in a finite domain with definite boundaries. When a wave encounters such a
boundary, it behaves as though it has met a perfect reflector. These boundary reflections would be
much worse if afd shotrec did not incorporate absorbing boundaries (Clayton and Engquist, 1977).
Obviously, absorbing boundaries are not perfect but they do help. The boundary related artifacts
can be seen to be generally quite steeply dipping in (x, t). This is because the absorbing boundary
conditions are optimized for a small range of wavefront incidence angles around normal incidence.
That is, a wavefront that encounters the boundary at normal incidence is completely absorbed while
one making an angle of, say, 30◦ is partially reflected.
Chapter 4
77
78 CHAPTER 4. ELEMENTARY MIGRATION METHODS
migration algorithms that are generally distinguished by the kind of approximations made. As a
result there are many different, overlapping ways to categorize migration algorithms. For example,
it can be shown that an assumption that waves are traveling only upward at the earth’s surface,
though physically incorrect, allows the solution to proceed without knowledge of ∂z ψ. (Or equiv-
alently, this assumption allows ∂z ψ to be calculated from ψ.) This approach is called the one-way
wave assumption and is very common. Given that one-way waves will be considered, there are: finite
difference algorithms that offer great flexibility with spatial variations of velocity but suffer from grid
dispersion and angle limitations, Kirchhoff methods that can easily accommodate arbitrary varia-
tions in seismic recording geometry but require raytracing to guide them, and Fourier (or spectral)
techniques that offer high fidelity but have difficulty coping with rapid variations in either velocity
or recording geometry. All of these techniques attempt, in some manner, to downward continue
measurements made at z = 0 into the subsurface. There are also useful distinctions about whether
the estimation of R(x, y, z) is made directly from z = 0, direct methods, or whether R(x, y, z) is
estimated from z − ∆z, recursive methods. Finally, perhaps the most common categorization comes
under the confusing labels of time migration and depth migration. The former is an earlier tech-
nology that remains viable because it is very insensitive to velocity errors even though it is known
to be accurate only if ∂x v ∼ 0. On the other hand, depth migration is the only currently available
imaging technology that is accurate when ∂x v varies strongly; however, it requires a very accurate
specification of the velocity model.
These last remarks hint at the general chicken-and-egg nature of modern migration methods.
A true inversion technique would derive the velocity model from the data as part of the inversion.
Migration requires the velocity model as input and merely attempts to reposition, focus, and adjust
amplitudes. These really all turn out to be the same thing. In order to be a useful process, it should
be true that migration requires only an approximate, or background, velocity model and that more
detailed velocity information can be extracted from a successful migration than was input to it.
This is generally the case though it can be very difficult to achieve in structurally complex areas.
Especially for depth migration, the construction of the velocity model is the central difficulty in
achieving a successful result.
This chapter will discuss “elementary” migration methods. This refers to techniques designed
to migrate stacked seismic data (post-stack migration) in simple velocity variations. More complex
topics including pre-stack migration will be discussed in the next chapter.
At best, equation (4.1) can only be an expression for normal incidence reflectivity. More generally,
reflectivity must be acknowledged to be a function of the angle of incidence as well as spatial position.
4.1. STACKED DATA 79
For the case of plane elastic waves in layered media, the Zoeppritz equations (Aki and Richards,
1980) are an exact prescription of this variation. The Zoeppritz equations are highly complex and
nonlinear and the estimation of their approximate parameters from seismic data is called the study
of amplitude variation with offset or AVO. (A nearly synonymous term is amplitude variation with
angle or AVA.) Such complex reflectivity estimates are key to true lithology estimation and are a
subject of much current research. Ideally, AVO analysis should be conducted simultaneously with
migration or after migration. For the current discussion, it is assumed that only the normal-incidence
reflectivity, r(x, y, z), is of interest and that a stacked seismic section provides a bandlimited estimate
of r(x, y, z).
The estimate of reflectivity must always be bandlimited to some signal band of temporal fre-
quencies. This is the frequency range over which signal dominates over noise. Unless some sort of
nonlinear or model-based inversion is done, this signal band is determined by the power spectrum
of the source, the data fold, the types of coherent and random noise, and the degree of anelastic
attenuation (Q loss). The optimal stacked section will have a zero-phase, white, embedded wavelet.
This means that it should obey the simple convolutional model
where w(t) is a zero phase wavelet whose amplitude spectrum is white over some limited passband.
If w(t) has residual phase or spectral color then these will adversely affect the resolution of the final
migrated image. They should be dealt with before migration.
Spherical spreading correction: A time dependent amplitude scaling that approximately corrects
for spherical divergence.
Deconvolution: The source signature and the average Q effect is estimated and removed.
Sort to (x, h) coordinates: The data is considered to be gathered in (s, g) (source,geophone) coor-
dinates and then sorted to midpoint, x, and half-offset, h (Figure 4.1).
Statics correction: Near surface static time delays are estimated and removed.
NMO removal: The data is taken through stacking velocity analysis and the resulting velocities are
used to remove normal moveout time delays.
Residual statics correction: Residual time delays (small) are sought and removed. Usually this is
a surface-consistent step.
Trace balancing: This might be considered optional but some step is needed to correct for source
strength and geophone coupling variations.
80 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Common
offset
24
22
s/r s/r s/r s/r s/r s/r s/r s/r
Source coordinate
20
18 Common
16
midpoint
14
12
VelocityWlens
10
8
6
4
0 5 10 15 20 25 30
Receiver coordinate
Figure 4.1: A 2D seismic line is depicted in the Figure 4.2: Most, but not all, zero-offset ray-
(s, g) plane. Each shot is represented as having paths are also normal-incidence raypaths.
six receivers with a central gap. Common mid-
point, x, and common offset, h coordinates are
depicted.
CMP stacking: All traces with a common midpoint coordinate are summed. Events which have
been flattened on CMP gathers are enhanced. Other events and random noise are reduced.
The ZOS is said to be signal enhanced because the stacking process has been designed to select
against multiples. The simplest (and dominant) raypath which returns energy to its source is called
the normal-incidence raypath. Quite obviously, the normal-incidence reflection will send energy
back up along the path that it traveled down on. However, it is possible to have zero-offset paths
that are not normal-incidence paths (Figure 4.2).
The relationship between (s, g) and (x, h) is depicted graphically in Figure 4.1.
Mathematically, a 2D pre-stack seismic dataset is a 3D function, ψ(s, g, t). This dataset can be
written as an inverse Fourier transform of its spectrum through
ψ(s, g, t) = φ(ks , kg , f )e2πi(ks s+kg g−f t) dks dkg df (4.5)
V∞
4.1. STACKED DATA 81
35
30
25
∆ g (meters)
20
15
10
2500 m/s
2000 m/s
5 1500 m/s
1000 m/s
500 m/s
0
50 100 150 200
maximum frequency (Hertz)
Figure 4.3: These curves show the maximum spatial sample rate required to avoid aliasing
of kg on shot records (assuming point receivers) as a function of the maximum signal
frequency. A curve is drawn for each of five near surface velocities. This is based on
expression (4.15).
where the notation V∞ indicates the integration covers the entire relevant, possibly infinite, portion
of (ks , kg , f ) space. Imposing the coordinate transform of equations 4.4 gives
ψ(x, h, t) = φ(ks , kg , f )e2πi(ks (x−h)+kg (x+h)−f t) dks dkg df (4.6)
V∞
or
ψ(x, h, t) = φ(ks , kg , f )e2πi((ks +kg )x+(ks −kg )h−f t) dks dkg df. (4.7)
V∞
kx = ks + kg and kh = ks − kg (4.8)
where the factor of 1/2 comes from the Jacobian of the coordinate transformation given in equation
4.9. Both of equations 4.5 and 4.10 are proper inverse Fourier transforms which shows that the
wavenumber relationships given in equations 4.8 and 4.9 are the correct ones corresponding to the
(s, g) to (x, h) coordinate transformation.
Quite a bit can be learned from equations 4.9 about how the spectral content of the pre-stack
data gets mapped into the post-stack data. The maximum midpoint wavenumber, kxmax comes
82 CHAPTER 4. ELEMENTARY MIGRATION METHODS
from the sum of ksmax and kgmax . As will be shown formally later, the maximum wavenumber
is directly proportional to horizontal resolution (i.e. the greater the wavenumber, the better the
resolution). The meaning of the maximum wavenumbers is that they are the largest wavenumbers
that contain signal. They typically are less than the Nyquist wavenumbers but they cannot be
greater. Thus ksmax ≤ ksN yq = 1/(2∆s) and kgmax ≤ kgN yq = 1/(2∆g). It is customary to sample
the x coordinate at ∆x = .5∆g so that kxN yq = 1/(2∆x) = 2kgN yq which means that the stacking
process will generate kx wavenumbers up to this value from combinations of ks and kg according
to equation 4.8. However, because spatial antialias filters (source and receiver arrays) are not very
effective, it must be expected that wavenumbers higher than Nyquist will be present in both the ks
and kg spectra. The only way to generate unaliased wavenumbers up to kxN yq in the stacking process
is if ∆s = ∆g so that kxN yq = ksN yq + kgN yq . This means that there must be a shotpoint for every
receiver station which is very expensive acquisition. Normal 2D land shooting puts a shotpoint for
every n receiver stations where n ≥ 3. This means that kx wavenumbers greater than a certain limit
will be formed from completely aliased ks and kg wavenumbers. This limiting unaliased wavenumber
is
1 1 1 1 n+1
kxlim = + = + = . (4.11)
2∆g 2∆s 2∆g 2n∆g 2n∆g
For n=3, kxlim = 23 kxN yq so that shooting every third group will result in a conventional stack with
the upper third of the kx spectrum being completely useless.
Even if n = 1, there will still be aliased contributions to kx if the ks and kg spectra are aliased
as they normally are. For example, kxN yq can be formed with unaliased data from ksN yq + kgN yq
but it can also be formed from the sum of ks = 0 and kg = 2kgN yq and there are many more such
aliased modes. Similarly, wavenumbers less that kxlim can be formed by a wide variety of aliased
combinations. Further analysis is helped by having a more physical interpretation for ks and kg so
that their potential spectral bandwidth can be estimated.
As explained in section 2.9, the apparent velocity of a wave as it moves across a receiver array is
given by the ratio of frequency to wavenumber. For a source record (common source gather ), this
has the easy interpretation that
f v0
= (4.12)
kg sin θ0
which means that
f sin θ0
kg = . (4.13)
v0
So, if fmax is the maximum temporal frequency and (sin θ0 )max = 1 then
fmax
kgmax = (4.14)
v0min
where v0min is the slowest near surface velocity found along the receiver spread. So, to avoid any
aliasing of kg , assuming point receivers, the receiver sampling must be such that kgN yq ≥ kgmax
which leads to the antialiasing condition
v0min
∆g ≤ . (4.15)
2fmax
Assuming the equality in expression (4.15) allows the maximal sampling curves in Figure 4.3 to
be calculated. Alternatively, if coarser sampling than suggested in these curves is planned, then
a portion of the kg spectrum will be aliased. An array can be designed to suppress the aliased
4.1. STACKED DATA 83
source
receivers
α θ
receiver
sources
α θ
B
Figure 4.4: (A) Wavenumber kg is measured on a common source gather and estimates
the emergence angle, θ0 , of a monochromatic wave. (B) Wavenumber ks is measured on a
common receiver gather and estimates the takeoff angle α0 of the ray that arrives at the
common receiver.
If equation (4.10) is substituted into equation (4.17) and the order of integration reversed, it results
that
1
ψ0 (x, t) = e2πikh h dh φ(kx , kh , f )e2πi(kx x−f t) dkx dkh df. (4.18)
2
V∞ all h
The integral in square brackets is δ(kh ) which, in turn, collapses the kh integral by evaluating it at
84 CHAPTER 4. ELEMENTARY MIGRATION METHODS
R
Surface depth=0
Q P
Reflector S
depth=z
w
Figure 4.5: The zero-offset Fresnel zone is defined as the width of a disk on a subsurface
reflector from which scattered energy arrives at R with no more than a λ/2 phase difference.
kh = 0 to give
1
ψ0 (x, t) = φ(kx , 0, f )e2πi(kx x−f t) dkx df. (4.19)
2
V∞
So the CMP stacking process passes only the kh = 0 wavenumber which is a very severe f -kfilter
applied to a CMP gather. It is for this reason that f -kfiltering applied to CMP gathers does
not typically improve the stack. However, f -k filtering of source and receiver gathers can have a
dramatic effect. The kh = 0 component is the average value over the CMP gather and corresponds
to horizontal events. Of course, this is done after normal moveout removal that is designed to flatten
those events deemed to be primaries.
So far, the discussion has been about wavenumber spectra with nothing being said about the
effect of stacking on temporal frequencies. Bearing in mind the goal to estimate bandlimited reflec-
tivity, the ideal stacked section should have all spectral color removed except that attributable to
the reflectivity. Most pre-stack processing flows rely on one or more applications of statistical decon-
volution to achieve this. Typically these deconvolution techniques are not capable of distinguishing
signal from noise and so whiten (and phase rotate) both. If the pre-stack data input into the stacking
process has been spectrally whitened, it will always suffer some attenuation of high frequencies due
to the
√ stacking out of noise. This is because the stacking process improves the signal-to-noise ratio
by f old. Moreover, this effect will be time-variant in a complicated way because of two competing
effects. First, the fold in any stacked section is actually time-variant because of the front-end mute.
Until the time at which nominal fold is reached, the stacking effect is therefore time variant. Second,
the signal-to-noise ratio in any single pre-stack trace must decrease with both increasing time and
increasing frequency due to attenuation (i.e. Q effects). As a result, it is virtually guaranteed that
the post-stack data will need further spectral whitening. It is also highly likely that some sort of
wavelet processing will be required to remove residual phase rotations.
1200
20Hz.
1000
Fresnel zone (meters)
800
Original Fresnel zone
600
60Hz.
w
100Hz. λ/2
400 140Hz.
180Hz.
Figure 4.6: The width of the Fresnel zone (based Figure 4.7: A 3D migration collapses the Fres-
on equation (4.21)) is shown versus depth for a nel zone to a disk of diameter λ/2 while a 2D
variety of frequencies and an assumed velocity migration only collapses the Fresnel zone in the
of 3000 m/s. inline direction
disk whose center is marked by the ray P and edge by the ray Q. Ray P is a simple zero-offset ray
while ray Q is defined to be λ/4 longer than ray P. Assuming constant velocity, this means that
scattered energy that travels down-and-back along path Q will be λ/2 out of phase at the receiver,
R, compared with that which travels along path P. All paths intermediate to P and Q will show a
lesser phase difference. Therefore, it is expected that path Q will interfere destructively with P and
is taken to define the diameter of the central disk from which scattered energy shows constructive
interference at R. The equation defining the Fresnel zone diameter is therefore
λ
z 2 + (w/2)2 − z = (4.20)
4
which is easily solved for w to give
√
w = 2 (λ/4 + z)2 − z 2 = 2zλ 1 + λ/(8z) = 2zv/f 1 + v/(8f z) (4.21)
where,
√ in the third expression, λf = v, has been used. If z λ, then the approximate form
w ∼ 2zλ is often used. The size of a Fresnel zone can be very significant in exploration as it is
often larger than the target (Figure 4.6).
Migration can be conceptualized as lowering the source and receiver to the reflector and so shrinks
the Fresnel zone to its theoretical minimum of λ/2. That the Fresnel zone does not shrink to zero is
another way of understanding why reflectivity estimates are always bandlimited. However, seismic
data is a bit more complicated than this in that it possesses signal over a bandwidth while the
Fresnel zone concept is monochromatic. Roughly speaking, the effective Fresnel zone for broadband
data can be taken as the average of the individual Fresnel zones for each frequency.
The 3D Fresnel disk collapses, under 3D migration, from a diameter of w to a diameter of λ/2.
However, much seismic data is collected along lines so that only 2D migration is possible. In this
case, The Fresnel zone only collapses in the inline direction and remains at its unmigrated size in
the crossline direction (Figure 4.7). Thus the reflectivity shown on a 2D migrated seismic line must
be viewed as a spatial average in the crossline direction over the width of the Fresnel zone. This is
one of the most compelling justifications for 3D imaging techniques, even in the case of sedimentary
86 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Figure 4.8: Two normal-incidence rays from a dipping reflector beneath a v(z) medium.
The delay between the rays allows the ray parameter to be estimated.
comes from a popular, though somewhat incorrect, view that the process just “moves events around”
on a seismic section. This encourages a fundamental, and very pervasive, confusion that time and
depth sections are somehow equivalent. A better view is that an unmigrated time section and a
migrated depth section are independent (in fact orthogonal) representations of a wavefield. Raytrace
migration emphasizes this latter view and makes the role of the velocity model very clear.
Figure 4.8 illustrates the link between reflector dip and the measurable traveltime gradient (for
normal-incidence rays) on a ZOS. The figure shows two rays that have normal incidence on a dipping
reflector beneath a v(z) medium. More general media are also easily handled, this just simplifies the
discussion. The two rays are identical except for an extra segment for the lower ray in the bottom
layer. From the geometry, this extra segment has a length ∆x sin δ where ∆x is the horizontal
distance between the rays. If there were receivers on the interface between layers three and four,
then for the upcoming rays, they would measure a traveltime delay, from left to right, of
2∆x sin δ
∆t = (4.25)
v4
where the factor of two accounts for two-way travel. This implies a horizontal traveltime gradient
(horizontal slowness) of
∆t sin δ
=2 . (4.26)
∆x v4
Inspection of the geometry of Figure 4.8 shows that equation (4.26) can be rewritten as
∆t
= 2pn (4.27)
∆x
where pn is the ray parameter of the normal-incidence ray. In the v(z) medium above the reflector,
Snell’s law ensures that pn will remain equal to twice the horizontal slowness at all interfaces including
the surface at z = 0. Thus, a measurement of ∆t/∆x at the surface of a normal-incidence seismogram
completely specifies the raypaths provided that the velocity model is also known. In fact, using
equation (4.26) allows an immediate calculation of the reflector’s dip.
Of course, in addition to the reflector dip, the spatial position’s of the reflection points of the
normal rays are also required. Also, it must be expected that more general velocity distributions
88 CHAPTER 4. ELEMENTARY MIGRATION METHODS
θ
Datum (recording surface)
φ
θ v
φ
θ v
φ θ v
φ v
θ
θ
δ v
Figure 4.9: The migration of a normal-incidence ray is shown. Knowledge of the ray
parameter allows the ray to be traced down through the velocity model.
than v(z) will be encountered. In v(x, y, z) media, the link between reflector dip and horizontal
slowness is not as direct but knowledge of one still determines the other. Normal-incidence raytrace
migration provides a general algorithm that handles these complications. This migration algorithm
will be described here for 2D but is easily generalized to 3D. Before presenting the algorithm, a
formal definition of a pick is helpful:
A pick is defined to be a triplet of values (x0 , t0 , ∆t/∆x)measured from a ZOS having
the following meaning
x . . . inline coordinate at which the measurement is made.
t0 . . . normal-incidence traveltime (two way) measured at x.
∆t/∆x . . . horizontal gradient of normal-incidence traveltime measured at (x, t0 ).
To perform a raytrace migration, the following materials are required:
1. (xi , t0ij , ∆t/∆xij ) . . . a set of picks to be migrated. The picks are made at a the locations xi
and at each location a number of t0ij and ∆t/∆xij values are measured.
2. v(x, z) . . . a velocity model is required. It must be defined for the entire range of coordinates
that will be encountered by the rays.
3. Computation aids. . . one or more of: a calculator, compass, protractor, computer.
Finally, the algorithm is presented as a series of steps with reference to Figure 4.9. For all
locations xi , i = 1, 2, . . . n then each ray must be taken through the following steps:
Step 1 . . . Determine the emergence angle at the surface of the ijth ray according to the formula
v0 (xi ) ∆t
sin θ0ij = (4.28)
2 ∆x ij
Step 2 . . . Denote the current layer by the number h. Project the ray down to the next interface
and determine the incident angle φhij .
4.2. FUNDAMENTAL MIGRATION CONCEPTS 89
Step 3 . . . Determine the length of the ijth raypath in the current (hth ) layer, lhij and compute the
layer traveltime
lhij
δthij = (4.29)
vh
Step 4 . . . Determine the refraction angle into the (h + 1)th layer from Snell’s law:
vh+1
sin(θ(h+1)ij ) = sin(φhij ) (4.30)
vh
Step 5 . . . Repeat steps 2 → 4 for each velocity layer until the stopping condition is fulfilled
n
t0ij
tij = δthij = (4.31)
2
h=1
which simply says that the raypath traveltime must be half the measured traveltime.
Step 6 . . . Draw in a reflecting segment perpendicular to the raypath at the point that the stopping
condition is satisfied.
least-time
crest ray
ray
vslow
vfast
true distorted
image
Figure 4.10: The least-time ray does not always come from the crest of an anticline.
a A
b B
sin a sin b sin A sin B
v = v v = v
Figure 4.11: (Left) When a ray encounters a dipping velocity interface, Snell’s law predicts
the transmitted ray using a relation between angles with-respect-to the interface normal.
(Right) In a time migration algorithm, Snell’s law is systematically violated because the
vertical angles are used. Both techniques give the same result if the velocity interface is
horizontal.
point of the least-time ray. Generally, a traveltime crest will have zero ∆t/∆x which means that
the corresponding ray emerges vertically. It turned out that the migration schemes were producing
a distorted image in cases like this and the distortion placed the crest of the migrated anticline at
the emergence point of the least time ray. Just like in the sedimentary basin case, the algorithm was
leaving flat events alone even though, in this case, it should not.
The understanding and resolution of this problem resulted in the terms time migration and depth
migration. The former refers to any method that has the bias towards flat events and that works
correctly in v(z) settings. The latter refers to newer methods that were developed to overcome these
problems and therefore that can produce correct images even in strong lateral velocity gradients.
Robinson (1983) presented the raytrace migration analogs to both time and depth migration and
these provide a great deal of insight. The raytrace migration algorithm presented previously in
section 4.2.2 is a depth migration algorithm because it obeys Snell’s law even when ∂x v is significant.
Figure 4.11 shows how Snell’s law must be systematically altered if a time migration is desired.
Instead of using the angles that the incident and transmitted rays make with respect to the normal
4.2. FUNDAMENTAL MIGRATION CONCEPTS 91
in Snell’s law, the time migration algorithm uses the angles that the rays make with respect to the
vertical. It is as though the dipping interface is rotated locally to the horizontal at the incident
point of the ray. Thus, any vertical ray (i.e. ∆t/∆x = 0) will pass through the velocity interface
without deflection regardless of the magnitude of the velocity contrast.
This explains why time migration gets the wrong answer and it also explains why the technology
remains popular despite this now well-understood shortcoming. The fact that flat events are unaf-
fected by time migration regardless of the velocity model means that the process has a “forgiving”
behavior when the velocity model has errors. Velocity models derived automatically from stacking
velocities tend to be full of inaccuracies. When used with time migration, the inaccuracies have
little effect but when used with depth migration they can be disastrous. Even today, the majority
of seismic sections come from v(z) settings and they consist of mostly flat events. Such data can be
migrated with relatively little care using time migration and well-focused sections are obtained. Of
course, these are still time sections and a pleasing time migration does not mean that the migration
velocities can be used for depth conversion. In general, the velocity errors that are “ignored” by
time migration will cause major problems in a subsequent depth conversion. Much more care must
be taken if an accurate depth section is required.
The terms time migration and depth migration arose, in part, because the early versions of
both technologies tended to produce only time and depth sections respectively. However, this is no
longer the case and either technology can produce both time and depth sections. Therefore, the
terms should be treated simply as jargon that indicates whether or not an algorithm is capable of
producing a correct image in the presence of strong lateral velocity gradients. The mere presence of
a time or depth axis on a migrated seismic display says nothing about which technology was used.
In this book, time migration will be used to refer to any migration technique that is strictly valid
only for v(z) while depth migration will refer to a technique that is valid for v(x, z). Thus, when
employed in constant velocity or v(z) settings, both techniques are valid and there is no meaningful
distinction. Under these conditions, a migrated time section can always be converted to a depth
section (or vice-versa) with complete fidelity.
The first time-migration algorithms were almost always finite-difference methods. Today, a better
understanding of the subject allows virtually any migration method to be recast as a time migration.
Therefore, it is very important to have some knowledge of the algorithm, or to have a statement
from the developer, to be sure of the nature of a method. The issue has been clouded even further
by the emergence of intermediate techniques. Sometimes, it is worth the effort to build a synthetic
seismic section to test a migration algorithm whose abilities are uncertain. One thing remains clear,
depth migrations always require much more effort for a successful result. Not only is more computer
time required; but, much more significantly, vastly greater human time may be necessary to build
the velocity model.
A B
x x x x x x x
x x x x x
z d t
l t t t t t
l
l
l
l
δ
x C
x x x x x
za
l l l l l
Figure 4.12: (A) Normal-incidence raypaths are shown for a dipping reflector in a constant
velocity medium. (B) The traveltime section for (A) plots the arrival times for each raypath
beneath its emergence point. (C) When (B) is plotted at apparent depth, za = vt/2, an
apparent dip can be defined.
4.2. FUNDAMENTAL MIGRATION CONCEPTS 93
The reflector dip, δ, can be related to the lengths of the raypaths and their emergence points by
l5 − l1
sin δ = . (4.32)
x5 − x1
The time section corresponding to this raypath diagram is shown in Figure 4.12B. Since velocity is
constant, the arrival times are directly proportional to the path lengths through tk = 2lk /v and the
ZOS image of the reflector is seen to have a time dip of
∆t t5 − t1 2 l5 − l1 sin δ
= = =2 (4.33)
∆x x5 − x1 v x5 − x1 v
which is in agreement with equation (4.27) that the time dip on a ZOS gives twice the normal-
incidence ray parameter. A further conceptual step can be taken by mapping the traveltimes to
apparent depth, za , by za = vt/2 as shown in Figure 4.12C. At this stage, the unmigrated display
has both coordinate axes in distance units so that it makes sense to define an apparent dip through
l 5 − l1
tan α = . (4.34)
x5 − x1
Finally, comparison of equations (4.32) and (4.34) gives the relation
Equation (4.35) is called the migrator’s equation and is of conceptual value because is illustrates
the geometric link between unmigrated and migrated dips. However, it should not be taken too far.
It is not generally valid for non-constant velocity and cannot be applied to a time section. Though
this last point may seem obvious, it is often overlooked. It is not meaningful to talk about a dip
in degrees for an event measured on a time section because the axes have different units. A change
of time-scale changes the apparent dip. On an apparent depth section, an apparent dip can be
meaningfully defined and related to the true dip through equation (4.35). The apparent dip can
never be larger than 45◦ because sin δ ≤ 1 which is a re-statement of the fact that the maximum
time dip is ∆t/∆x = 1/v as discussed in section 2.11.1. Though intuitively useful, the apparent dip
concept is not easily extended to variable velocity. In contrast, the limitation on maximum time dip
is simply restated in v(z) media as ∆t/∆x = 1/vmin where vmin is the minimum of v(z).
Comparison of Figures 4.12A and 4.12C shows that, for a given surface location xk , the normal-
incidence reflection point and the corresponding point on the ZOS image both lie on a circle of radius
lk centered at xk . This relationship is shown in Figure 4.13 from which it is also apparent that
• The ZOS image lies at the nadir (lowest point) of the circle.
• The ZOS image and the reflector coincide on the datum (recording plane).
These observations form the basis of a simple wavefront migration technique which is kinemati-
cally exact for constant velocity:
ii) Replace each point in the (x, za ) picture with a wavefront circle of radius za . Amplitudes along
the wavefront circle are constant and equal to the amplitude at (x, za ).
94 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Datum xk Datum
δ α
R
ef
le
lk
ct
R
ef
or
le
ct
fo
lk
or
rm
s
Zo
ta
s
ng
im
ag
en
e
t
to
Zo
w
si
av
ma
ef
ge
ro
nt
s
Figure 4.13: The normal-incidence reflection Figure 4.14: The wavefront migration of a dip-
point and the corresponding point on the ZOS ping reflector is shown. Each point on the ZOS
image both lie on a circle of radius lk centered image is replaced by a wavefront circle of radius
at xk . lk . The wavefronts interfere constructively at
the location of the actual reflector.
iii) The superposition of all such wavefronts will form the desired depth section.
Figure 4.14 shows the migration of a dipping reflector by wavefront superposition. The actual
mechanics of the construction are very similar to a convolution. Each point in the ZOS image is used
to scale a wavefront circle which then replaces the point. A 2D convolution can be constructed by
a similar process of replacement. The difference is that the replacement curve in a 2D convolution
is the same for all points while in the migration procedure it varies with depth. Thus the migration
process can be viewed as a nonstationary convolution.
Wavefront migration shows that the impulse response of migration is a wavefront circle. That is,
if the input ZOS contains only a single non-zero point, then the migrated section will be a wavefront
circle. The circle is the locus of all points in (x, z) that have the same traveltime to (xk , 0). Two
equivalent interpretations of this are either that the circle is the only earth model that could produce
a single output point or that the circle is the curve of equal probability. The first interpretation is
a deterministic statement that only one earth model could produce a ZOS image with a single
live sample. The second interpretation is a stochastic one that suggests why the generalization
from one live sample to many works. By replacing each point with it’s curve of equal probability,
the migrated section emerges as the most likely geology for the superposition of all points. The
constructive interference of a great many wavefronts forms the migrated image.
Exercise 4.2.1. Describe the appearance of:
• a migrated section that contains a high amplitude noise burst on a single trace.
• a migration that results from a ZOS image whose noise level increases with time.
• the edge of a migrated section.
Datum
v
v
point sources
v
t v
Point diffractors
(scatterers)
v
t + ∆t Reef
Figure 4.15: Huygens’ principle reconstructs the Figure 4.16: The seismic response of a contin-
wavefront at t + ∆t by the superposition of the uous geological structure is synthesized as the
small wavefronts from secondary point sources superposition of the responses of many point
placed on the wavefront at t. diffractors.
θ is the”migration dip”
point diffractors
reflector
A B
(scatterers)
Figure 4.17: Each point on a reflector is consid- Figure 4.18: (A) The ray paths for the left side of
ered to be a point diffractor that scatters energy a zero-offset point diffractor. (B) The raypaths
to all surface locations. for a CMP gather assuming a zero-dip reflector.
These raypaths are identical for v(z) leading to
the conclusion that the traveltime curves for the
point diffractor are the same as for the NMO
experiment.
0 0
A B
0.5 0.5
seconds
1 1
1.5 1.5
-3000 -2000 -1000 0 1000 2000 3000 -3000 -2000 -1000 0 1000 2000 3000
meters meters
Figure 4.19: (A) This diffraction chart shows the traveltime curves for five point diffractors
in a constant velocity medium. All five curves are asymptotic to the line x=vt/2. (B) The
second hyperbola has been migrated by wavefront superposition.
4.2. FUNDAMENTAL MIGRATION CONCEPTS 97
Zo
s
im
500 ag
e
Di
pp
met er s
ing
ref
lec
tor
1000
1500
-1000 -500 0 500 1000
met er s
Figure 4.20: The ZOS image of the dipping reflector is formed by replacing each point on
the reflector with the ZOS image of a point diffractor. The superposition of these many
hyperbolae forms the dipping reflector’s ZOS image.
that the Dix equation proof that the NMO hyperbola is approximately characterized by vrms means
that the diffraction traveltimes are approximately
4(x − x0 )2
t2x ≈ t20 + 2
(4.37)
vrms
where t0 = 2z/vave .
Figure 4.14 shows how wavefront superposition migrates the ZOS image of a dipping reflector.
The ZOS image of a single point diffractor is the hyperbola given by equation (4.36) so it is instructive
to see how wavefront migration can collapse such hyperbolae. Figure 4.19A shows a diffraction
chart that gives the traveltime curves for five point diffractors in a constant velocity medium. The
wavefront migration of this chart should convert it into five impulses, one at the apex of each
hyperbola. Figure 4.19B shows the wavefront migration of the second hyperbola on the chart. Each
point on the second hyperbola has been replaced by a wavefront circle whose radius is the vertical
time of the point. Since the chart is drawn with a vertical coordinate of time rather than depth,
these wavefront curves may not appear to be exactly circular. The geometry of hyperbola and
circles is such that the wavefront circles all pass through the apex of the second hyperbola. Thus
the amplitude of the wavefront superposition will be large at the apex and small elsewhere. The
same process will also focus all other hyperbolae on the chart simultaneously.
The diffraction curves can be used to perform the inverse of migration or modelling. Figure 4.20
shows the construction of the ZOS image of a dipping reflector using Huygens’ principle. Each point
on the dipping reflector is replaced by the ZOS image of a point diffractor. The superposition of
these many hyperbolae constructs the ZOS image of the dipping reflector. For a given hyperbola,
its apex lies on the geology and it is tangent to the ZOS image. This directly illustrates that the
wavefront migration of dipping reflectors is completely equivalent to the migration of diffraction
responses. In fact, the collapse of diffraction responses is a complete statement of the migration
problem and a given algorithm can be tested by how well it achieves this goal.
98 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Receivers
Receivers
.
d
−
ix
feπ
k
φ
∞
2
1
,t)=
V (x
0
ψ
Wavefield snapshot
Wavefield at time 0
Reflector
Figure 4.21: The exploding reflector model at Figure 4.22: A snapshot of the ERM wavefield
the instant of explosion establishes a wavefield at some instant of time as it approaches the re-
that is isomorphic with the reflector. ceivers
(x − x0 )2
t2x = t20 + (4.38)
v̂ 2
where v̂ = v/2 is called the exploding reflector velocity and t0 = z/v̂. This trivial recasting allows
the interpretation that the point diffractor is a seismic source of waves that travel at one-half of
the physical velocity. As shown in Figure 4.21 the exploding reflector model adopts this view and
postulates a model identical to the real earth in all respects except that the reflectors are primed with
explosives and the seismic wave speeds are all halved. Receivers are placed at the CMP locations
and at t = 0 the explosives are detonated. This creates a wavefield that is morphologically identical
to the geology at the instant of explosion.
If ψ(x, z, t) denotes the exploding reflector wavefield, then the mathematical goal of the migration
problem is to calculate ψ(x, z, t = 0). Thus ψ(x, z, t = 0) is identified with the geology and represents
the migrated depth section. The ERM wavefield is allowed to propagate only upward (the −z
direction) without reflections, mode conversions, or multiples. It refracts according to Snell’s law and
the half-velocities mean that the traveltimes measured at the receivers are identical to those of the
ZOS. In the constant velocity simulation of Figure 4.22, a snapshot of the ERM wavefield, ψ(x, z, t =
4.2. FUNDAMENTAL MIGRATION CONCEPTS 99
receivers
ge
ma
osWi gy
z lo
eo
G
Figure 4.23: The ERM seismogram is shown superimposed on the migrated depth section.
const), is shown as the wavefield approaches the receivers. The raypaths are normal-incidence
raypaths (because of the isomorphism between the geology and ψ(x, z, t = 0)) and spreading and
buried foci are manifest. The figure emphasizes three Huygen’s wavelets that emanate from the three
corners of the geologic model. The migrated depth section is also a snapshot, but a very particular
one, so the term snapshot will be used to denote all such constant-time views.
In Figure 4.23, the ERM seismogram, ψ(x, z = 0, t), is shown in apparent depth superimposed
on top of the depth section. Wavefront circles are shown connecting points on the geology with
points on the seismogram. The ERM seismogram is kinematically identical (i.e. the traveltimes are
the same) with the normal-incidence seismogram and is thus a model of the CMP stack. It is also
kinematically identical with all normal-incidence, primary, events on a ZOS image. This allows the
migration problem to be formulated as a solution to the wave equation. Given ψ(x, z = 0, t) as a
boundary condition, a migration algorithm solves the wave equation for the entire ERM wavefield,
ψ(x, z, t), and then sets t = 0 to obtain the migrated depth section. This last step of setting t = 0
is sometimes called imaging though this term has also come to refer to the broader context of
migration in general. Precisely how the wavefield is evaluated to extract the geology is called an
imaging condition. In the post-stack case the imaging condition is a simple evaluation at t = 0.
Later, it will be seen that there are other imaging conditions for pre-stack migration.
Thus far, the ERM has given simple mathematical definitions to the migrated depth section, the
recorded seismogram, and the wavefield snapshot. In addition, the extrapolated seismogram can be
defined as ψ(x, z = ∆z, t). Both the ERM seismogram and the extrapolated seismogram are time
sections while the snapshot and the migrated section are in depth. The extrapolated seismogram
is a mathematical simulation of what would have been recorded had the receivers been at z = ∆z
rather than at z = 0. The construction of ψ(x, z = ∆z, t) from ψ(x, z = 0, t) is called downward
continuation and, synonymously, wavefield extrapolation.
As a summary, the ERM has defined the following quantities:
cmp stack
tN
z x
extrapolated
∆z section
migrated depth
section
Figure 4.24: This prism shows the portion of (x, z, t) space that can be constructed from
a finite-extent measurement of ψ(x, z = 0, t) (upper face). Also shown are the migrated
depth section, ψ(x, z, t = 0), (vertical slice) and an extrapolated section, ψ(x, z = ∆z, t),
(horizontal slice).
These quantities are depicted in Figure 4.24. It is significant that any extrapolated seismogram can
be evaluated at t = 0 to give a single depth sample of the migrated depth section. The process
of deducing a succession of extrapolated seismograms, and evaluating each at t = 0, is called a
recursive migration. It is recursive because the extrapolated seismogram ψ(x, z = zk , t) is computed
from the previous seismogram ψ(x, z = zk−1 , t). Examples of recursive migrations are the finite-
difference methods and the phase-shift techniques. An alternative to the recursive approach is a
direct migration that computes ψ(x, z, t = 0) directly from ψ(x, z = 0, t) without the construction
of any intermediate products. Examples of direct migration are (kx , f ) migration and Kirchhoff
migration.
A direct migration tends to be computationally more efficient than a recursive approach both
in terms of memory usage and computation speed. However, the direct methods must deal with
the entire complexity of the velocity structure all at once. In contrast, the recursive approach
forms a natural partitioning of the migration process into a series of wavefield extrapolations. The
extrapolation from z = zk−1 to z = zk need only deal with the velocity complexities between these
two depths. If the interval zk−1 → zk is taken sufficiently small, then the vertical variation of velocity
can be ignored and v can be considered to depend only upon the lateral coordinates. In this way
the migration for a complex v(x, z) variation can be built from the solution of many v(x) problems.
Code Snippet 4.3.2. This is similar to Code Snippet 4.3.1 but differs by using event diph to
create the dipping event. The result is shown in Figure 4.26.
1 v=2000;dx=10;dt=.004;%basic model parameters
2 x=0:dx:2000;%x axis
3 t=0:dt:2;%t axis
4 seis=zeros(length(t),length(x));%allocate seismic matrix
5 seis=event_hyp(seis,t,x,.4,700,v,1,3);%hyperbolic event
6 seis=event_dip(seis,t,x,[.75 1.23],[700 1500],1);%linear event
7 [w,tw]=ricker(dt,40,.2);%make ricker wavelet
8 seis=sectconv(seis,t,w,tw);%apply wavelet
102 CHAPTER 4. ELEMENTARY MIGRATION METHODS
meters meters
0 500 1000 1500 2000 0 500 1000 1500 2000
0 0
0.5 0.5
seconds
seconds
1 1
1.5 1.5
2 2
Figure 4.25: The hyperbolic response of a point Figure 4.26: Similar to Figure 4.25 except that
diffractor and a simple linear event are shown. the linear event was created by the superposition
The display is slightly clipped to make the of many hyperbolae. The display is clipped to
diffraction more visible. This was created by match that of Figure 4.25. See Code Snippet
Code Snippet 4.3.1. 4.3.2.
End Code
Code Snippet 4.3.1 illustrates the use of event hyp and event dip to create the model section
shown in Figure 4.25. The first three lines establish the basic model geometry and define the
coordinate axes; then, line 4 initializes the seismic matrix to zero. On line 5, event hyp is invoked
to make the hyperbolic diffraction response shown in the top left of Figure 4.25. The seismic matrix
is input to event hyp and is replaced by the output. The coordinate vectors (second and third input
parameters) are required to defined the geometry. The fourth and fifth input parameters are the
(x, t) coordinates of the apex of the diffraction and the sixth parameter is the velocity. The seventh
parameter sets the amplitude at the apex of the hyperbola and the eighth is a flag that determines
how amplitude decays down the limbs of the hyperbola. There are four possibilities: (1) no amplitude
decay, (2) decay given by a(x) = a(0)t0 /tx (where a(x) is the amplitude at offset x, t0 is the zero-
offset traveltime, and tx is the traveltime at offset x), (3) decay given by a(x) = a(0)(t0 /tx )3/2 , and
(4) decay given by a(x) = a(0)(t0 /tx )2 . Line 6 invokes event dip to create the linear dipping event
seen in the center of Figure 4.25. The fourth input parameter specifies the time of the event at
its beginning and end while the fifth parameter gives the corresponding lateral positions. The last
parameter gives the event amplitude. Once the events are created, the final two lines make a Ricker
wavelet (dominant frequency of 40 Hz) and convolve it with the seismic section.
The linear event shown in Figure 4.25 is physically unrealistic because it lacks diffractions from
its endpoints. Real wavefields are never discontinuous. A much more realistic event can by created
using event diph . This function synthesizes a linear event by hyperbolic superposition exactly
as illustrated in Figure 4.20. The result is the seismic section shown in Figure 4.26. In order to
do so, the input parameters must describe the dipping event in (x, z) rather than (x, t). Code
Snippet 4.3.2 shows how this is done and differs from Code Snippet 4.3.1 only on line 6. The fourth
input parameter for event diph is the velocity while the next four parameters are respectively: the
starting and ending x coordinates, the starting z coordinate, and the dip in degrees. The lateral
extent of the resulting event is generally much greater than the prescribed coordinates. If the section
4.3. MATLAB FACILITIES FOR SIMPLE MODELLING AND RAYTRACE MIGRATION 103
meters meters
0 500 1000 1500 2000 0 500 1000 1500 2000
0 0
0.5 0.5
seconds
seconds
1 1
1.5 1.5
2 2
Figure 4.27: The response of five reflectors with Figure 4.28: The same dataset as Figure 4.28 is
dips of 0◦ , 20◦ , 40◦ , 60◦ , and 80◦ . shown with a strongly clipped display to reveal
the underlying hyperbolae.
is migrated, then the migrated reflector will be found between the specified coordinates.
A more complex seismic section is shown in Figure 4.27. The geologic model is a fan of five
dipping reflectors that originate at the point (x, z) = (250m, 200m) and extend to the right until
x = 1500m. In Figure 4.28 the same seismic section is shown with a strongly clipped display to
reveal the underlying hyperbolae.
Code Snippet 4.3.3. This code illustrates the use of event pwlinh to create a simple model of a
reef. The result is shown in Figure 4.29.
End Code
Code Snippet 4.3.3 illustrates the use of these hyperbolic summation tools to create a simple
model of a reef. The resulting seismic response is shown in Figures 4.29 and 4.30. The reef is a
simple trapezoidal structure, 400 m wide and 50 m high, on top of a flat reflector 600 m below the
recording plane. The reflection coefficient of the reef is modelled as -.1 (the acoustic impedance
within the reef is assumed to be lower than the surrounding material), +.1 on the base reflector,
and +.2 beneath the reef.
104 CHAPTER 4. ELEMENTARY MIGRATION METHODS
meters meters
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
0 0
0.5 0.5
seconds
seconds
1 1
1.5 1.5
Figure 4.29: The seismic response of a simple Figure 4.30: The same dataset as Figure 4.29 is
reef model. This was created with Code Snippet shown with a strongly clipped display to reveal
4.3.3 the underlying hyperbolae.
Code Snippet 4.3.4. This code uses event diph2 to construct the response of a trapezoidal struc-
ture. The parameter ndelx (line 5) controls the spacing of the hyperbolae. Figures 4.31, 4.32, 4.33
and 4.34 correspond to values of ndelx of 15, 10, 5, and 1 respectively.
End Code
The function event diph2 is similar to event diph except that it allows control over the spacing
of hyperbolae through the input parameter ndelx. In event diph the hyperbolae spacing is never
greater than the grid spacing while in event diph2 the spacing is ndelx times greater than in
event diph . Code Snippet 4.3.4 illustrates the use of this function to create a series of figures
illustrating the gradual formation of the seismic response of a trapezoidal structure. The sequence
of Figures 4.31, 4.32, 4.33 and 4.34 shows the complete seismic response gradually forming as the
density of hyperbolae increases.
Also present in the synsections toolbox are two functions that make “standard” synthetics
called makestdsyn and makestdsynh . These differ in that one uses hyperbolic superposition for
linear events and the other does not. Their use can be illustrated in the script makesections that
can be run as a demonstration.
4.3. MATLAB FACILITIES FOR SIMPLE MODELLING AND RAYTRACE MIGRATION 105
0 0
0.5 0.5
seconds
seconds
1 1
1.5 1.5
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
meters meters
Figure 4.31: A trapezoidal structure is modelled Figure 4.32: Similar to Figure 4.31 except that
with only a few sparse hyperbolae. Relative to every 10th hyperbola is shown.
Figure 4.34 only every 15th hyperbola is shown.
0 0
0.5 0.5
seconds
seconds
1 1
1.5 1.5
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
meters meters
Figure 4.33: Similar to Figure 4.31 except that Figure 4.34: A trapezoidal structure is com-
every 5th hyperbola is shown. pletely modelled by hyperbolic superposition.
Compare with Figures 4.31, 4.32, and 4.33.
106 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0
0
0.1
200 0.2
0.3
400 0.4
seconds
meters
0.5
600 0.6
0.7
800 0.8
0.9
1000
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.35: The velocity model created by Code Figure 4.36: The fourth-order exploding reflec-
Snippet 4.3.5. tor seismogram created on line 35 of Code Snip-
pet 4.3.5.
30
31 %create the exploading reflector model
32 dt=.004; %temporal sample rate
33 dtstep=.001; %modelling step size
34 tmax=2*zmax/vlow; %maximum time
35 [seisfilt,seis,t,xrec]=afd_explode(dx,dtstep,-dt,tmax, ...
36 vel,length(x),x,zeros(size(x)),[10 15 40 50],0,2,0);
End Code
Code Snippet 4.3.5 illustrates the creation of a model of a small channel beneath a layered
medium. The creation of the velocity model uses afd vmodel and was discussed in detail in section
3.3. The channel is defined on lines 25-29 as 20 m wide and 50 m deep. With the 5 m grid (line 1)
this means that the channel is 4 samples wide and 10 samples deep. The resulting velocity model is
shown in Figure 4.35. Lines 32-36 take this velocity model into afd explode to create an exploding
reflector seismogram using the fourth-order Laplacian approximation (equation (3.20)). As with
afd shotrec , two temporal sample rates are prescribed. The coarser one (dt on line 32) specifies
the sample rate of the final seismogram while the finer one (dtstep on line 33) is the time step of
the modelling. The resulting seismogram, filtered with Ormsby parameters of [10 15 40 50] is shown
in Figure 4.36. Two small horizontal lines are superimposed to indicate the position of the channel.
The lines are at the time of the top and bottom of the channel (.312 and .355 seconds) and have the
same width as the channel. The channel produces an extensive diffraction pattern showing scattered
energy over the entire model.
Figure 4.37 shows a second exploding reflector seismogram of the channel model that was cal-
culated with the second-order approximate Laplacian (equation (3.19)). In comparison with the
fourth-order seismogram of Figure 4.36, this result is considerably less accurate though it required
only one-half of the computation time. Even less accurate is the second-order solution shown in
Figure 4.38 that was done with a grid sample size of 10 m rather than the 5m of the previous two
figures. This result is almost unrecognizable as the response of the same structure. An obvious
conclusion is that the results from finite difference modelling should not be trusted without a careful
consideration of the algorithms and parameters.
108 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0 0
0.1 0.1
0.2 0.2
0.3 0.3
0.4 0.4
seconds
seconds
0.5 0.5
0.6 0.6
0.7 0.7
0.8 0.8
0.9 0.9
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.37: Similar to the seismogram of Fig- Figure 4.38: Similar to the seismogram of Figure
ure 4.36 except that the Laplacian was approx- 4.37 except that the spatial sample rate was 10m
imated with second-order finite differences. rather than 5m.
0 0
200 200
400 400
meters
meters
600 600
800 800
1000 1000
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.39: An anticline beneath a high veloc- Figure 4.40: The reflectivity, corresponding to
ity wedge. Black is 4000 m/s and white is 2000 Figure 4.39, is shown. This was calculated with
m/s. afd reflect .
4.3. MATLAB FACILITIES FOR SIMPLE MODELLING AND RAYTRACE MIGRATION 109
0 0
0.1 0.1
0.2 0.2
0.3 0.3
0.4 0.4
seconds
seconds
0.5 0.5
0.6 0.6
0.7 0.7
0.8 0.8
0.9 0.9
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.41: The fourth-order exploding reflec- Figure 4.42: Similar to Figure 4.41 except that
tor seismogram corresponding to the model of the spatial grid size was 5 m.
Figure 4.39. This model was created using a
spatial grid size of 10 m.
As a slightly more complex example that will be useful in examining depth migration, consider
the response of an anticline beneath a high-velocity wedge as shown in Figure 4.39. Code Snippet
4.3.6 creates this velocity model and the fourth-order exploding reflector seismogram shown in Figure
4.41. This result was created with a spatial grid size of 10 m. It might be tempting to trust this
result as the proper physical response of the model, but that could lead to erroneous conclusions.
For example, it might be concluded that the response of the anticline has a series of reverberations
following the primary response. However, before leaping to conclusions, it is prudent to recreate
the model with a finer spatial grid size. Figure 4.42 shows the result of re-running Code Snippet
4.3.6 with the spatial grid size set to 5 m. The reverberations have vanished and the amplitude
variation of the primary response of the anticline has changed considerably. Once again, it pays to
be cautious when working with finite difference synthetics. The only way to be certain is to fully
understand the algorithm and how to change the parameters to increase the accuracy. When the
accuracy parameters are increased by a factor of two and the response does not change significantly,
then it is safe to conclude that the model is showing an accurate physical response. Of course, this
can be a tedious process because increasing the accuracy will always increase the computation time.
Code Snippet 4.3.6. This creates an exploding reflector model of an anticline beneath a high
velocity wedge. Lines 1-25 build the velocity model and lines 27-34 create a fourth-order exploding
reflector seismogram. The results are shown in Figures 4.39, and 4.41.
1 dx=5; %cdp interval
2 xmax=2500;zmax=1000; %maximum line length and maximum depth
3 xpinch=1500; % wedge pinchout coordinates
4 zwedge=zmax/2; % wedge maximum depth
5 x=0:dx:xmax; % x coordinate vector
6 z=0:dx:zmax; % z coordinate vector
7 vhigh=4000;vlow=2000; % high and low velocities
8 vel=vlow*ones(length(z),length(x));%initialize velocity matrix
9 % define the wedge as a three point polygon
10 dx2=dx/2;
11 xpoly=[-dx2 xpinch -dx2];zpoly=[-1 -1 zwedge];
110 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0
0
0.1
0.2 200
0.3
0.4 400
seconds
meters
0.5
0.6 600
0.7
0.8 800
0.9
1000
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.43: The seismogram of Figure 4.42 is Figure 4.44: The picks of Figure 4.43 are shown
shown with picks on the image of the anticline. migrated on top of the reflectivity section.
with the command rayvelmod(vel,dx) where vel is the velocity matrix and dx is the grid size.
Then, the command eventraymig(figno), where figno is the MATLAB figure number of Figure
4.39, causes the picks to be migrated and the resulting raypaths displayed in Figure 4.44. At the
termination of each raypath, a small line segment, perpendicular to the raypath, is drawn that
indicates the implied reflector dip. (These do not appear perpendicular in Figure 4.44 because the
(x, z) axes do not have the same scale. The command axis equal will display any figure with equal
scales on all axes.)
The relative accuracy of the migrations in Figure 4.44 is instructive. The picks have all migrated
to positions near the anticline but some have fallen short while others have gone too far. Among the
obvious reasons for this are the difficulty in determining which phase (peak, trough, zero crossing,
etc.) of the input waveform should be picked and then making a consistent pick at two points. A
slight error in either case can result in a very large error in the final position. Since the material
beneath the anticline has a high velocity, a pick that arrives at the anticline with a little time left
(see section 4.2.2) will continue a significant distance. Thus, small errors in the pick time can make
a large error in the result. Also, an error in picking ∆t/∆x will cause the initial trajectory of the ray
to be incorrect. Since sin θ0 = .5v0 ∆t/∆x, these emergence angle errors are also more significant in
high-velocity material.
Migration by normal raytracing can reveal a lot about the nature of a seismic dataset. Also
instructive is the complementary process of normal-incidence raytrace modelling. This process is
implemented in the function normray that is logical reverse of normraymig . The latter requires
the pick specification of (x0 , t0 , ∆t/∆x) while the former needs the specification of the normal ray:
(xn , zn , θn ). Here, (xn , zn ) are the coordinates of the normal incidence point and θn is the structural
dip (in degrees) at the normal incidence point. Though logically similar, it is convenient to use
separate raytracing engines for these two tasks because they have different criteria for stopping
the ray. In migration, the ray is terminated when it has used the available traveltime while in
modelling, it is stopped when it encounters the recording surface (z = 0). These raytracing engines
are shootrayvxz and shootraytosurf respectively.
As with the migration tools, it is tedious to invoke normray at the command line for each pick.
Therefore, a convenience function, eventraymod is provided that automatically models any picks
found in the global variable PICKS. These picks are expected to have been made on a depth section
112 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0
0
0.1
200 0.2
0.3
400 0.4
seconds
meters
0.5
600 0.6
0.7
800 0.8
0.9
1000
0 500 1000 1500 2000 2500 0 500 1000 1500 2000 2500
meters meters
Figure 4.45: The reflectivity section of Figure Figure 4.46: The picks of Figure 4.45 are shown
4.40 is shown with picks on the anticline and modelled on top of the seismic section of Figure
normal rays to the surface. 4.42.
though no check is made to ensure this. Figure 4.45 shows the reflectivity section of Figures 4.40 and
4.44 with a series of picks, (xn , zn , θn ), made on the anticline. Also shown are the normal incidence
raypaths (drawn by normray ) to the surface. Figure 4.46 shows the modelled picks, (x0 , t0 , ∆t/∆x),
on top of the seismic section of Figures 4.42 and 4.43.
4.4.1 f -k migration
Stolt (1978) showed that the migration problem can be solved by Fourier transform. Here, Stolt’s
solution will be developed from a formal solution of the wave equation using Fourier transforms. It
will be developed in 2D and the direct 3D extension will be stated at the end.
Let ψ(x, z, t) be a scalar wavefield that is a solution to
1 ∂2ψ
∇2 ψ − =0 (4.40)
v̂ 2 ∂t2
where v̂ is the constant ERM velocity. It is desired to calculate ψ(x, z, t = 0) given knowledge of
ψ(x, z = 0, t). The wavefield can be written as an inverse Fourier transform of its (kx , f ) spectrum
as
ψ(x, z, t) = φ(kx , z, f )e2πi(kx x−f t) dkx df (4.41)
V∞
4.4. FOURIER METHODS 113
where cyclical wavenumbers and frequencies are used and the Fourier transform convention uses a +
sign in the complex exponential for spatial components and a − sign for temporal components. (The
notation for the integration domain is explained in section 4.1.3.) If equation (4.41) is substituted
into equation (4.40), the various partial derivatives can be immediately brought inside the integral
where they can be readily evaluated. The result is
2 2
∂ φ(z) 2 f 2
+ 4π − kx φ(z) e2πi(kx x−f t) dkx df = 0 (4.42)
∂z 2 v̂ 2
V∞
where the (kx , f ) dependence in φ(z) has been suppressed for simplicity of notation. The derivation of
equation (4.42) does not require that v̂ be constant; however, the next step does. If v̂ is constant1 ,
then the left-hand-side of equation (4.42) is the inverse Fourier transform of the term in curly
brackets. The uniqueness property of Fourier transforms (that there is a unique spectrum for a
given function and vice-versa) guarantees that is a function vanishes everywhere in one domain, it
must do so in another. Put another way, the zero function has a zero spectrum. Thus, it results
that
∂ 2 φ(z)
+ 4π 2 kz2 φ(z) = 0 (4.43)
∂z 2
where the wavenumber kz is defined by
f2
kz2 = − kx2 . (4.44)
v̂ 2
Equation (4.44) is called the dispersion relation for scalar waves though the phrase is somewhat
misleading since there is no dispersion in this case.
Equations (4.43) and (4.44) are a complete reformulation of the problem in the (kx , f ) domain.
The boundary condition is now φ(kx , z = 0, f ) which is the Fourier transform, over (x, t), of ψ(x, z =
0, t). Equation (4.43) is a second-order ordinary differential equation for fixed (kx , f ). Either of the
functions e±2πikz z solve it exactly as is easily verified by substitution. Thus the unique, general
solution can be written as
where A(kx , f ) and B(kx , f ) are arbitrary functions of (kx , f ) to be determined from the boundary
condition(s). The two terms on the right-hand-side of equation (4.45) have the interpretation of a
downgoing wavefield, A(kx , f )e2πikz z , and an upgoing wavefield, B(kx , f )e−2πikz z . This can be seen
by substituting equation (4.45) into equation (4.41) and determining the direction of motion of the
individual Fourier plane waves as described in section 2.9. It should be recalled that z increases
downward.
Given only one boundary condition, φ(kx , z = 0, f ), it is now apparent that this problem cannot
be solved unambiguously. It is a fundamental result from the theory of partial differential equations
that Cauchy boundary conditions (e.g. knowledge of both ψ and ∂z ψ) are required on an open
surface in order for the wave equation to have a unique solution. Since this is not the case here, the
migration problem is said to be ill-posed. If both conditions were available, A and B could be found
as the solutions to
φ(z = 0) ≡ φ0 = A + B (4.46)
1 Actually v̂(z) could be tolerated here. The necessary condition is that v̂ must not depend upon x or t.
114 CHAPTER 4. ELEMENTARY MIGRATION METHODS
and
∂φ
(z = 0) ≡ φz0 = 2πikz A − 2πikz B. (4.47)
∂z
When faced with the need to proceed to a solution despite the fact that the stated problem does
not have a unique solution, a common approach is to assume some limited model that removes the
ambiguity. The customary assumption of one-way waves achieves this end. That is, ψ(x, z, t) is
considered to contain only upgoing waves. This allows the solution
Then, the ERM wavefield can be expressed as the inverse Fourier transform
ψ(x, z, t) = φ0 (kx , f )e2πi(kx x−kz z−f t) dkx df. (4.49)
V∞
Equation (4.50) gives a migrated depth section as a double integration of φ0 (kx , f ) over f and
kx . Though formally complete, it has the disadvantage that only one of the integrations, that over
kx , is a Fourier transform that can be done rapidly as a numerical FFT. The f integration is not
a Fourier transform because the Fourier kernel e−2πf t was lost when the imaging condition (setting
t = 0) was invoked. Inspection of equation (4.50) shows that another complex exponential e−2πikz z
is available. Stolt (1978) suggested a change of variables from (kx , f ) to (kx , kz ) to obtain a result
in which both integrations are Fourier transforms. The change of variables is actually defined by
equation (4.44) that can be solved for f to give
f = v̂ kx2 + kz2 . (4.51)
Performing the change of variables from f to kz according to the rules of calculus transforms equation
(4.50) into
ψ(x, z, t = 0) = φm (kx , kz )e2πi(kx x−kz z) dkx dkz (4.52)
V∞
where
∂f (kz ) v̂kz
φm (kx , kz ) ≡ φ0 (kx , f (kz )) = φ0 (kx , f (kz )) (4.53)
∂kz kx2 + kz2
Equation (4.52) is Stolt’s expression for the migrated section and forms the basis for the f -
k migration algorithm. The change of variables has recast the algorithm into one that can be
accomplished with FFT’s doing all of the integrations. Equation (4.53) results from the change of
variables and is a prescription for the construction of the (kx , kz ) spectrum of the migrated section
from the (kx , f ) spectrum of the ERM seismogram.
Many of the important properties of post-stack migration can be discerned from Stolt’s result.
First, notice that kz defined though
f2
kz = − kx2 (4.54)
v̂ 2
4.4. FOURIER METHODS 115
0.1
0.08 wavelike
kz=.06
0.06 f/v=.06
0.06
f/v
kz=.04
0.04 f/v=.04
0.04
kz=.02 kz
0.02 f/v=.02
0.02
evanescent
0 0
-0.06 -0.02 0 0.02 0.06 -0.06 -0.02 0 0.02 0.06
kx kx
Figure 4.47: The space (f /v, kx ) is shown con- Figure 4.48: The space (kx , kz ) is shown with
toured with values of kz from equation (4.54). contours of f /v as given by equation (4.51).
The dashed lines are the boundary between the
wavelike and evanescent regions.
f/v
10 30 50 70
kx
kz
10
30
50
70
kx
Figure 4.49: The mapping from (kx , f ) (top) to (kx , kz ) (bottom) is shown. A line of
constant f is mapped, at constant kx , to a semi-circle in (kx , kz ). The compressed apparent
dip spectrum of (f /v, kx ) space unfolds into the uniform fan in (kx , kz ) space. The numbers
on each space indicate dips in degrees.
116 CHAPTER 4. ELEMENTARY MIGRATION METHODS
is real-valued when |f /kx | ≥ v̂ and is otherwise imaginary. Only for real kz will equation (4.45)
correspond to traveling waves in the positive and negative z directions. For complex kz e±2πikz z
becomes a real exponential that either grows or decays; however, on physical grounds, growing
exponentials must be excluded. Given a value for v̂ this dual behavior of kz divides (kx , f ) space into
two regions as was discussed from another perspective in section 2.11.1. The region where |f /kx | ≥ v̂
is called the wavelike or body wave regions and where |f /kx | < v̂ it is called the evanescent region.
Equation (4.54) is sometimes called a dispersion relation for one-way waves since the choice of the
plus sign in front of the square root generates upward traveling waves in e±2πikz while the minus
sign generates downward traveling waves.
The geometric relationships between the spaces of (f /v, kx ) and (kx , kz ) are shown from two
different perspectives in Figures 4.47 and 4.48. In (f /v, kx ) space, the lines of constant kz are
hyperbolae that are asymptotic to the dashed boundary between the wavelike and evanescent regions.
In (kx , kz ) space, the curves of constant f /v are semi-circles. At kx = 0, kz = f /v so these hyperbolae
and semi-circles intersect when the plots are superimposed.
The spectral mapping required in equation (4.53) is shown in Figure 4.49. The mapping takes a
constant f slice of (kx , f ) space to a semi-circle in (kx , kz ) space. Each point on the f slice maps at
constant kx which is directly down in the figure. It is completely equivalent to view the mapping as
a flattening of the kz hyperbolae of Figure 4.47. In this sense, it is conceptually similar to the NMO
removal in the time domain though here the samples being mapped are complex valued. That the
spectral mapping happens at constant kx is a mathematical consequence of the fact that kx is held
constant while the f integral in equation (4.50) is evaluated. Conceptually, it can also be viewed as
a consequence of the fact that the ERM seismogram and the migrated section must agree at z = 0
and t = 0.
On a numerical dataset, this spectral mapping is the major complexity of the Stolt algorithm.
Generally, it requires an interpolation in the (kx , f ) domain since the spectral values that map to
grid nodes in (kx , kz ) space cannot be expected to come from grid nodes in (kx , f ) space. In order to
achieve significant computation speed that is considered the strength of the Stolt algorithm, it turns
out that the interpolation must always be approximate. This causes artifacts in the final result.
This issue will be discussed in more detail in section 4.4.2.
The creation of the migrated spectrum also requires that the spectrum be scaled by v̂kz / kx2 + kz2
as it is mapped (Equation (4.53)). In the constant velocity medium of this theory, sin δ = v̂kx /f (δ
is the scattering angle) from which it follows that cos δ = v̂kz /f = kz / kx2 + kz2 . Thus the scaling
factor is proportional to cos δ and therefore ranges from unity to zero as δ goes from zero to 90◦ .
This scaling factor compensates for the “spectral compression” that is a theoretical expectation of
the ERM seismogram. Recall the migrator’s formulae (equation (4.35)) that relates apparent angles
in (f /v, kx ) space to real angles in (kx , kz ) space. If there is uniform power at all angles in (kx , kz )
space, then the migrator’s formula predicts a spectral compression in (f /v, kx ) (with consequent
increasing power) near 45◦ of apparent dip. As shown in Figure 4.49 it is as though (kx , kz ) space is
an oriental fan that has been folded to create (f /v, kx ) space. Migration must then unfold this fan.
f -k migration is called a steep-dip method because it works correctly for all scattering angles from
0◦ to 90◦ .
The f -k migration algorithm just described is limited to constant velocity though it is exact
in this case. Its use of Fourier transforms for all of the numerical integrations means that it is
computationally very efficient. Though it has been used for many years in practical applications
its restriction to constant velocity is often unacceptable so that it has gradually being replaced by
more flexible, though usually slower, methods. Today, one major virtue still remains and that is the
conceptual understanding it provides. The description of the construction of the migrated spectrum
will provide the basis for a realistic theory of seismic resolution in section 4.7.
4.4. FOURIER METHODS 117
0 0
200 200
400 400
600 600
800 800
meters
meters
1000 1000
1200 1200
1400 1400
1600 1600
1800 1800
2000 2000
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 200 400 600 800 1000 1200 1400 1600 1800 2000
meters meters
Figure 4.50: The result of the f -k migration of Figure 4.51: The result of the f -k migration of
the data of Figure 4.25. the data of Figure 4.26.
• Forward (kx , f ) transform. This is done on the unmigrated data after all pre-conditioning such
as spectral whitening, bandpass filtering, and noise reduction.
• Spectral mapping. The spectral mapping and scaling described by equation (4.53) requires a
careful interpolation process.
• Inverse (kx , f ) transform. This must be the exact inverse of the transform done in the first
step.
This process is formalized in the MATLAB function fkmig . This function has the external inter-
face: [seismig,tmig,xmig] = fkmig(seis,t,x,v,params). Here the first four input variables are
simply the seismic matrix, its time coordinate vector, its x coordinate vector, and a scalar giving the
velocity of migration. The final input variable is a vector of parameters that control various aspects
of the migration. This vector, params has length thirteen and affords control over the spatial and
temporal zero pads, the maximum dip to migrate, the maximum frequency to migrate, the type of
spectral interpolation, and the kind of progress messages that are written to the screen. Consult
the online help for fkmig for a complete description. Usually, the default actions for all of these
behaviors are acceptable and params can be omitted. Occasionally, it will be desired to program
one or more elements of params. This can be done while still defaulting the others by first creating
a vector of thirteen NaN’s (e.g. params=nan*1:13;) and then setting a few elements of params to
specific values.
As a first example, consider the migration of the synthetic sections shown in Figures 4.25 and 4.26.
This is very simply accomplished with the code in Code Snippet 4.4.1 with the results in Figures 4.50
and 4.51. These figures are displayed with slight clipping to make visible their more subtle details.
In Code Snippet 4.4.1, the variable seis refers to the data of Figure 4.25 while seis2 refers to Figure
4.26. Figure 4.25 contains the image of a dipping reflecting segment without diffractions while Figure
4.26 contains a similar image except that diffractions are present. The resulting migrations make
118 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0 0
200
400
600
500
800
meters
meters
1000
1200
1000
1400
1600
1800
2000 1500
0 200 400 600 800 1000 1200 1400 1600 1800 2000 0 500 1000 1500 2000 2500 3000
meters meters
Figure 4.52: The result of the f -k migration of Figure 4.53: The result of the f -k migration of
the data of Figure 4.27. the data of Figure 4.29.
clear the advantage (indeed, the necessity) of modelling with diffractions. Only when diffractions
are included are the edges of the reflecting segment imaged with crisp definition.
Code Snippet 4.4.1. This code uses fkmig to migrate the sections shown in Figures 4.25 and
4.26. The results are shown in Figures 4.50, and 4.51.
1 seismig=fkmig(seis,t,x,v);
2 plotimage(seismig,t*v/2,x);
3 xlabel(’meters’);ylabel(’meters’);
4 seismig2=fkmig(seis2,t,x,v);
5 plotimage(seismig2,t*v/2,x);
6 xlabel(’meters’);ylabel(’meters’);
End Code
As further examples, the f -k migrations of the data of Figures 4.27, 4.29, 4.31, and 4.34 can be
migrated in similar fashion. The results are shown in Figure 4.52, 4.53, 4.54, and 4.55.
It is apparent that fkmig creates high quality migrations of constant-velocity synthetics, but it
is instructive to examine the code to see precisely how this is done. Code Snippet 4.4.2 contains an
excerpt from the source code of fkmig that illustrates the steps of basic geophysical importance.
The forward f -k transform is accomplished on line 2 by calling fktran . Prior to this, the seismic
matrix, seis, has had zero pads attached in both x (column) and t (row). The variables tnew
and xnew are time coordinate and x coordinate vectors that describe the padded seismic matrix.
Similarly, nsampnew and ntrnew are the number of samples-per-trace and the number of traces after
padding. Function fktran is a geophysical wrapper around MATLAB’s built-in two-dimensional
FFT fft2. This means that, in addition to calling fft2 to compute the f -k spectrum, fkspec, of
seis, it does some useful things like calculating the frequency and wavenumber coordinate vectors
f and kx. Also, because the input matrix is real valued, it only returns the non-negative temporal
frequencies. This is possible because the negative frequencies can be deduced from the positive ones
by a symmetry argument. Technically, if φ(kx , f ) is the f -k transform of the real-valued wavefield
ψ(x, t) then it can be shown that φ(kx , f ) has the symmetry φ(kx , f ) = φ̄(−kx , −f ), where the
overbar indicates the complex conjugate. Both positive and negative wavenumbers are required
since, after the temporal Fourier transform, the matrix is no longer real valued. Working only with
4.4. FOURIER METHODS 119
0 0
500 500
meters
meters
1000 1000
1500 1500
0 500 1000 1500 2000 2500 3000 0 500 1000 1500 2000 2500 3000
meters meters
Figure 4.54: The result of the f -k migration of Figure 4.55: The result of the f -k migration of
the data of Figure 4.31. the data of Figure 4.34.
the non-negative temporal frequencies is more efficient because less memory and fewer calculations
are required. Also, it requires more care to formulate an algorithm correctly to process both positive
and negative frequencies. If the processed spectrum does not have the correct conjugate symmetry
then its inverse transform will will result in a complex-valued seismic matrix. It is easier to process
only the non-negative frequencies and calculate the negative ones as needed from the symmetry
condition. The function fktran has a companion inverse ifktran (invoked on line 34) that creates
the negative temporal frequencies from the non-negative ones and then calls fft2.
Continuing with Code Snippet 4.4.2, line 3 calculates the exploding reflector velocity that is
required in the theory. The major computation loop is over kx (the columns of fkspec) and extends
from line 9 to line 32. Each iteration through the loop converts one column of fkspec, representing
φ0 of equation (4.53), into a column vector representing φm . Prior to the loop, lines 5-6 compute
the vertical coordinate vector, kz, for the matrix representing φm (df is the frequency sample rate).
On line 7, kz2 is pre-computed as this is needed in every loop. The first action in the loop (line 11)
is to apply frequency and dip masks to the relevant column of fkspec and the result is stored in
the temporary vector tmp. The calculation of these masks is not shown in this example. They are
simply vectors of the same size as a column of fkspec whose entries range between zero and one.
The frequency mask, fmask, is pre-computed outside the loop but the dip mask, dipmask must be
recalculated with every iteration. These masks are designed to attenuate the frequencies and dips
that are not of interest. For example, it is customary to migrate only frequencies below a certain
maximum, fmax . Thus a simple fmask could could be unity for f ≤ fmax and zero for f > fmax .
However, from a signal processing perspective, such an abrupt cutoff is not desirable so a better
fmask might be
1, f ≤ fmax − fwid
!
fmask = 21 + 12 cos π f −fmax +fwid
fwid , fmax − fwid < f ≤ fmax . (4.55)
0, f > f
max
This defines a raised-cosine ramp of width fwid from unity at f = fmax − fwid to zero at f = fmax .
The dip mask must be recomputed at every loop iteration because the frequency corresponding to
120 CHAPTER 4. ELEMENTARY MIGRATION METHODS
a given dip changes as kx changes according to f = kx v̂/ sin θ. Like the frequency mask, the dip
mask should use a raised-cosine ramp from zero at θmax to unity at θmax − θwid . Thus dipmask
attenuates low frequencies while fmask attenuates high frequencies.
Next, on line 13 in Code Snippet 4.4.2, the frequencies that map to each output kz called
fmap, are calculated via equation (4.51). Since the spectral mapping generally lowers frequencies
(see Figure 4.49) many of these frequencies will, in general, be greater than the maximum desired
frequency to migrate. Thus, on line 14, MATLAB’s powerful find command is used to identify
the entries in the vector fmap that are lower than the maximum frequency to migrate, fmapmig.
(fmapmig corresponds to fmax in equation (4.55).) At this stage, the vector ind is a vector of
indices into fmap that points to those frequencies that will be migrated.
On line 16 the current column of fkspec is set to zero in preparation for the actual migration
that happens in the if-block from line 17 to line 28. Lines 18-25 compute the cosine scale factor that
occurs in equation (4.53) with a special case for f = 0 to avoid division by zero. Line 27 does the
actual spectral mapping and applies the cosine scale factor. The spectral mapping is accomplished by
a sinc-function interpolation, that is optimized for complex-valued spectra, by the function csinci .
On lines 29-31, a progress message is printed. Though not strictly necessary, this is considered
good form because there are few things more disconcerting than waiting for a computer to finish a
long calculation without any indication that progress is occurring. A message is printed only after
every kpflag wavenumbers have been processed. This is controllable through the input params
vector.
Finally, after the completion of the loop on line 34, the migrated spectrum is inverse transformed.
Generally, the zero pads are then removed though this is not shown.
Code Snippet 4.4.2. This is an excerpt from the code of fkmig that shows the steps of major
importance. The actual fkmig source code contains 319 lines that handle the many non-technical
details required for a useful code.
1 %forward fk transform
2 [fkspec,f,kx] = fktran(seis,tnew,xnew,nsampnew,ntrnew,0,0);
3 ve = v/2; %exploding reflector velocity
4 %compute kz
5 dkz= df/ve;
6 kz = ((0:length(f)-1)*dkz)’;
7 kz2=kz.ˆ2;
8 %now loop over wavenumbers
9 for j=1:length(kx)
10 % apply masks
11 tmp=fkspec(:,j).*fmask.*dipmask;
12 %compute f’s which map to kz
13 fmap = ve*sqrt(kx(j)ˆ2 + kz2);
14 ind=find(fmap<=fmaxmig);
15 %now map samples by interpolation
16 fkspec(:,j) = zeros(length(f),1); %initialize output spectrum to zero
17 if( ˜isempty(ind) )
18 %compute cosine scale factor
19 if( fmap(ind(1))==0)
20 scl=ones(size(ind));
21 li=length(ind);
22 scl(2:li)=ve*kz(ind(2:li))./fmap(ind(2:li));
23 else
24 scl=ve*kz(ind)./fmap(ind);
25 end
26 %complex sinc interpolation
27 fkspec(ind,j) = scl.*csinci(tmp,f,fmap(ind),[lsinc,ntable]);
4.4. FOURIER METHODS 121
0 0
20 0.02
40 0.04
1/meters
60 0.06
Hz
80 0.08
100 0.1
120 0.12
−0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08 −0.1 −0.08 −0.06 −0.04 −0.02 0 0.02 0.04 0.06 0.08
1/meters 1/meters
Figure 4.56: The f -k spectrum of the data of Figure 4.57: The (kx , kz ) spectrum of the data
Figure 4.31. of Figure 4.54.
28 end
29 if( floor(j/kpflag)*kpflag == j)
30 disp([’finished wavenumber ’ int2str(j)]);
31 end
32 end
33 %inverse transform
34 [seismig,tmig,xmig]=ifktran(fkspec,f,kx);
End Code
As a last example, consider the computation of the discrete f -k spectra of one of the preceding
examples before and after migration. This should produce a graphical confirmation of the mapping
of Figure 4.49. This is very simply done using fktran . Specifically, the case of the synthetic of
Figure 4.31 and its migration in Figure 4.54 is shown. If seis is the unmigrated seismic matrix,
then the command [fkspec,f,kx]=fktran(seis,t,x) computes the complex-valued f -k spectrum
and plotimage(abs(fkspec),f,kx) produces the result shown in Figure 4.56. In a similar manner,
the (kx , kz ) spectrum after migration can be computed and is shown in Figure 4.57. Comparison of
these figures shows that the spectral mapping has been performed as described.
Considering the first velocity layer only, this result is valid for any depth within that layer provided
that v̂ is replaced by the first layer velocity, v̂1 . If the thickness of the first layer is ∆z1 , then the
ERM wavefield just above the interface between layers 1 and 2 can be written as
ψ(x, z = ∆z1 , t) = φ0 (kx , f )e2πi(kx x−kz1 ∆z1 −f t) dkx df. (4.56)
V∞
where
f2
kz1 = − kx2 . (4.57)
v̂12
Equation (4.56) is an expression for downward continuation or extrapolation of the ERM wavefield
to the depth ∆z1 . The extrapolated wavefield is distinguished from the surface recorded wavefield
by the presence of the term e2πikz1 ∆z1 under the integral sign that is a specific form of the Fourier
extrapolation operator, e2πikz ∆z . Any extrapolated wavefield is a temporal seismogram more akin to a
surface recording than to a migrated depth section. In the phase-shift method, as with any recursive
technique, the migrated depth section is built little-by-little from each extrapolated seismogram.
The extrapolated wavefield is a simulation of what would have been recorded had the receivers
actually been at z = ∆z rather than z = 0 . Since any extrapolated section intersects the depth
section at (z = ∆z, t = 0) each extrapolation can contribute one depth sample to the migration (see
Figure 4.24). This process of evaluating the extrapolated section at t=0 was discussed in section
4.2.6 as the post-stack imaging condition.
For the wavelike portion of (kx , f ) space, the extrapolation operator has unit amplitude and a
phase of 2πkz ∆z. For evanescent spectral components, it is a real exponential e±2π|kz |∆z . Forward
wavefield propagation must obey physical law and propagate evanescent spectral components using
the minus sign in the exponent, e.g. e−2π|kz |∆z . Therefore, inverse wavefield extrapolation, as is done
for migration, should use e+2π|kz |∆z . However, this inversion of evanescent spectral components is a
practical impossibility because they have decayed far below the noise level in forward propagation.
The practical approach is to use e−2π|kz |∆z for the evanescent spectral components for both forward
and inverse extrapolation. Even more practical is to simply zero the evanescent spectral components
on inverse extrapolation. In all that follows, it is assumed that e2πikz ∆z has one of these two practical
extensions implemented for evanescent spectral components.
The wavefield extrapolation expression (equation (4.56)) is more simply written in the Fourier
domain to suppress the integration that performs the inverse Fourier transform:
Now consider a further extrapolation to estimate the wavefield at the bottom of layer 2 (z =
∆z1 + ∆z2 ). This can be written as
where T (kx , f ) is a correction factor for the transmission loss endured by the upgoing wave as it
crossed from layer 2 into layer 1. If transmission loss correction is to be incorporated, then the
data must not have had such amplitude corrections applied already. This is extremely unlikely
because seismic data is generally treated with a statistical amplitude balance that will compensate
for transmission losses. Also, correcting for transmission losses at each extrapolation step can be
numerically unstable because the correction factors are generally greater than unity. Consequently,
it is customary to set t(kx , f ) = 1. With this and incorporating equation (4.58), equation (4.59)
4.4. FOURIER METHODS 123
Upper Datum
L
θ
θ
∆z
Lower Datum
Figure 4.58: A geometric interpretation of the extrapolation phase shift. A dipping reflec-
tor emits a monochromatic signal that is received at two different datums. The extrapola-
tion phase shift is the phase difference between this signal as it arrives at a specific (x0 , z0 )
location on the upper datum and the signal as it arrives at the lower datum at the same
horizontal coordinate, i.e. at (x0 , z0 + ∆z).
becomes 2
φ(kx , z = ∆z1 + ∆z2 , f ) = φ0 (kx , f )e2πi m=1 kzm ∆zm
(4.60)
which has the immediate generalization to n layers
n n
φ(kx , z = m=1 ∆zm ,f )=φ0 (kx ,f )e
2πi m=1 kzm ∆zm . (4.61)
In the limit as ∆zm → dz, the summation in the extrapolator phase becomes an integral with the
result z
φ(kx , z, f ) = φ0 (kx , f )e2πi 0 kz (z )dz . (4.62)
2
where kz (z) = f 2 /v̂(z)2 − kx2 . This result is known as a first-order WKBJ solution and can be
derived as an approximate solution to the scalar wave equation for v̂(z) (Aki and Richards, 1980).
The second-order WKBJ solution arises when transmission losses are considered.
Equation (4.61) expresses what can be realized in a practical implementation of recursive wave-
field extrapolation while equation (4.62) is the theoretical limit of the process. The theoretically
correct phase shift from z1 to z2 is 2π times the area under the kz (z) curve between these depths.
It turns out that this result is true in one, z two or three dimensions
z but the form of kz (z) changes.
In one dimension, kz (z) = f /v̂(z) and 0 kz (z )dz = f 0 v̂(z )−1 dz = f τ (z) where τ (z) is the
vertical traveltime. In three dimensions, kz (z) = f 2 /v̂(z)2 − kx2 − ky2 and the integral of kz (z) has
no immediate simplification.
An interpretation of the extrapolation phase shift is suggested by incorporating sin θ = v̂kx /f
with the result
2
f2 2
∆z v̂kx
2πkz ∆z = 2π∆z − kx = 2πf 1− = 2πf τ cos θ (4.63)
v̂ 2 v̂ f
where τ = ∆z/v̂ has been used. Thus for constant velocity, the extrapolation phase shift is
2 The name WKBJ is the initials of Wentzel, Kramers, Brilloin, and Jeffreys who all derived it independently. The
Datum
∆L
θ
∆z
∆z
∆L
θ
∆z
∆z
∆z
D
ip
pi
ng
re
fle
ct
or
Figure 4.59: The recursive application of the constant-velocity phase shift can model v(z)
media as in equation (4.61).
which it the same result as equation (4.63). For a recursive sequence of many constant-velocity
phase shifts, the geometric interpretation is as shown in Figure 4.59.
The extrapolation operator is often separated into two terms, one that accomplishes a bulk time
shift and another that accomplishes dip-dependent focusing. The bulk time delay operator is simply
the extrapolation operator evaluated at kx = 0 (i.e. zero dip) and is called the static shift operator
or thin lens operator. The phase shift it applies is simply µs = 2πf ∆z/v̂. Denoting the total phase
shift by µ, this suggests that the focusing phase shift is simply µf = µ − µs . Thus,
" #
f2 2πf ∆z 2πf ∆z k x
2 v̂ 2
µf = µ − µs = 2π∆z − kx2 − = 1− 2 −1 (4.65)
v̂ 2 v̂ v̂ f
where
2πf ∆z
µs = (4.66)
v̂
4.4. FOURIER METHODS 125
.
seconds
.
.
.
- -
meters
Figure 4.60: A diffraction chart showing the response of five point diffractors as recorded
on the surface z = 0.
and, in summary,
total phase shift = µ = µs + µf . (4.67)
The focusing phase shift is angle-dependent and vanishes for kx = 0. That the static phase shift
accomplishes a bulk time delay can be seen as a consequence of the phase shift theorem of digital
signal theory (e.g. Karl (1989) page 87). This well-known result from time-series analysis says that
a time shift is accomplished in the Fourier domain by a phase shift where the phase is a linear
function of frequency. The slope of the linear phase function determines the magnitude of the time
shift. Rather than merely quote this result, it is instructive to actually demonstrate it. The static
phase shift can be applied to the ERM seismogram by
ψs (x, t) = φ0 (kx , f )e−µs +2πi(kx x−f t) dkx df. (4.68)
V∞
Since the static phase shift is independent of kx , the kx integral can be done directly to give
∞
ψs (x, t) = ψ̂0 (x, f )e−2πi(f t+f ∆z/v̂) df (4.69)
−∞
where ψ̂0 (x, f ) is the temporal Fourier transform of the ERM seismogram. Letting τ = ∆z/v̂, this
becomes ∞
ψs (x, t) = ψ̂0 (x, f )e−2πif (t+τ ) df = ψ0 (x, t + τ ). (4.70)
−∞
The last step follows because the integral is an inverse Fourier transform for the time coordinate
t + τ . This result is simply a proof of the phase shift theorem in the present context.
. .
seconds
seconds
. .
. .
A B
. .
- - - -
meters meters
Figure 4.61: The first hyperbola of diffraction chart of Figure 4.60 is shown convolved
with the its time reverse. (A) The focusing term only is applied. (B) Both the focusing
and the thin-lens term are applied
. .
seconds
seconds
. .
. .
A B
. .
- - - -
meters meters
Figure 4.62: The second hyperbola of diffraction chart of Figure 4.60 is shown convolved
with the time reverse of the first hyperbola. (A) The focusing term only is applied. (B)
Both the focusing and the thin-lens term are applied
Conceptually, the process just described is very simple though it might seem tedious to compute
for all samples. Could this cross correlation process be the space-time equivalent of extrapolation?
In fact, it is and to realize this it helps to recall from time series analysis that the cross correlation
of a(t) with b(t) can be done by time reversing b(t) and convolving. That is a(t) ⊗ b(t) = a(t) • b(−t).
This result carries over into two dimensions so that the two dimensional cross correlation of the ERM
seismogram with the first diffraction curve can be done by reversing the curve in both space and
time and then doing a two dimensional convolution. In this case the diffraction curve is symmetric in
space (it will not always be though) to the reversal in space does nothing. However, the time reversal
flips the first diffraction curve upside down to produce the concave operator that was envisioned.
Figures 4.61 and 4.62 illustrate these concepts with the extrapolation of the ERM seismogram of
Figure 4.60 to the depth of the first point diffractor. The extrapolation operator is the time reverse
of the first diffraction curve. In Figure 4.61A, only the focusing term is applied to the first diffraction
curve with the result that the curve is focused at its apex. In Figure 4.61B, the thin lens term is also
included so that the focused diffraction curve is shifted to time zero. In Figure 4.62A, the focusing
term is applied to the second diffraction curve and only partial focusing is achieved. When the thin
lens term is included in Figure 4.62B, it is apparent that the partially-focused second hyperbola is
now equivalent to the first hyperbola.
This perspective of wavefield extrapolation is very powerful and general. The cross correlation
argument shows that to extrapolate a wavefield to a particular depth, the response of a point
diffractor at that depth as viewed from the current datum must be constructed. Then the wavefield
is cross correlated with the diffraction response. The convolution argument is completely equivalent
and graphically illustrates the effects of the two parts of the operator: focusing and static shift. It
also shows how the extrapolation operator is similar to, but simpler, than the migration operator.
time coordinate for interpretation. Called τ , this migrated time may be created from z by doing a
simple one dimensional stretch using the migration velocity model as discussed in section 4.2.1. It
is always correct to do this following migration with v(x, z) with the stretch being defined by
z
dz
τ (x, z) =
(4.71)
0 v̂(x, z )
A time display, ψ(x, τ ), created from a migrated depth section, ψ(x, z), using equation (4.71) is
called a migrated time display. Generally this has a meaning distinct from that of time migration.
The two are equivalent only when lateral velocity variations are absent.
Time migration seeks to create ψ(x, τ ) directly without first creating ψ(x, z). To see how this
might be done, recall that the extrapolation phase shift can be written in the (kx , f ) domain as
where the static phase shift, µs is given by equation (4.66) and the focusing phase shift is given
by (4.65). If this extrapolator is applied in a recursive scheme, it is µs that progressively moves
data to earlier times so that each depth sample can be estimated through the imaging condition
ψ(x, z, t = 0). If an extrapolation were done with only µf then diffractions would focus at their
apex time but never move to tome zero and flat events (i.e. kx = 0) would never move at all. This
is exactly the behavior desired of a time migration. In fact, a time migration can be computed by
recursive phase shift using the extrapolation equation
where µf m is
" #
k 2 v̂ 2
µf = 2π∆τ 1 − x 2m − 1 (4.74)
f
where ∆τ = ∆z/v̂m has been used. When computing a time migration by recursive extrapolation,
the imaging condition must be changed because data not longer moves to t = 0 as it is focused.
Instead, focused events are found at t = τ so that the time migration is given by φ(x, τ, t = τ ).
Recursive extrapolation with equation (4.73), or some approximation to it, is a form of time
migration while recursive extrapolation with equation (4.72) is a depth migration. There is no
essential difference between the two for v(z) because the depth migration can be stretched to time
with equation (4.71) or the time migration can be stretched to depth with
τ
z(τ ) = v̂(x, τ )dτ (4.75)
0
To see this last point more clearly, if a recursive time migration is formulated as was done for depth
migration that lead to equation (4.62), it will be seen that the time migration omits the accumulated
phase shift operator z
z dz
ei 0 µs (z )dz = exp i2πf
. (4.76)
0 v(z )
A stretch of the time migration to depth effectively applies this operator. Time migration is thus
equivalent to depth migration for v(z) because the static delay, ∆z/v̂(z) = ∆τ , does not depend on
x. Since data moves all around during migration following various raypaths, any accumulated time
delays in ψ(x, τ, t = τ ) will be complicated functions of dip and position if v̂ varies with x. When v̂
4.4. FOURIER METHODS 129
x x
z vhigh vlow z vhigh vlow
∆z ∆z
A x
point diffractor
B x
t t
true image
point
tra
apparent image ve
lti
point m
ec
ur
ve
x
z vhigh vlow
∆z
C x
t
Figure 4.63: (A) A point diffractor sits beneath a low-velocity wedge (top) and its ERM
seismogram shows an apparent image point displaced to the left (bottom). (B) After
application of the focusing phase shift, the diffraction curve is partially focused but the
apparent image point has not moved. (C) The static phase shift moves the apparent image
point towards the true image point.
varies only with z, then all data at fixed τ have the same accumulated delay regardless of migration
raypath. When v̂ varies laterally, time migration actually refracts data in a way that systematically
violates Snell’s law as was discussed in section 4.2.3.
Figure 4.63A shows the ERM response of a point diffractor in a medium with v̂(x, z). The
apparent image point, which is the crest of the traveltime curve, is displaced to the left because
velocity is higher there. In Figure 4.63B, the ERM seismogram has been extrapolated to half the
depth of the diffractor using the focusing phase shift only. Since µf vanishes for kx = 0 it cannot
shift the crest of the traveltime curve even if it can apply the lateral velocity variations. All it
achieves is a partial focusing of the diffraction response. In Figure 4.63C, the static phase shift has
been applied and has had the effect of moving the crest of the traveltime curve to the right towards
the true image point. This happens because the time advance of the static shift is greater for lower
velocities than for higher ones. Thus it is apparent that extrapolation using µs and µf can focus
the diffractor at the proper location while µf alone can never do so.
130 CHAPTER 4. ELEMENTARY MIGRATION METHODS
An interpretation of this statement is that the integral of a function φ over the interval a ≤ x ≤ b is
found by evaluating a different function φ at the end points of the interval. The theorem assures us
that if φ exists then so does φ under a fairly general set of conditions. These functions are of course
mathematically related and φ is said to be the integral of φ or, equivalently, φ is the derivative of
φ.
Gauss’s theorem generalizes this basic result to dimensions greater than one. That is, it says that
if the integral over a volume of a certain function is desired, then it can be obtained by evaluating
another related function over the surface that bounds the volume. In higher dimensions, the notion
of direction is more complicated than the + or - needed in one dimension and vector functions
express this. Gauss’ theorem is usually written
$
# #
∇ · A dvol = A# · #n dsurf. (4.78)
V ∂V
Here A# is a vector function that might physically represent something like fluid flow or electric field
and #n is the outward pointing normal to the surface bounding the volume of integration (Figure 4.64).
4.5. KIRCHHOFF METHODS 131
n
volume
Equation (4.78) generalizes equation (4.77) in several ways. First, a vector function A# generalizes
the notion of direction where, in one dimension, the sign of φ was sufficient. Second, the derivative
has been generalized from φ to ∇# · A,
# and is called the divergence of A.
# Finally, the dot product of
#
A with #n integrated over the bounding surface generalizes the simple difference φ(b) − φ(a). In the
one dimensional case, the outward pointing normal points in the +x direction at b and in the −x
direction at a and the surface integral degenerates to a simple sum of two end members.
In many important cases, the vector function A# can be calculated as the gradient of a scalar
# = ∇φ.
potential A # In this case Gauss’ theorem becomes
$
2 ∂φ
∇ φ dvol = dsurf (4.79)
V ∂V ∂n
# · ∇φ
where ∇2 φ = ∇ # and ∂φ = ∇φ# · #n have been used. The meaning of ∂φ
= ∂n φ is that it is the
∂n ∂n
#
component of the vector ∇φ that is normal to the surface.
Now, return to equation 4.77 and consider the case when φ = φ1 φ2 . Then φ = φ2 φ1 + φ1 φ2 and
b
[φ2 φ1 + φ1 φ2 ] dx = φ1 φ2 ba (4.80)
a
or
b b
φ2 φ1 dx = φ1 φ2 ba − φ1 φ2 dx (4.81)
a a
which is the formula for integration by parts. An analogous formula in higher dimensions arises by
# 1 into equation (4.78). Using the identity ∇
substituting A = φ2 ∇φ # · a∇b
# = ∇a# · ∇b
# + a∇2 b leads
immediately to $
∂φ1
# 2 · ∇φ
∇φ # 1 + φ2 ∇2 φ1 dvol = φ2 dsurf. (4.82)
V ∂V ∂n
132 CHAPTER 4. ELEMENTARY MIGRATION METHODS
1 ∂ 2 ψ(#x, t)
∇2 ψ(#x, t) = (4.85)
v2 ∂t2
where the velocity v may depend upon position or it may be constant. To eliminate the time
dependence in equation (4.85), let ψ be given by a single Fourier component ψ(#x, t) = ψ̂(#x)e−2πif t .
Then equation (4.85) becomes the Helmholtz equation
where k02 = 4π 2 f 2 /v02 with v0 constant over all space and δ(#x − #x0 ) represents a source at #x = #x0 .
The analytic solution to equation (4.87) is well known (e.g. Morse and Feshbach (1953) page 810)
and can be build from a linear combination of the two functions
e±ik0 r
g ± (#x; #x0 ) = (4.88)
r
where r = |#x − #x0 | and three spatial dimensions are assumed. In two dimensions, the solution must
be expressed with Hankel functions that have the asymptotic form
± 2π ±ikr+iπ/4
g (#x; #x0 ) ∼ e , r → ∞. (4.89)
kr
4.5. KIRCHHOFF METHODS 133
Since g(#x; #x0 )e−2πif t is the time-dependent Green’s function, it is apparent that g + = r−1 eik0 r
corresponds to a wavefield traveling out from r = 0 while g − = r−1 e−ik0 r is a wavefield traveling
inward towards r = 0. In modelling, g + is commonly used and is called the causal Green’s function
while, in migration, it turns out that g − is appropriate and it is called the anticausal Green’s function.
Now, apply Greens theorem (equation 4.84) using ψ̂ and g − to get
$ " #
− 2 2 − − ∂ ψ̂(#x) ∂g − (#x; #x0 )
g (#x; #x0 )∇ ψ̂(#x) − ψ̂(#x)∇ g (#x; #x0 ) dvol = g (#x; #x0 ) − ψ̂(#x) dsurf.
V ∂V ∂n ∂n
(4.90)
Substituting equations (4.86) and (4.87) into the left hand side of this expression leads to
$ " #
2 ∂ ψ̂(#
x ) ∂g −
(#
x ; #
x 0 )
k − k02 g − (#x; #x0 )ψ̂(#x) dvol+ ψ̂(#x)δ(#x−#x0 ) dvol = g − (#x; #x0 ) − ψ̂(#x) dsurf
V V ∂V ∂n ∂n
(4.91)
Assuming that the point #x0 is interior to the volume V , the delta function collapses the second
integral on the left and this expression can be rewritten as
$ " −
#
∂ ψ̂(#
x ) ∂g (#
x ; #
x 0 )
ψ̂(#x0 ) = Λ(#x0 ) + g − (#x; #x0 ) − ψ̂(#x) dsurf (4.92)
∂V ∂n ∂n
where
2
Λ(#x0 ) ≡ k − k02 g − (#x; #x0 )ψ̂(#x) dvol. (4.93)
V
Equation (4.92) estimates the wavefield ψ̂ at the point #x0 interior to V as a volume integral plus
a surface integral over ∂V . The surface integral is what is desired since we can hope to know ψ̂
over the boundary of V . However, the volume integral involves the unknown ψ̂ and is essentially
not computable. The function Λ(#x0 ) expresses this volume integral and can be seen to vanish if the
reference medium v0 is equivalent to the actual medium v over the entire volume. Since g has been
chosen as a constant velocity Green’s function, Λ can only vanish precisely for constant velocity.
However, in the variable velocity case, approximate, ray theoretic Green’s functions can be used to
help minimize Λ (Docherty, 1991). To the extent that the reference medium does not equal the true
medium, then Λ expresses the error in a ψ̂ that is computed without Λ. In any case, the next step
is to drop Λ and substitute in g − = e−ik0 r /r into equation (4.92) with the result
$ " −ik0 r #
e ∂ ψ̂(#x) ∂ e−ik0 r
ψ̂(#x0 ) = − ψ̂(#x) dsurf. (4.94)
∂V r ∂n ∂n r
Multiplying both sides of this result by e−2πif t and recalling that ψ(#x0 , t) = ψ̂(#x0 )e−2πif t gives
$ " −ik0 r #
−2πif t e ∂ ψ̂(#x) ik0 ψ̂(#x)e−ik0 r ∂r e−ik0 r ∂r
ψ(#x0 , t) = e + + ψ̂(#x) 2 dsurf (4.96)
∂V r ∂n r ∂n r ∂n
134 CHAPTER 4. ELEMENTARY MIGRATION METHODS
where the time derivative in the second term results from ∂t ψ = −2πif ψ. This is a famous result and
is known as Kirchhoff’s diffraction integral. (In most textbooks this integral is derived for forward
modelling with the result that all of the terms are evaluated at the retarded time t − r/v0 instead of
the advanced time.) It expresses the wavefield at the observation point #x0 at time t in terms of the
wavefield on the boundary ∂V at the advanced time t + r/v0 . As with Fourier theory, it appears
that knowledge of both ψ and ∂n ψ are necessary to reconstruct the wavefield at an internal point.
source xs S x receiver
θ
S∞ r S∞
x
scatterpoint
Sz
image source
Figure 4.65: The geometry for Kirchhoff migration is shown. The integration surface is
S0 + Sz + S∞ and it is argued that only S0 contributes meaningfully to the estimation of
the backscattered field at vx0 .
S∞ , though it may contribute, can never be realized due to finite aperture limitations and its neglect
may introduce unavoidable artifacts. With these considerations, equation (4.98) becomes
$ " #
−1 ∂ψ 1 ∂r ∂ψ
ψ(#x0 , t) = + dsurf (4.99)
S0 r ∂z t+r/v0 v0 r ∂z ∂t t+r/v0
where the signs on the terms arise because #n is the outward normal and z is increasing downward
so that ∂n = −∂z .
Now, ∂z ψ must be evaluated. Figure 4.65 shows the source wavefield being scattered from the
reflector at #x0 which is called the scatterpoint. A simple model for ψ is that it is approximately
the wavefield from a point source, placed at the image source location, that passes through the
scatterpoint to the receiver. This can be expressed as
1 [A]t−r/v
ψ(#x, t) ∼ A(t − r/v) = (4.100)
r r
where A(t) is the source waveform at the scatterpoint. Using the chain rule gives
" #
∂ψ ∂r ∂ψ ∂r −1 ∂A [A]t−r/v
= = − . (4.101)
∂z ∂z ∂r ∂z vr ∂t t−r/v r2
When this is substituted into equation 4.99, the two terms in square brackets become similar both
involving the time derivative of the advanced wavefield. These will combine if v0 is now taken to be
the same as v. Thus $
2 ∂r ∂ψ
ψ(#x0 , t) = dsurf (4.103)
S0 vr ∂z ∂t t+r/v
136 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Finally, consider ∂z r. Since r = (x − x0 )2 + (y − y0 )2 + (z − z0 )2 this can be written
∂r ∂ z
= (x − x0 )2 + (y − y0 )2 + (z − z0 )2 = = cos θ (4.104)
∂z ∂z r
where θ is the vertical angle between the receiver location and ray to the scatterpoint. With this,
the final formula for the scattered wavefield just above the reflector is
$
2 cos θ ∂ψ
ψ(#x0 , t) = dsurf (4.105)
S0 vr ∂t t+r/v
This result, derived by many authors including Schneider (1978) and Scales (1995), expresses mi-
gration by summation along hyperbolic travelpaths through the input data space. The hyperbolic
summation is somewhat hidden by the notation but is indicated by [∂t ψ]2r/v . Recall that this no-
tation means that the expression in square brackets is to be evaluated at the time indicated by the
subscript. That is, as ∂t ψ(#x, t) is integrated over the z = 0 plane, only those specific traveltimes
values are selected that obey
2r 2 (x − x0 )2 + (y − y0 )2 + z02
t= = (4.107)
v v
which is the equation of a zero-offset diffraction hyperbola. When squared, this result is a three
dimensional version of equation (4.36).
In addition to diffraction summation, equation (4.106) requires that the data be scaled by
4 cos θ/(vr) and that the time derivative be taken before summation. These additional details
were not indicated by the simple geometric theory of section 4.2.4 and are major benefits of Kirch-
hoff theory. It is these sort of corrections that are necessary to move towards the goal of creating
bandlimited reflectivity. The same correction procedures are contained implicitly in f -k migration.
The considerations taken in deriving equation (4.106) suggest why ERM migration does not
achieve correct amplitudes. As mentioned following equation (4.105), a model linking the backscat-
tered wavefield to reflectivity was required. A more physical model will illustrate the shortcomings
of the exploding reflector model. Such a model has been advanced by Berkhout (1985) and others.
They consider that the source wavefield propagates directly to the reflector, undergoing transmission
losses and geometrical spreading but without multiples and converted modes. At the reflector it is
scattered upward with an angle-dependent reflection coefficient and then propagates upward to the
receivers, with further geometrical spreading and transmission losses but again without multiples
and mode conversions. This allows an interpretation of equation (4.105) as equivalent to the wave-
4.6. FINITE DIFFERENCE METHODS 137
field propagated from the source to the scatterpoint and scaled by the reflection coefficient. Thus, a
better way to estimate the reflectivity is to produce a model of the source wavefield as propagated
down to the scatterpoint and divide the result of equation (4.105) by this modeled wavefield. This is
a very general imaging condition called the deconvolution imaging condition that works for pre-stack
and zero offset data. The ERM imaging condition is kinematically correct but does not achieve the
same amplitudes as the deconvolution imaging condition.
Kirchhoff migration is one of the most adaptable migration schemes available. It can be eas-
ily modified to account for such difficulties as topography, irregular recording geometry, pre-stack
migration, and converted wave imaging. When formulated as a depth migration, it tends to be a
slow method because great care must be taken in the raytracing (Gray, 1986). When formulated as
a time migration, straight-ray calculations using RMS velocities can be done to greatly speed the
process. Another advantage of Kirchhoff methods is the ability to perform ”target-oriented” migra-
tions. That is, equation (4.106) need only be evaluated for those points in (x,z) which comprise the
target. Since the cost of computation is directly proportional to the number of output points, this
can greatly reduce the run times and makes migration parameter testing very feasible.
If finite difference extrapolation is equivalent to using a truncated Taylor series, then what
relation does wavefield extrapolation by phase shift have with Taylor series? Recall that if φ(z) is a
solution to the scalar wave equation in the Fourier domain, then it can be extrapolated by
As an illustration of the connection between Taylor series and the phase-shift extrapolator, consider
using Taylor series applied to equation (4.113) to derive something equivalent to equation (4.112).
From equation (4.113) the various derivatives can be estimated with the result that the nth derivative
is
dn φ(z)
= [2πikz ]n φ(z). (4.114)
dz n
Then, using this result, a Taylor series expression for the extrapolation of φ(z) to φ(z + ∆z) is
1 1
φ(z + ∆z) = φ(z) + [2πikz ]φ(z)∆z + [2πikz ]2 φ(z)∆z 2 + [2πikz ]3 φ(z)∆z 3 + . . . (4.115)
2 6
or
1 1
φ(z + ∆z) = φ(z) 1 + 2πikz ∆z + [2πikz ∆z]2 + [2πikz ∆z]3 + . . . . (4.116)
2 6
At this point, recall the expression for the series expansion of an exponential is ex = 1 + x + x2 /2 +
x3 /6 + . . . which affords the conclusion that the infinite series in square brackets
of equation (4.116)
may be summed to obtain 1 + 2πikz ∆z + 12 [2πikz ∆z]2 + 16 [2πikz ∆z]3 + . . . = e2πikz ∆z . Thus, if
infinitely many terms are retained in the series, equation (4.116) is equivalent to equation (4.112).
It may be concluded that wavefield extrapolation by phase shift is equivalent to extrapolation with
an infinite-order Taylor series and there is no upper limit on the allowed size of ∆z (in constant
velocity). Alternatively, wavefield extrapolation by finite difference approximations is equivalent to
extrapolation with a truncated Taylor series and the ∆z step size will have a definite upper limit of
validity.
while the backward difference is centered at z − ∆z/2. This suggests forming a centered difference
by averaging the forward and backward operators
1) Develop a suitable rational approximation to the one way dispersion relation for scalar waves.
A rational approximation is one which eliminates the square root though it may involve mul-
tiplication, division, and powers of frequencies and wavenumbers.
2) Construct a space-time domain differential equation from the approximate dispersion relation
by the replacement rules
i ∂ −i ∂ −i ∂ −i ∂
f→ ; kx → ; ky → ; kz → ; (4.121)
2π ∂t 2π ∂x 2π ∂y 2π ∂z
3) Choose an appropriate form for the finite difference operators (i.e. forward, backward,or central
differences).
4) Develop the difference equation which corresponds to the differential equation found in step 2.
As an illustration, the 15 degree finite difference time migration algorithm will be developed.
This was one of the first migration algorithms to be developed and is described in Claerbout (1976).
Recall that in the Fourier domain, time migration by recursive extrapolation proceeds by using the
140 CHAPTER 4. ELEMENTARY MIGRATION METHODS
0.4
C
0
kz A 15
o
- 0.4 B
- 0.8
-1 - 0.5 0 0.5 1
kx
Figure 4.66: For a particular f and v̂, (A) the dispersion relation (equation (4.125)) of
the parabolic wave equation; (B) the exact dispersion relation (equation (4.123)) of the
focusing phase shift; (C) the full dispersion relation of the scalar wave equation.
Since the general form of an extrapolation phase shift is phase = 2πkz ∆z, a finite difference time-
migration requires a “wave equation” whose dispersion relation approximates
" #
f kx2 v̂ 2
k̃z = 1− 2 −1 . (4.123)
v̂ f
where k̃z is used rather than kz to reserve the latter term for the real vertical wavenumber. It is not
acceptable to perform the substitutions of relation (4.121) into this result because the square root of
a differential operator results and that is difficult to work with. Therefore, a rational approximation
to the square root in equation (4.123) is desired. The second term in the square root is sin2 θ where
θ is the scattering angle. For small angles (i.e. energy traveling nearly vertically) it is expected that
an approximate form can be obtained by expanding the square root and keeping only the first few
terms
k 2 v̂ 2 k 2 v̂ 2 3k 2 v̂ 2
1 − x2 = 1 − x 2 + x2 + . . . . (4.124)
f 2f 8f
Truncating equation (4.124) at two terms and substituting the result into equation (4.123) gives
f k 2 v̂ 2 k 2 v̂ 2
k̃z ≈ 1− x 2 −1 =− x . (4.125)
v̂ 2f 2f
To construct a partial differential equation from this approximate dispersion, first multiply both
sides by 2f /v̂ and replace the spectral multiplications with partial derivatives according to relation
4.6. FINITE DIFFERENCE METHODS 141
2 ∂ 2 ψ(x, τ, t) ∂ 2 ψ(x, τ, t)
+ = 0. (4.127)
v̂ 2 ∂τ ∂t ∂x2
Equation (4.127) is an approximate wave equation which can be used to migrate stacked data from
zero offset time to migrated time. Though this equation is rarely used anymore, its finite difference
implementation illustrates most of the features of the method and is simpler than higher order
methods.
The finite difference approximation to equation (4.127) is well described by Claerbout (1976)
and that approach is followed here. Let
1 2 ∂2
T = δ x ≈ and ψjk = ψ(x, k∆τ, j∆t). (4.128)
∆x2 ∂x2
Then the x derivatives of equation (4.127) are approximated as
v̂ 2
δτ δt ψjk = − T ψjk (4.129)
2∆x2
where δτ and δt are, as yet, unspecified finite difference operators. The δt operator is implemented
as a forward difference but the right-hand-side of equation (4.129) is also modified to ensure that
both sides of the equation are centered at the same grid location. The result is
k ∆tv̂ 2 k
δτ ψj+1 − ψjk = − 2
T ψj+1 + ψjk (4.130)
4∆x
where the right-hand-side represents the average of T applied to two grid points. This process
of maintaining grid balance is essential in producing a stable algorithm and is a Crank-Nicholson
method. Next the τ operator is implemented in a similar way as a balanced forward difference with
the result
k+1 k+1 ∆τ ∆tv̂ 2 k+1
k
ψj+1 − ψj+1 − ψj − ψjk = − 2
k
T ψj+1 + ψj+1 + ψjk+1 + ψjk . (4.131)
8∆x
Finally, to downward continue, equation (4.131) must be solved for ψjk+1 . This can be written
as k+1
[I − aT ] ψjk+1 = [I + aT ] ψj+1 + ψjk − [I − aT ] ψj+1
k
(4.132)
where a = ∆τ ∆tv̂ 2 /8∆x2 . Equation (4.132) is solved numerically for the unknown ψjk+1 where it
is assumed that all of the quantities on the right-hand-side are known. This is an example of an
implicit finite difference method because it requires the numerical inversion of the matrix operator
I + aT .
142 CHAPTER 4. ELEMENTARY MIGRATION METHODS
j
ψ 11 ψ 12 ψ 13 ___ ψ n-1 1
ψ 22 ψ n-2 1
k ψ 33
Figure 4.67: The parabolic wave equation is solved on a 2D grid using finite differences.
The first row contains the known data and the last column is set to zero. The solution
proceeds row-by-row and begins in the upper right corner.
The differencing procedure can be viewed on a (j, k) grid as shown in Figure 4.67. The x axis
is orthogonal to this figure and is not shown. The first row contains the measured CMP stack and
the last column (corresponding to maximum time) is set to zero. Consider the four grid points in
1
the upper right corner. Three of these are prescribed by initial conditions and the fourth, ψn−1 can
1
then be calculated by equation (4.132). Subsequently, ψn−2 can be solved for and so on until the
entire second row in calculated. Computation then moves to the third row and continues until the
entire grid is completed. Each row corresponds to an extrapolated τ section and the final migrated
section corresponds to t = τ which is the diagonal.
From this example, many of the well known characteristics of finite difference migration algo-
rithms are easily deduced such as
• They use an approximate form of the exact dispersion relation for scalar waves. This means
that they are all dip or scattering angle limited.
• They are very flexible in handling velocity variations since the velocity may simply be varied
with each grid cell.
• The use of finite difference derivatives means that even the approximate dispersion relation is
not exactly realized.
• Finite difference derivatives are not unique and the accuracy of a method depends strongly on
the sophistication of the approximations used. Generally, finite difference methods need many
more samples per wavelength that Fourier methods (6-10 versus 2-3).
0
A B
frequencyW(Hz)
.5 40
seconds
1.0 80
120
1.0 2.0 -.06 0 .06
kilometers wavenumberW(m-1)
Figure 4.68: (A) An ERM seismogram that models a regular grid of point diffractors. (B)
The (kx , f ) spectrum of (A).
0 0
A B
frequencyW(Hz)
.5 40
seconds
1.0 80
120
1.5 3.0 -.06 0 .06
0 1.0 2.0
kilometers wavenumberW(m-1)
Figure 4.69: (A) The f -k migration of the ERM seismogram of Figure 4.68A. (B) The
(kx , f ) spectrum of (A).
0
.1 40
A
sec
Hz
B 80
.2
120 0
1.2 1.4 0 -.05 0 .05
1.0 km 1.0 m-1
100Wms 40 40
Hz
sec
sec
Hz
100Wm 80
1.1 1.1 80
120 120
.25 .4 2.6 2.8 -.05 0 .05 -.05 0 .05
km km 0 m-1
m-1
1.3 40
Hz
sec
80
1.4
120
1.2 1.4 -.05 0 .05
km m-1
Figure 4.70: (A) Enlargements of the boxed focal points of Figure 4.69A. (B) Local (kx , f )
spectra of the focal points in (A).
A
10 2.5 7.5 5 5 7.5 2.5 10
∆x km km km km km km km km
Z θ
6Wkm
vT o o
60o 23o 51 40o 40o 51o 23 60o
Figure 4.71: A ray from a scatterpoint on a re- Figure 4.72: A fixed length seismic line induces
flector requires a certain aperture A and record an aperture effect that varies laterally along the
length T to be captured. Also, the CMP sam- line. In turn, this limits the scattering angle
ple rate ∆x must be sufficiently small to avoid spectrum.
spatial aliasing.
146 CHAPTER 4. ELEMENTARY MIGRATION METHODS
the (kx , kz ) plane as was discussed in section 4.4.1 and illustrated in Figure 4.49. As can be deduced
from the figure, the kx bandwidth after migration is limited by
f2
kz2 = − kx2 (4.134)
(v/2)2
λdom
λmin = λdom = 2.85λmin = = 1.425λmin (4.135)
2
α
δx = (4.136)
kxmax
αv α
δx = = λmin (4.137)
2fmax sin (θmax ) 2 sin (θmax )
1.4 1.4v vT
δx ∼ λmin = = (4.138)
sin (θmax ) fmax sin (θmax ) 2
In this expression, kxmax is the limiting wavenumber to be expected from a migration with no
limitation on scattering angle. In any practical setting, there is always a limit on scattering angle
and it is generally spatially variant. The limit may be a result of the effects of finite aperture,
finite record length, and discrete spatial sample size or any of many other possibilities including:
the migration algorithm may be ”dip limited”, lateral or complex vertical velocity variations can
create shadow zones, attenuation effects are dependent on raypath length and hence affect the larger
scattering angles more. Whatever, the cause, a limitation of the range of scattering angles which can
be collected and focused to a particular point translates directly into a resolution limit as expressed
by equation (1). The size of the smallest resolvable feature, say δx, is inversely proportional to
kxlim . For definiteness, let kxlim = α/(2δx), where α is a proportionality constant near unity, and
solve for δx to get
αv
δx = . (4.139)
4fmax sin (θmax )
Since the aperture limit is x and z variant (and record length and spatial aliasing limits are z variant)
δx must also vary with position. An interpretation of equation (4.139) is that it gives the smallest
discernible feature on a reflector whose dip is normal to the bisector of the scattering angle cone at
any position (Figure 4.71).
For constant velocity, the limits imposed on zero-offset scattering angle are easily derived using
straight ray theory and have been shown (Lynn and Deregowski, 1981) to be
A
tan θA = (4.140)
z
and
vT
cos θT = . (4.141)
2z
In these expressions, A is the available aperture, T is the record length, z is depth, and v is the
presumed constant velocity. θA and θT are limitations on scattering angle imposed by aperture and
record length respectively (Figure 4.71). Aperture is defined as the horizontal distance from an
analysis point to the end of the seismic line or the edge of a 3-D patch and is thus dependent on
4.7. PRACTICAL CONSIDERATIONS OF FINITE DATASETS 147
120
90
angleW(degrees)
Rec
ord
angleW(degrees)
Wlen 80
60 gth Apertu
SpatialWaliasing re
Aperture 40 asing
30 SpatialWali
RecordWlength
0 0 3.0
1.0 3.0 5.0 1.0 2.0
DepthW(km) depthW(km)
Figure 4.73: Constant velocity scattering an- Figure 4.74: Constant gradient scattering an-
gle chart showing the aperture, record length, gle chart showing the aperture, record length,
and spatial aliasing limits for A = 2500m, T = and spatial aliasing limits for A = 2500m, T =
3.0sec, v = 3500m/s, ∆x = 20m and f = 60Hz. 3.0sec, v = 1500 + .6zm/s, ∆x = 20m and
f = 60Hz.
azimuth and position (Figure 4.72). Alternatively, the record length limit has no lateral variation.
Taken together, these expressions limit the scattering angle spectrum to a recordable subset.
A third limiting factor is spatial aliasing which further constrains the possible scattering angle
spectrum to that which can be properly imaged (migrated). (Liner and Gobeli (1996) and Liner
and Gobeli (1997) give an analysis of spatial aliasing in this context.) The Nyquist requirement is
that there must be at least two samples per horizontal wavelength to avoid aliasing
λ
λx = ≥ 2∆x. (4.142)
sinθx
Here, ∆x is the spatial sample size (CMP interval), λ and λx are wavelength and its apparent
horizontal component, and θx is most properly interpreted as the emergence angle of a dipping event
on a zero offset section. Consistent with zero offset migration theory, the exploding reflector model
Lowenthal and Sherwood (1976) can be used to relate wavelength to velocity through λ = v/(2f )
where f is some frequency of interest. This leads to an angle limited by spatial aliasing given by
v
sin θx = (4.143)
4f ∆x
In the constant velocity case, the emergence angle of a ray and the scattering angle at depth are
equal and thus equation (4.143) expresses the constant velocity limit on scattering angle imposed
by spatial aliasing. For vertical velocity variation, v(z), the result still applies provided that θx is
simply interpreted as emergence angle and v as near surface velocity. The emergence angle can be
related to the scattering angle at depth using Snell’s law. This is done by recalling that the ray
parameter, p = sin(θ(z))/v(z), is conserved (Slotnick, 1959), which leads to
v(z)
sin θx = . (4.144)
4f ∆x
This expression generalizes spatial aliasing considerations to monotonically increasing but otherwise
arbitrary v(z) and θx is interpreted as the scattering angle from depth, z.
148 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Returning to constant velocity, equations (4.140), (4.141), and (4.143) can be used to create a
scattering angle resolution chart for an assumed recording geometry, position on the line, frequency
of interest and constant velocity. A typical case is shown in Figure 4.73 where it is seen that the
aperture limit is concave upward and tends asymptotically to zero at infinite depth. The record
length limit has the opposite curvature and reaches zero degrees at the depth z = vT /2. Both limits
admit the possibility of 90◦ only for z = 0. The spatial aliasing limit is depth independent but
requires a frequency of interest which can conservatively be taken as the maximum (not dominant)
signal frequency.
Charts such as Figure 4.73 can be used as an aid in survey design but tend to give unrealistic
parameter estimates due to the assumption of straight raypaths. In most exploration settings,
velocity increases systematically with depth and thus raypaths bend upward as they propagate from
the scatterpoint to the surface. Intuitively, this should lead to shorter aperture requirements and
allow the possibility of recording scattering angles beyond 90◦ in the near surface. The spatial
aliasing limit has already been discussed in this context and the aperture and record length limits
will now be derived exactly for the case of a constant velocity gradient, that is when v(z) = vo + cz.
The derivation requires solution of the Snell’s law raypath integrals for the linear gradient case
(Slotnick, 1959). If pA is the ray parameter required to trace a ray from a scatterpoint to the end
of the spatial aperture, then z
p v(z )
A(z) = A dz . (4.145)
0 1 − p2A v(z )2
Similarly, let pT be the ray parameter for that ray from scatterpoint to the surface which has
traveltime (two-way) equal to the seismic record length, then
z
1
T (z) =
2
dz . (4.146)
2
0 v(z ) 1 − pT v(z )
and
γ −1 − cosh(cT /2)
cos (θT ) = , (4.150)
sinh(cT /2)
where
v0
γ= . (4.151)
v(z)
4.7. PRACTICAL CONSIDERATIONS OF FINITE DATASETS 149
When equations (4.144), (4.149), and (4.150) are used to create a scattering angle resolution
chart, the result is typified by Figure 4.74. The parameters chosen are the same as for Figure 4.73
and the linear velocity function was designed such that it reaches 3500 m/s (the value used in Figure
4.73) in the middle of the depth range of Figure 4.74. It can be seen that the possibility of recording
angles beyond 90◦ is predicted for the first 1000 m and the aperture limit is everywhere more broad
than in Figure 4.73. The record length limit forces the scattering angle spectrum to zero at about
3700 m compared to over 5000 m in the constant velocity case. This more severe limit is not always
the case, in fact a record length of 6 seconds will penetrate to over 12000 m in the linear velocity
case and only 10500 m in the constant case. Also apparent is the fact that the spatial aliasing limit
predicts quite severe aliasing in the shallow section though it gives exactly the same result at the
depth where v(z) = 3500 m/s.
4.7.2 Examples
Figures 4.75A, 4.75B, and 4.75C provide further comparisons between the linear v(z) theory and
constant velocity results. In Figure 4.75A, the aperture limits are contrasted for the same linear
velocity function (v = 1500 +.6 z m/s) and constant velocity (v 3500 m/s) used before. The dark
curves show the linear velocity results and the light curves emanating from 90◦ at zero depth are
the constant velocity results. Each set of curves covers the range of aperture values (from top to
bottom): 1000, 4000, 12000, and 20000 meters. The dramatic effect of the v(z) theory is especially
obvious for larger apertures which admit angles beyond 90◦ for a considerable range of depths.
Figure 4.75B is similar to Figure 4.75A except that the record length limit is explored. For each set
of curves, the record lengths shown are (from top to bottom): 2.0, 4.0, 6.0, and 8.0 seconds. Figure
4.75C shows spatial aliasing limits for a frequency of 60 Hz. and a range of ∆x values (from top to
bottom): 10 20 40 60 80 100 and 150 meters.
Next consider Figure 4.76 that shows a synthetic demonstration of the aperture effect. Here,
a number of unaliased point diffractor responses have been arranged at constant time. Thus the
record length limit is constant and the aliasing limit does not apply. When migrated, the resulting
display clearly shows the effect of finite spatial aperture on resolution. Comparison with Figure 4.72
shows the direct link between recorded scattering angle spectrum and resolution. Approximately,
the focal points appear as dipping reflector segments oriented such that the normal (to the segment)
bisects the captured scattering angle spectrum.
Figure 4.77 is a study designed to isolate the effects of temporal record length on resolution. The
unmigrated section is constructed such that all five point diffractors are limited by record length and
not by any other effect. Upon migration, the focal points are all symmetric about a vertical axis but
show systematic loss of lateral resolution with increasing time. As shown in Figures 4.73, 4.74, and
4.75B, the record length limit always forces the scattering angle spectrum to zero at the bottom of
the seismic section. Equation (4.139) then results in a lateral resolution size that approaches infinity.
This effect accounts for the often seen data ’smearing’ at the very bottom of seismic sections.
Figure 4.78 shows the migration of a single diffraction hyperbola with three different spatial
sample intervals to illustrate the resolution degradation that accompanies spatial aliasing. In Figure
4.78B, the migration was performed with a spatial sample rate of 4.5 m which represents a comfort-
ably unaliased situation. In Figures 4.78C and 4.78D the sample intervals are 9 m (slightly aliased)
and 18 m (badly aliased). The slightly aliased situation has not overly compromised resolution but
the badly aliased image is highly degraded. Note that the vertical size of the image is unaffected.
In Figure 4.79, the effect of maximum temporal frequency is examined. A single diffraction
hyperbola was migrated with three different maximum frequency limits . In Figure 4.79B, the focal
point and its local (kx , f ) spectrum are shown for an 80 Hz. maximum frequency while Figures 4.79C
150 CHAPTER 4. ELEMENTARY MIGRATION METHODS
180
140
A B
angle (degrees)
angle (degrees)
100
120 8s
4s
ec
ec
20 km
6
60
se
60
c
12 km
2 sec
4 km 20
0 1 km
0 5 10 15 20 0 5 10 15 20 25
depth (km) depth (km)
20 m 40 m
90
C
angle (degrees)
60
computed for 60 Hz
80 m
30
160 m
0 5 10 15 20
depth (km)
Figure 4.75: The constant and linear v(z) theories are compared. In each case, the linear
v(z) theory is presented with bold lines while the constant velocity theory uses thin lines.
The constant velocity used was 3500 m/s and the linear function was v = 1500 + .6z
m/s. (A) The aperture limits are contrasted. The apertures used are labeled on the bold
curves but the light curves use sequentially the same A values. (B) The record limits are
contrasted. The record lengths used are labeled on the bold curves but the light curves
use sequentially the same T values. (C) The spatial aliasing limits are contrasted. The
∆x values used are labeled on the bold curves but the light curves use sequentially the
same values.
4.7. PRACTICAL CONSIDERATIONS OF FINITE DATASETS 151
.8
seconds
1.0
1.2
1.4
A
.5 km 1.5 2.5
.8
seconds
1.0
1.2
500Wm
1.4 B
.5 km 1.5 2.5
Figure 4.76: (A) A series of point scatterers all at the same time so that their record
length effect is constant. (B) The migration of (A) shows the spatial variability of focal
point geometry caused by the aperture effect. Compare with Figure 4.72
0 0
A B 200Wm
.5 .5
seconds
1.50 1.5
1.0 2.0 3.0 1.1 1.5 1.9
km km
Figure 4.77: An illustration of the record length effect. (A) A series of diffractions that
are all truncated by the bottom of the section to minimize the aperture effect. (B) The
migration of (A) shows a systematic loss of lateral resolution with increasing time. At the
bottom of the section, the lateral size of a focal point is theoretically infinite. Note the
horizontal scale change between (A) and (B).
152 CHAPTER 4. ELEMENTARY MIGRATION METHODS
Hz
B
sec
.56 80
.5 4.5m
seconds
.48 20
Hz
sec
C
1.0 .56 9m
80
.48 20
Hz
1.5 D
sec
0 1.0 km 2.0 3.0 80
.56 18m
-.05 .05
100Wm m-1
Figure 4.78: The spatial aliasing effect is demonstrated. (A) A single diffraction hyperbola.
(B) A zoom of the focal point and the (kx , f ) spectrum after migration of (A) with ∆x =
4.5 m. (C) A zoom of the focal point and the (kx , f ) spectrum after migration of (A) with
∆x = 9 m. (D) A zoom of the focal point and the (kx , f ) spectrum after migration of (A)
with ∆x = 18 m.
and 4.79D are similar except that the maximum frequencies were 60 Hz. and 40 Hz. respectively. It is
clear from these figures that limiting temporal frequency affects both vertical and lateral resolution.
As expected from equation (4.139), a reduction of fmax from 80 Hz to 40 Hz causes a doubling of
the focal point size.
In summary, the theory of (kx , f ) migration predicts a simple model for the resolving power
of seismic data. The result is a spatial bandwidth that depends directly on frequency and sine of
scattering angle and inversely on velocity. Finite recording parameters (aperture and record length)
place space and time variant limits on the observable scattering angle spectrum. Thus the resolution
of a seismic line is a function of position within the aperture of the line. The scattering angle limits
imposed by aperture, record length, and spatial sampling can be derived exactly for the case of
constant velocity and for velocity linear with depth. The linear velocity results are more realistic
and lead to considerably different, and usually cheaper, survey parameters than the constant velocity
formulae.
4.7. PRACTICAL CONSIDERATIONS OF FINITE DATASETS 153
0 0
frequency (Hz)
A
.5 40
seconds
1.0 80
1.5 120
0 1.0 km 2.0 3.0 -.06 0 .06
m-1
focalWpoints f-kWspectra
.48 20
Hz
sec
fmax ˜ 80 Hz B
.56 80
.48 20
fmax ˜ 60 Hz C Hz
sec
80
.56
20
.48
Hz
fmax ˜ 40 Hz D
sec
80
.56
-.05 .05
100Wm m-1
Figure 4.79: The resolution effect of decreasing Fmax is demonstrated. (A) A single
diffraction hyperbola and its (kx , f ) spectrum. (B) A zoom of the focal point and the
(kx , f ) spectrum after migration of (A) with fmax = 80 Hz. (C) A zoom of the focal point
and the (kx , f ) spectrum after migration of (A) with fmax = 60 Hz. (D) A zoom of the
focal point and the (kx , f ) spectrum after migration of (A) with fmax = 40 Hz.
Bibliography
Aki, K., and Richards, P. G., 1980, Quantitative Seismology: W.H. Freeman.
Ames, W. F., 1992, Numerical Methods for Partial Differential Equations: Academic Press.
Berkhout, A. J., 1985, Seismic Migration: Developments in solid earth geophysics:, volume 14
Elsevier.
Chun, J. H., and Jacewitz, C. A., 1981, Fundamentals of frequency domain migration: Geophysics,
46, 717–733.
Clayton, R., and Engquist, B., 1977, Absorbing boundary conditions for acoustic and elastic wave
equations: Bull. Seis. Soc. Am., 67, 1529–1540.
Dix, C. H., 1955, Seismic velocities from surface measurements: Geophysics, 20, 68–86.
Docherty, P., 1991, A brief comparison of some Kirchhoff integral formulas for migration and inver-
sion: Geophysics, 56, 1164–1169.
Durran, D. R., 1999, Numerical Methods for Wave Equations in Geophysical Fluid Dynamics:
Springer.
Etter, D. E., 1996, Introduction to Matlab for Engineers and Scientists: Prentice-Hall.
Futterman, W. I., 1962, Dispersive body waves: J. Geophys. Res., 73, 3917–3935.
Gazdag, J., 1978, Wave-equation migration with the phase-shift method: Geophysics, 43, 1342–1351.
Gray, S. H., 1986, Efficient traveltime calculations for Kirchhoff migration: Geophysics, 51, 1685–
1688.
Gray, S. H., 1997, Trueamplitude seismic migration: A comparison of three approaches: Geophysics,
62, 929–936.
Gustafson, K. E., 1987, Partial Differential Equations and Hilbert Space Methods: John Wiley and
Sons.
154
BIBLIOGRAPHY 155
Hanselman, D. C., and Littlefield, B., 1998, Mastering MATLAB 5: A Comprehensive Tutorial and
Reference: Prentice Hall.
Higham, D. J., and Higham, N. J., 2000, MATLAB Guide: Society for Industrial and Applied
Mathematics.
Karl, J. H., 1989, An Introduction to Digital Signal Processing: Academic Press, Inc.
Kjartansson, E., 1980, Constant Q-wave propagation and theory: J. Geophys. Res., 84, 4737–4748.
Liner, C. T., and Gobeli, R., 1996, Bin size and linear v(z): Expanded abstracts 66th Annual Mtg.
Soc. Expl. Geophy.
Liner, C. T., and Gobeli, R., 1997, 3-d seismic survey design and linear v(z): Expanded abstracts
67th Annual Mtg. Soc. Expl. Geophy.
Lines, L. R., Slawinski, R., and Bording, R. P., 1999, A recipe for stability of finite-difference
wave-equation computations: Geophysics, 64, 967–969.
LowenthalD., L. L. R. R., and Sherwood, J. W., 1976, The wave equation applied to migration:
Geophysical Prospecting, 24, 380–399.
Lynn, H. B., and Deregowski, S., 1981, Dip limitations on migrated sections as a function of line
length and recording time: Geophysics, 46, 1392–1397.
Mitchell, A. R., and Griffiths, D. F., 1980, The Finite Difference Method in Partial Differential
Equations: John Wiley and Sons.
Morse, P. M., and Feshbach, H., 1953, Methods of Theoretical Physics: McGraw-Hill.
Press, W., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P., 1992, Numerical Recipes in C:
The Art of Scientific Computing: Cambridge University Press.
Redfern, D., and Campbell, C., 1998, The Matlab 5 Handbook: Springer-Verlag.
Robinson, E. A., 1983, Migration of Geophysical Data: International Human Resources Development
Corp. (IHRDC).
Schneider, W. S., 1978, Integral formulation for migration in two or three dimensions: Geophysics,
43, 49–76.
Schuster, G. T., 1997, Green’s function for migration: Expanded Abstracts 67th Annual Mtg. Soc.
Expl. Geophy., pages 1754–1757.
Taner, M. T., and Koehler, F., 1969, Velocity spectra-digital computer derivation and applications
of velocity functions: Geophysics, 34, 859–881.
Taner, M. T., Koehler, F., and Sheriff, R. E., 1979, Complex seismic trace analysis: Geophysics, 44,
1041–63.
156 BIBLIOGRAPHY
Wiggins, J. W., 1984, Kirchhoff integral extrapolation an migration of nonplanar data: Geophysics,
49, 1239–1248.
Zauderer, E., 1989, Partial Differential Equations of Applied Mathematics: John Wiley and Sons.