Method and apparatus for enhancing seismic data

A method and apparatus for stacking a plurality of seismic midpoint gathers to provide an enhanced pictorial representation of seismic events is disclosed. The approximate propagation velocity corresponding to a selected event in a common midpoint gather, is determined by summing the common midpoint gather using first and second weights to provide respective first and second weighted sums over offset based on an estimated velocity corresponding to the event, and developing from the sums a velocity error value indicative of the approximate error between the estimated velocity and the actual velocity. The common midpoint gather is then re-stacked in accordance with the determined propagation velocity to provide an enhanced pictorial representation of the seismic event. The first and second weighted sums are taken over a time window centered upon an estimated zero offset travel time for the event. The first and second weights can be selected to provide rapid, slow or intermediate convergence upon the true velocity. The velocity error value is determined as a function of the deviation of the peak of the first weighted sum from the center of the time window, relative to the deviation of the peak of the second weighted sum from the center of the time window. Alternatively, the velocty error value is determined as a function of the deviation of the peak of the crosscorrelation of the first and second weighted sums from the center of the time window.

FIELD OF THE INVENTION 
The present invention relates to the field of seismic prospecting, and in 
particular, the field of seismic signal processing for producing a 
pictorial representation of subsurface formations. 
BACKGROUND OF THE INVENTION 
One well-known technique of seismic prospecting involves the placement of a 
linear array of geophones along the surface of the earth, or a body of 
water, producing a "shot", or source of seismic signals, through the use 
of an explosive charge or vibratory stimulus, and receiving signals 
reflected from the subsurface formations, or seismic "events", at each of 
the geophones in the linear array. The locations of the geophones and the 
signal source are then moved, in a well known manner, and the process is 
repeated. In this way, raw data in the form of a plurality of "common shot 
gathers" are collected along a seismic line of interest. 
To analyze the raw data the data is gain-adjusted in a well-known manner to 
remove the influence of distance between various source and receiver 
pairs, and deconvolved to produce a relatively narrow pulse in response to 
each seismic event. The geometry of the common shot gathers is also 
converted in a well-known way to a plurality of "common midpoint gathers", 
one of which is illustrated in FIG. 1. The common midpoint gather is 
defined as those source and receiver pairs which have a common midpoint. 
In FIG. 1, S.sub.i denotes the source location, R.sub.i denotes the 
associated receiver location, x is the source-receiver offset, equal to 
R.sub.i -S.sub.i, and the common midpoint value is defined as (R.sub.i 
+S.sub.i)/2. The formula relating travel time of the seismic signal from 
the source S.sub.i (x) to the receiver R.sub.i (x) is given by v.sup.2 
t.sup.2 (x)=x.sup.2 +4d.sup.2, where v is velocity of propagation, t is 
travel time, and d is the distance from zero offset (x=0) to the seismic 
reflection. By defining 2d=vt(0), t(0) being the travel time of a signal 
from source S.sub. 0 to colocated receiver R.sub.0, the following formula 
results: 
EQU v.sup.2 t.sup.2 (x)=x.sup.2 +v.sup.2 t.sup.2 (0). 
This formula, which defines the normal moveout relationship, is the most 
commonly used method of determining the signal arrival time differences of 
seismic data as a function of the offset distance of receiver from the 
source. The normal moveout relationship is hyperbolic between offset x and 
time t(x). 
It is well known in the art that common midpoint data are noisy to the 
point where t(x) cannot be measured directly, and in order to reveal the 
unknown velocity v, corresponding to a particular event, a velocity 
analysis must be undertaken. Furthermore, since the velocity varies with 
both depth d, and position x, many such velocity analyses must be 
performed along the seismic line of interest to establish these 
variations. Clearly, the amount of time required to perform each analysis 
will greatly impact the time required to analyze an entire line. 
Heretofore, each velocity analysis has been performed in accordance with a 
technique such as that proposed by Tanner and Koehler, in "Velocity 
Spectra-Digital Computer Derivation and Applications of Velocity 
Functions", Geophysics, v. 34, pp. 859-899 (1969), as follows: A set of 
trial velocities are used to "stack" the data, and then either a manual or 
automatic search is made to determine that velocity which gives the best 
response to a particular event. This is explained with reference to FIGS. 
2-4. FIG. 2 illustrates common midpoint seismic data, which for the 
purposes of simplicity of explanation, is a noise-free single event, shown 
as a function of offset x and time t. In FIG. 2, the seismic "wavelet" 
corresponding to the event is shown as having a normal hyperbolic moveout 
relationship as defined above. The prior art method of velocity analysis 
is to "correct" the data with a suite of velocities, "stack" the data for 
each such velocity, and examine the result of the stacked data, choosing 
the velocity that results in the stack of highest amplitude. The 
assumption is that the highest amplitude is produced by the closest 
velocity. 
With reference to FIGS. 3A and 3B, a first trail velocity is chosen and the 
data is corrected, i.e. the effect of the normal moveout relationship is 
removed, making each recording at the various shotreceiver offsets look 
like they were recorded at x=0. This is done by shifting the collection of 
seismic wavelets up in time using the normal moveout relationship and the 
estimated trial velocity, as graphically illustrated in FIG. 3A. The 
shifted wavelets are then summed along offset x to produce summed trace or 
stack 2, FIG. 3B. 
FIG. 4 shows the shifted common midpoint data using a more correct estimate 
of velocity, resulting in a stack 4 having a higher amplitude. 
Having determined the velocity, for a particular event, within a reasonable 
tolerance, the stack of the common midpoint gather, or "zero offset trace" 
is used to portray the associated seismic event in pictorial form, an 
example of which will be referred to below. The stacked data is so used 
since the individual reflection signals are usually too noisy to 
accurately portray useful information. By stacking the common midpoint 
gathers using a good velocity estimate, the step of summing greatly 
enhances the information, by increasing the signal-to-noise ratio. The 
closer the velocity used to stack the data is to the actual velocity, the 
better the enhancement of the information will be. 
It will be appreciated that the traditional processing technique described 
above requires that each set of common midpoint data be stacked using a 
relatively large number of trial velocities depending upon the required 
tolerance in the velocity estimate, and any a priori knowledge of the 
propagation velocity at nearby or similar formations. Moreover, it will be 
appreciated that the stack is not performed using the exact (within 
tolerance) velocity required, but is rather based on an interpolation 
using the suite of trial velocities. 
OBJECTS AND SUMMARY OF THE INVENTION 
It is therefore an object of the present invention to overcome the 
shortcomings of prior art seismic processing techniques. 
It is a further object of the present invention to provide a new seismic 
processing technique which is more rapid and accurate than prior art 
techniques. 
It is a further object of the present invention to provide a new seismic 
processing technique to enhance the pictorial representation of common 
midpoint stack data. 
It is a further object of the present invention to provide a seismic 
processing technique for enhancing the pictorial representation of seismic 
data in which seismic events are more accurately portrayed. 
It is a further object of the present invention to provide an improved 
technique for pictorially representing seismic events at locations 
substantially corresponding to their actual locations. 
In accordance with a first aspect of the present invention, a method for 
stacking a plurality of seismic midpoint gathers, to provide an enhanced 
pictorial representation of seismic events, includes the steps of 
providing a common midpoint gather, determining the approximate 
propagation velocity corresponding to a selected event, based on the 
common midpoint gather, by summing the common midpoint gather using first 
and second weights to provide respective first and second weighted sums 
over offset, based on an estimated velocity corresponding to the event, 
and developing from the sums a velocity error value indicative of the 
approximate error between the estimated velocity and the actual velocity. 
The common midpoint gather is then re-stacked in accordance with the 
determined propagation velocity to thereby provide an enhanced pictorial 
representation of the seismic event. More specifically, the first and 
second weighted sums are taken over a time window centered upon an 
estimated zero offset travel time for the event. The first and second 
weights can be selected to provide rapid, slow or intermediate convergence 
upon the true velocity. The velocity error value is determined as a 
function of the deviation of the peak of the first weighted sum from the 
center of the time window, relative to the deviation of the peak of the 
second weighted sum from the center of the time window. Alternatively, and 
in accordance with the preferred embodiment, the velocity error value is 
determined as a function of the deviation of the peak of the 
crosscorrelation of the first and second weighted sums from the center of 
the time window. 
Additionally, an update zero offset travel time can be determined by 
developing from the first and second sums a time error value indicative of 
the approximate error between the estimated zero offset the travel time 
and the actual zero offset travel time. 
Furthermore, the propagation velocity, as determined above, (also referred 
to as "apparent RMS velocity" or "stacking velocity" versus time) can be 
converted to "true interval velocity" versus depth, or "true average 
velocity" versus time. The re-stacked common midpoint gather may then be 
converted, in accordance with the true interval velocity, or true average 
velocity, and zero offset time, into a pictorial representation of the 
seismic event at a location substantially corresponding to the actual 
spacial location of the event. 
In accordance with another aspect of the present invention, an apparatus 
for stacking a plurality of seismic midpoint gathers to provide an 
enhanced pictorial representation of a seismic event, includes a stack 
processor, summation processor, error processor and feedback means. The 
stack processor stacks the common midpoint gather which includes the 
event, using an estimated propagation velocity, to thereby estimate the 
zero offset travel time for the event. The summation processor sums the 
common midpoint gather using first and second weights to provide 
respective first and second weighted sums over offset, using the estimated 
propagation velocity, over a time window centered upon the estimated zero 
offset travel time for the event. The error processor responds to the 
first and second weighted sums, and approximates the error in the 
estimated propagation velocity, to thereby obtain an updated estimate of 
the propagation velocity. The feedback means applies the updated estimate 
of the propagation velocity to the stack processor to thereby restack the 
common midpoint gather in accordance therewith, to thereby provide an 
enhanced pictorial representation of the seismic event. 
The apparatus may further include iteration means for checking to see 
whether the approximate error in propagation velocity is greater than a 
predetermined tolerance, and if so, for applying the updated estimate of 
propagation velocity to the summation processor, to produce third and 
fourth weighted sums which are applied to the error processor, to thereby 
produce a further updated estimate of propagation velocity, for 
application to the stack processor. 
As mentioned in connection with the first aspect of the invention, the 
summation processor can preferably select the first and second weights to 
give the apparatus desired convergence properties. The error processor 
includes means for determining the deviation of the peak of the first 
weighted sum from the center of the time window, relative to the deviation 
of the peak of the second weighted sum from the center of the time window, 
or in accordance with the preferred embodiment, it can crosscorrelate the 
first and second weighted sums and determine the deviation of the peak of 
the crosscorrelated sums from the center of the time window. 
As also mentioned in connection with the first aspect of the present 
invention, the error processor may also approximate the error in the 
estimated zero offset travel time, to produce an updated zero offset 
travel time, and the feedback means can apply the updated zero offset 
travel time to the stack processor to thereby restack the common midpoint 
gather in accordance therewith. Also, the apparatus may include means for 
converting the velocity, as determined above, into true interval velocity 
versus depth or true average velocity versus time, and for converting the 
re-stacked common midpoint gather into a pictorial representation of the 
seismic event at a location substantially corresponding to the actual 
spacial location of the event, using the true interval or average 
velocities, and zero offset times.

DETAILED DESCRIPTION OF THE INVENTION 
The principles of the invention will broadly be described with first 
reference to FIG. 5. Common midpoint data, such as those described with 
reference to FIG. 1, are provided, for a series of adjacent midpoints 
along a seismic line under study, on magnetic tape 10 or other suitable 
data storage device. A velocity estimate, preferably based upon other 
velocity analyses performed on nearby lines, or from sonic logs taken from 
a nearby well, or from a reasonable guess, is manually input to the 
velocity store and calculator 12, by way of keyboard 11, or the like, as 
shown. Based on this initial velocity estimate, the common midpoint data 
from tape 10 are stacked in stack processor 14 and the output of the 
processor 14 is applied to graphics apparatus 16. In this manner, 
zero-offset traces of the data are produced by apparatus 16, an example of 
which is illustrated in FIG. 6, wherein each trace along the time axis 
represents the summed trace, or stack, as shown in FIGS. 3B and 4, of a 
single common midpoint gather of the type shown in FIG. 1. In FIG. 6, each 
trace is darkened in a well-known manner by the graphics apparatus as it 
is deflected along the common midpoint position (or distance) axis, to 
thereby contrast the seismic events against the background. It will be 
noted that when the stack processor stacks the midpoint gathers and places 
adjacent stacks, or zero-offset traces next to one another in the manner 
shown in FIG. 6, seismic events are pictorially represented in two 
dimensions, namely along the common midpoint position axis, and in time. 
Based upon the set of stacks shown in FIG. 6, produced by graphics 
apparatus 16, an event or events of interest are selected for further 
analysis. Based upon the distance of each selected event along the time 
axis, an estimate of the zero offset time t.sub.0, defined as the time at 
which the event occurs at x=0 (time from S.sub.0 to R.sub.0 in FIG. 1), is 
made and manually input to t.sub.0 store and calculator 18 by way of 
keyboard 11, or the like. Alternatively, an event can be selected based on 
the raw common midpoint data, without stacking, and an initial t.sub.0 
estimated therefrom; this alternative is indicated by the data path, shown 
in phantom, from the magnetic tape 10 to the t.sub.0 store and calculator 
18. 
Having estimated the propagation velocity and zero offset travel time 
t.sub.0, the common midpoint data are applied, along with the velocity and 
the t.sub.0 estimates, to summation processor 20 which functions to 
produce a pair of weighted sums of the selected common midpoint gather, 
using the estimated velocity, over a window centered upon the zero offset 
travel time, as better explained with reference to FIG. 7. As shown in 
FIG. 7 the selected common midpoint gather, shifted in time based upon an 
estimated velocity for the selected event, is summed in summation 
processor 20 over a time window centered upon t.sub.0, and two sums are 
produced. The first sum shown is based on a weighted distribution of the 
individual traces in the common midpoint gather which weights the close-in 
traces (small values of offset x) more heavily than, or at least as 
heavily as, the far-out traces (larger values of offset x). The second 
weighted sum weights the far-out traces more heavily than the close-in 
traces. The value G.sub.1 is defined as the time deviation between the 
peak of the first weighted sum from t.sub.0, and the value G.sub. 2 is 
defined as the time deviation between the peak of the second weighted sum 
from t.sub.0. The magnitude of the value of G.sub.2 -G.sub.1 is indicative 
of the error in the velocity used to stack the common midpoint gather. 
Alternatively, the first and second weighted sums can be crosscorrelated 
to produce the crosscorrelated sums as shown, and the quantity G.sub.c, 
defined as the time deviation between the peak of the crosscorrelated sums 
and t.sub.0, is indicative of the estimated velocity error. The 
relationships which yield the velocity error as a function of G.sub.1 and 
G.sub.2, or G.sub.c, as well as a full mathematical explanation of this 
principle, are given in detail in Appendix I, below. 
The first and second weighted sums are delivered from the summation 
processor 20 to error processor 22 which calculates the estimated velocity 
error in accordance with G.sub.1 and G.sub.2, or in accordance with 
G.sub.c alone. The use of the crosscorrelated sums is preferred insofar as 
crosscorrelation by itself produces a signal to noise enhancement. Also, 
when G.sub.c is used, the maximum amplitude of only a single signal (the 
crosscorrelated sums) must be picked, thus further facilitating the 
implementation of this technique by a computer, without human 
intervention. 
In addition to providing an estimate of the error in velocity, the 
technique in accordance with the present invention also provides a value 
associated with the estimated error in the selected zero offset time 
t.sub.0, also as a function of G.sub.1 and G.sub.2, or solely in 
accordance with the value of the G.sub.c. The relationships which yield 
the t.sub.0 error are also provided in Appendix I. 
After having calculated the estimated velocity error in processor 22, this 
error can be used to update the stacking velocity to a value more closely 
approximating the actual velocity. The velocity error can be checked at 
decision processor 24 to see whether it is within a predefined tolerance 
for the stacking velocity. If the initial velocity estimate is already 
within the predetermined tolerance, no further work need be done for the 
event under examination. However, if the velocity error is greater than 
the predetermined tolerance, the velocity error is fed back to the 
velocity store and calculator 12 to provide an updated velocity value. The 
updated velocity is applied from the velocity store and calculator 12 to 
the summation processor 20 whereupon a new pair of weighted sums, or their 
crosscorrelation, are produced in accordance with the new velocity, and 
the new velocity error can be checked, as described above, to see whether 
it is within the predetermined tolerance. In this iterative manner, it has 
been found that the updated velocity will rapidly converge to within good 
tolerances, with real data. The number of iterations required for 
convergence depends upon the closeness of the initial velocity estimate to 
the actual velocity, the signal-to-noise ratio of the seismic data, and 
the weighting schemes employed in the error processor 22, as explained 
below. 
After having determined the velocity to within a predetermined tolerance, 
corresponding to a particular event in a particular midpoint gather, a 
similar process can be performed on other events of interest within the 
same or different gathers, to thereby develop a set of velocities each 
corresponding to a particular event. The common midpoint data can be 
re-stacked in stack processor 14 using the updated velocities and the 
result delivered to the graphics apparatus 16, to portray the seismic 
events using the highly accurate stacking velocities, thereby enhancing 
the representation of those events. FIG. 8 is an illustration of one such 
set of common midpoint stacks using velocities determined to be within 5 
to 15 meters per second of the actual velocities, in accordance with the 
teachings herein. The enhancement of the pictorial representations of the 
events over that of FIG. 6 will be apparent to those skilled in the art. 
As noted above, the error processor 22 can also provide an indication of 
the error in the zero offset travel time t.sub.0, and such error signal 
can be used to revise the zero offset travel time used by summation 
processor 20 and by stack processor 14. Under normal circumstances, 
t.sub.0 will already be accurately determined by a visual inspection of 
the initial stack using the initial velocity estimate, and accordingly the 
value of G.sub.1 will ordinarily be small. However, the use of the t.sub.0 
updating technique can be useful where, for example, old common midpoint 
data are being used and there is some ambiguity as to the time axis used 
in collecting the data. 
In addition to providing more accurately portrayed events along the time 
axis, the events can be portrayed along a depth axis, by converting the 
"stacking" or "apparent RMS" velocities, as determined above, into "true 
interval" velocities versus depth, or "true average" velocities versus 
time, in velocity converter 27. The converted velocity for each event, 
from converter 27, the corrected stack from stack processor 14, and 
optionally the updated zero offset travel time, are then applied to a 
stack converter 26, which converts each zero-offset trace from a time to 
distance axis, to thereby portray each event substantially at its actual 
vertical location. The output of the stack converter 26 is applied to the 
graphics apparatus 16 to produce a visual representation of the events at 
their proper locations. An example of such visual representation is 
illustrated in FIG. 9. 
A suitable technique for converting from stacking velocities to true 
interval or true average velocities in stack converter 27 is given in 
Appendix II, below, and in an article presented by J. W. C. Sherwood et 
al., at the 18th Annual Offshore Technology Conference, Houston, Tex., May 
5-8, 1986, OTC Paper #5161, entitled "Depths and Interval Velocities From 
Seismic Reflection Data For Low Relief Structures", pp. 103-107, (and 
Figures appended thereto), the entire disclosure of which is hereby 
incorporated by reference. 
Returning to FIG. 7, the size of the window centered upon t.sub.0, over 
which the common midpoint data are summed should be large enough to 
contain at least one whole event, to thereby obtain a good signal-to-noise 
ratio, yet no so large as to be cluttered by a large number of events, 
which could hinder convergence to the true velocity, if different events 
within the window become confused for one another. However, under normal 
circumstances, the window can be large enough to include up to about three 
discrete events, with good results, and it has been found that a window of 
approximately 100 milliseconds is adequate in most cases. 
The rate at which the present technique converges to the true velocity is a 
function of the weights that are applied to the common midpoint data to 
produce the first and second weighted sums. With reference to FIGS. 10, 11 
and 12 three different weighting schemes are illustrated. In FIG. 10, the 
first weight W.sub.1, used to produce the first weighted sum, is unity, 
i.e. all signals within the common midpoint gather are weighted equally. 
The second weight, W.sub.2, used to produced the second weighted sum, 
weights the signals within the common midpoint gather in accordance with 
the square of the offset. That is, the far-out offsets are weighted more 
heavily than the close-in offsets. This weighting scheme has been found to 
converge relatively slowly, but has a high immunity against confusion of 
different events within a single time window. 
FIG. 11 illustrates another weighting scheme, whereby the first weight 
W.sub.1, is defined as unity for offsets less than or equal to one-half 
the maximum offset, and zero for offsets greater than one-half the maximum 
offset. Conversely, the second weight W.sub.2 is equal to unity for all 
offsets greater than one-half the maximum offset, and zero elsewhere. It 
has been found that this weighting scheme causes very rapid convergence, 
and is therefore very fast when implemented, but it is more subject to 
possible confusion between events within a single time window. 
FIG. 12 represents a compromise between the weighting schemes of FIGS. 9 
and 10, wherein the first weight W.sub.1 is given by 
EQU 1/2(1+cos (.pi.(x-x.sub.min)/(x.sub.max -x.sub.min))), 
where x equals offset, and the second weight is given by 1-W.sub.1. This 
weighting scheme is believed to converge rapidly while still maintaining 
significant immunity against confusion between events within a single 
window. 
Although explained with reference to the discrete components illustrated in 
FIG. 5 to facilitate the explanation of the present invention, the 
invention readily lends itself to programmed implementation on a digital 
computer, and in order to facilitate an understanding of a computer 
implementation, a generalized program is provided below in Appendix III. 
Thus, the present invention provides a technique for enhancing the 
representation of seismic events in a manner which is much more rapid and 
accurate than formerly possible. Because only those velocities needed to 
perform the iteration are used, the technique of the present invention 
will be faster than traditional techniques, and also more accurate since 
the common midpoint data are re-stacked at least once with a velocity 
within some tolerance of the true velocity, whereas in the traditional 
technique, an estimate of the true velocity is only obtained by 
interpolation. Also, in addition to being appreciably faster than 
techniques previously used, the present invention is more accurate and 
robust in the presence of seismic noise. 
Various changes in variations to the present invention will occur to the 
skilled in the art in view of the foregoing description. It is intended 
that all such changes and variations be encompassed so long as the present 
invention is employed, as defined by the claims appearing below. 
APPENDIX I 
Seismic data for one particular event in a common midpoint gather has the 
following approximate form: 
EQU s(x,t)=f(t-t.sub.x) 
where 
EQU t.sub.x.sup.2 =t.sub.0.sup.2 +x.sup.2 /v.sup.2 
The trajectory defined by the initial estimates of the time and velocity is 
given by: 
EQU (t.sub.x ').sup.2 =(t.sub.0 ').sup.2 +x.sup.2 /v'.sup.2 
In the following v' and t.sub.0 ' will be used to indicate initial velocity 
and zero offset time, and v and t.sub.0 to indicate true velocity and zero 
offset time. The true (or "observed") travel time will be denoted by 
t.sub.x, and the trajectory by t.sub.x '. 
Two weighted stacks of the data over offset x are produced along the 
trajectory described above, using weights W.sub.1 (x) and W.sub.2 (x). 
Mathematically this is expressed by: 
##EQU1## 
To simplify the notation, t' is defined as t+t.sub.x ' and .DELTA.t.sub.x 
as t.sub.x +t.sub.x '. By means of a Taylor series expansion of 
f(t-t.sub.x +t.sub.x ') in terms of f(t), the above equations can be 
simplified. The Taylor series expansion is given by: 
EQU f(t-.DELTA.t.sub.x)=f(t)-(df(t)/dt)..DELTA.t.sub.x +higher order terms (2) 
Dropping the higher order terms and placing Eq. 2 into Eqs. 1a and 1b gives 
the following: 
##EQU2## 
Since neither f(t) nor df(t)/dt are functions of x, they can be taken out 
of the summations, since the summations are over x. Letting Q.sub.1 and 
Q.sub.2 be the sums over x of the first and second weighting functions, 
the following results: 
##EQU3## 
Note that equations 4a and 4b are themselves Taylor series expansions to 
the first order of the following: 
##EQU4## 
Just as picking f(t-t.sub.x) gives t.sub.x, picking 5a will give 
##EQU5## 
and picking 5b will give 
##EQU6## 
The values of these picks will be called G1 and G2. Therefore: 
##EQU7## 
To obtain the updating formula, .DELTA.t.sub.x must be given in terms of 
update velocity and time. Note that .DELTA.t.sub.x, equal to t.sub.x 
-t.sub.x ', is a function of the initial and true velocity and time. Now 
t.sub.x can be related to t.sub.x ' via a Taylor series expansion (again 
dropping higher order terms): 
##EQU8## 
Here w has been used instead of v, where w is equal to 1/v, and is known as 
"slowness" in the art. 
From the equation for the trajectory, these derivatives can be computed and 
are: 
EQU (dt.sub.x ')/(dw)=x.sup.2 /(t.sub.x 'v') 
EQU and 
EQU (dt.sub.x ')/(dt.sub.0 ')=t.sub.0 '/t.sub.x ' 
With these values for the derivatives and by placing Eq. 7 into Eqs. 6a and 
6b, the following is obtained: 
##EQU9## 
Since neither (w-w') nor (t.sub.0 -t.sub.0 ") are functions of x, they can 
be taken out of the summation. Doing this results in: 
##EQU10## 
Introducing the notation X1=.SIGMA.W.sub.1 (x)x.sup.2 /t.sub.x ', 
X2=.SIGMA.W.sub.2 (X) x.sup.2 /t.sub.x ', T1=.SIGMA.W.sub.1 (X) t.sub.0 
'/t.sub.x ', and T2=.SIGMA.W.sub.2 (X) t.sub.0 '/t.sub.x ', results in: 
EQU Q1.multidot.G1=X1(w-w')/v'+T1(t.sub.0 -t.sub.0 ') (10a) 
EQU Q2.multidot.G2=X2(w-w')/v'+T2(t.sub.0 -t.sub.0 ') (10b) 
This equation can be inverted and solved to give the updates in terms of 
known and calculable parameters: 
EQU (w-w')/v'=(T1.multidot.Q2.multidot.G2-T2.multidot.Q1.multidot.G1)/D (11a) 
EQU (t.sub.0 -t.sub.0 
')=(X2.multidot.Q1.multidot.G1-X2.sqroot.Q2.multidot.G2)/D (11b) 
where D=T1.multidot.X2-T2.multidot.X1 
These are the updating equations that are used when the cross-correlation 
technique is not used. 
For the crosscorrelation method, the two weighed stacks given by Eqs. 5a 
and 5b are crosscorrelated: 
##EQU11## 
where .sym. represents the crosscorrelation operation. The mathematical 
definition of crosscorrelation is: 
##EQU12## 
Picking 12 results in G1-G2, which will be called G.sub.c. 
##EQU13## 
Equation 13 results from subtracting Eq. 10b from Eq. 10a after first 
dividing by Q.sub.1 and Q.sub.2. 
Since two equations have been combined, there are now two unknowns by only 
one equation; for a solution, another equation must be supplied. It will 
be assumed that: 
EQU .SIGMA..DELTA.t.sub.x =0 (14) 
This assumes that the initial time estimate, t.sub.0 ' was taken from an 
initial stack formed using the initial velocity estimate. 
Defining 
##EQU14## 
and 
##EQU15## 
and from Eqs. 14 and 7: 
##EQU16## 
Eqs. 15 and 13 can be inverted to give: 
##EQU17## 
Equations 18a and 18b are the updating relationships used when the 
crosscorrelation technique is used. 
APPENDIX II 
This appendix gives a general method of converting from stacking velocities 
(Vrms) versus time, to interval velocities (Vint) versus depth. 
T 1 Convert Vrms.sup.2, and to, to tx: 
EQU tx=.sqroot.to.sup.2 +X.sup.2 /V.sup.2 rms 
tx is a function of common midpoint (y), and offset (X) as well as layer. 
T 2 Develop a relationship between Vint, z and tx (known as 
ray-tracing). This relationship gives tx for known Vint, z values. 
Therefore this relationship must be inverted. 
EQU tx=R(Vint,z) (1) 
T 3 Invert the relationship. 
In general, R is non-analytical and difficult to invert. Therefore, an 
approximation to R (R') is used instead: 
EQU tx=R(Vint,z) 
EQU tx'=R(V'int,z') 
where V'int, z' are initial estimates of Vint, z. 
Therefore, tx-tx'=R (Vint, z)-R (V'int, z) which equals R'(Vint-V'int, 
z-z'), where Vint, z are updated estimates for Vint, z. R' is an 
invertible approximation for R. 
thus: 
EQU (Vint-V'int,z-z')=(R'.sup.T R').sup.- R'.sup.T (tx-t'x) (2) 
T 4 Iteration: Vint becomes V'int, z becomes z' and Eq. 1 is used to 
generate a new t'x. New Vint and z are obtained from Eq. 2. 
The process is repeated until (Vint-V'int) and (z-z') falls below a 
predetermined amount (say 0.005 z'int or 0.005 Vint'initial). 
The exact form of R' can be determined in accordance with principles known 
in the art, such as those discussed in the above-identified article by 
Sherwood et al. 
APPENDIX III 
This appendix presents the generalized or "pseudo" code of that portion of 
the invention in which velocity is analyzed. The remaining portions of the 
invention, such as stacking and depth conversion use wellknown code and 
are therefore not presented here. 
Definitions 
(a) constants 
NCMP is the number of common midpoint (CMP) gathers to be analyzed 
NVEL is the number of layers to be analyzed per CMP gather 
NSAMP is the number of time samples per trace 
NTR is the number of traces per CMP gather 
DT is the interval between time samples in milliseconds 
NITERS is the maximum number of iterations to perform 
TOLV is the fractional tolerance for convergence to occur 
NDOM is the number (odd) of samples defining the width of the time window 
(b) arrays 
Traces is a 2-D array containing all the traces and samples of a single CMP 
gather 
RANGES is a 1-D array consisting the shot-receiver distances of the CMP 
gather to be analyzed 
WT1 is a 1-D array containing the first weights to be applied to the data 
WT2 is a 1-D array containing the second weights to be applied to the data 
VP is a 2-D array containing the initial estimates of velocity for all CMP 
locations and layers 
TOP is a 2-D array containing the initial estimates of zero offset travel 
time for all layers and CMP locations 
V is a 2-D array containing the final estimates of velocity for all CMP 
locations and layers 
T is a 2-D array containing the final estimates of zero offset travel time 
for all layers and CMP locations 
The rest of the constants and arrays are implicitly defined by their usage 
is the pseudo code 
______________________________________ 
COPYRIGHT .COPYRGT. 1987 ARAMCO 
ALL RIGHTS RESERVED 
______________________________________ 
c start of the code: 
c 
DO ICMP = 1,NCMP 
c get a new CMP gather: 
CALL INPUT (TRACES,NTR,NSAMP,RANGES) 
DO IVEL = 1,NVEL 
c start new layer. Set up starting velocity and time, and 
c calculate tolerance: 
TSTART = TOP(IVEL,ICMP) 
VSTART = VP(INVEL,ICMP) 
TOVLI = TOLV*VSTART 
DO ITER = 1,NITERS 
c calculate weighted stacks (S1 & S2) and various parameters 
(TO, T1, T2, X0, X1, X2, Q1 & Q2): 
CALL GENWTS(TRACES,TSTART,VSTART,NTR, 
NSAMP, 
*NDOM,DT,WT1,WT2,S1,S2,T0,T1,T2,X0,X1,X2, 
Q1,Q2) 
c cross correlate S1 and S2 to get XC: 
CALL CROSSC(S1,S2,XC,NDOM) 
c pick XC to get PICKT: 
CALL PICK(XC,PICKT,NDOM) 
c calculate new velocity and time: 
CALL VELUP(VNEW,TNEW,VSTART,TSTART, 
PICKT,T0, 
* T1,T2,X0,X1,X2,W1,W2) 
c check to see if convergence has occurred: 
IF(ABS(VSTART-VNEW).LT.TOLV1) GO TO FINV: 
c convergence has not occured. Replace old estimates with 
new ones and repeat process: 
VSTART = VNEW 
TSTART = TNEW 
END DO ITER 
c convergence has occurred. Record the new velocities and 
start on new analysis location: 
FINV: TO(IVEL,ICMP) = VNEW 
V(IVEL,ICMP) = TNEW 
END DO IVEL 
END DO ICMP 
END PROGRAM 
c 
c subroutines (explanations or code): 
c 
(1) SUBROUTINE INPUT(TRACES,NCMP,NSAMP, 
RANGES) 
takes data from tape and places it into the 2-D array 
TRACES. 
takes shot receiver distance information and places it 
into 1-D array RANGES. 
(2) SUBROUTINE GENWTS(TRACES,TSTART, 
VSTART,NTR, 
* NSAMP,NDOM,DT,WT1,WT2,S1,S2,T0,T1,T2,X0, 
X1,X2,Q1,Q2) 
shifts data from traces into S1 and S2 after applying 
weights calculates T0,T1,T2,X0,X1,X2,W1,W2 
c code of this subroutine: 
IDOM = NDOM/2 
c set initial values of S1 and S2 and other parameters to zero: 
DO I = 1,NDOM 
S1(I) = 0.0 
S2(I) = 0.0 
END DO I 
TO = 0.0 
T1 = 0.0 
T2 = 0.0 
X0 = 0.0 
X1 = 0.0 
X2 = 0.0 
Q1 = 0.0 
Q2 = 0.0 
c iteratively calculate S1 and S2 and the other parameters: 
DO ITR = 1,NTR 
X = RANGES(ITR) 
c calculate trajectory: 
TSP = SQRT(TSTART**2 + (X/VSTART)**2) 
DLT = TSP/DT + 1 
T1 = TSP - IDOM 
T2 = T1 + NDOM - 1 
DO IT = T1,T2 
S1(IT) = S1(IT) + WT1(ITR)*SHIFT(TRACES 
(1,ITR),DLT) 
S2(IT) = S2(IT) + WT2(ITR)*SHIFT(TRACES 
(1,ITR),DLT) 
END DO IT 
T0 = T0 + TSTART/TSP 
T1 = T1 + WT1(ITR)*TSTART/TSP 
T2 = T2 + WT2(ITR)*TSTART/TSP 
X0 = X0 + X*X/TSP 
X1 = X1 + WT1(ITR)*X*X/TSP 
X2 = X2 + WT2(ITR)*X*X/TSP 
Q1 = Q1 + WT1(ITR) 
Q2 = Q2 + WT2(ITR) 
END DO ITR 
END SUBROUTINE 
SUBROUTINE SHIFT(TRACES,DLT) 
returns TRACES(DLT) 
because DLT is not an integer, interpolation between 
samples is usually necessary 
SUBROUTINE CROSSC(S1,S2,XC,NDOM) 
cross correlates S1 and S2 and places result into XC. 
S1, S2, and XC all contain NDOM samples 
SUBROUTINE PICK(XC,PICKT,NDOM) 
picks maximum amplitude of XC and puts PICKT equal 
to the time deviation of the peak from the center of 
the window 
SUBROUTINE VELUP(VNEW,TNEW,VSTART, 
TSTART,PICKT, 
* T0,T1,T2,X0,X1,X2,Q1,Q2) 
using the mathematical relationships developed in 
Appendix I, 
VNEW and TNEW are obtained from VSTART, 
TSTART, PICKT, 
and the parameters TO,T1,T2,X0,X1,X2,Q1 and Q2. 
c 
c this ends the psuedo code. 
c 
______________________________________