CHAPTER TWENTY-THREE: SEISMIC PETROPHYSICS 3 Synthetic Seismograms

Table of Contents
23.00: Introduction To This Chapter
23.01: Seismic Data Acquisition
23.02: Seismic Velocity Analysis
Techniques
23.03: Calculating Seismic Velocity From Logs
23.04: Calculating Interval Velocity From Seismic Data
23.05: Time to Depth Conversions
23.06: Seismic Velocity and Dipping Beds
23.07: Seismic Modeling Concepts
23.08: Step By Step Procedure for Modeling
23.09: Integrating the Sonic Log
23.10: Acoustic Impedance and Reflection Coefficients
23.11: Calculating Seismic Wavelets
23.12: Convolution of a Wavelet with the Reflection Coefficients
23.13: Case Histories: Synthetic Seismograms
1. Salt Solution
2. Devonian Reef
3. Granite Wash Sandstone
4. Stratigraphic Trap
23.14: In Conclusion
23.15: Exercises for Chapter Twenty-Three
23.16: Bibliography for Chapter Twenty-Three

Continue to Chapter Twenty-Four

Publication History: Portions of Sections 23.01 through 23.13 were originally published in "Review of Semi-Automatic Velocity Analysis Techniques" by E. R. Crain, ASEG Bulletin, Sydney, 1971 and in various technical bulletins prepared for Digitech Ltd, 1971-1973. Section 23.14 is based on training material prepared in 1971 by R. O. Lindseth, a pioneer in seismic inversion processing. Although these sources are more than 30 years old, the basic theory has not changed - only speed, quality, and implementation have improved dramatically.

This Chapter also formed part of Chapter Ten of Volume Two of The Log Analysis Handbook, a self published series of course notes covering geological and geophysical aspects of log analysis. First published in 1978, revised 1985, 1993. Completely revised and reorganized for this electronic edition Sep 2002.

CHAPTER TWENTY-THREE: SEISMIC PETROPHYSICS 3 Synthetic Seismograms

23.00 Introduction To This Chapter
This Chapter continues the discussion of Seismic Petrophysics with explanations and case histories of synthetic seismograms, seismic inversion, vertical seismic profiles, and amplitude versus offset studies. These are mainly the activities of geophysicists and seismic interpreters, but log data of the kinds described in Chapter Twenty-One are often used to calibrate or aid interpretation. Petrophysicists and geophysicists should both have a basic understanding of these techniques, so some brief description of seismic techniques are included to bridge the gap between the two disciplines.

Seismic petrophysics is a term used to describe the conversion of seismic data into meaningful petrophysical or reservoir description information. Until recently, this work was qualitative in nature, but as seismic acquisition and processing have advanced, the results are becoming more quantitative. Mapping of seismic attribute or inversion data to obtain lithology, porosity, or fluid content are now widespread. Calibrating this work to “ground truth” can convert the relatively fuzzy seismic attributes into useful reservoir exploration and development tools. A seismic petrophysics study aimed at quantifying porosity is shown in Figure 23.01.

FIGURE 23.01: Seismic petrophysics study for porosity

“Ground truth” takes the form of well logs, or log analysis results, which in turn are calibrated to core, well test, and production data. Unfortunately, it takes a fair amount of work to compare seismic results to log data. The logs will usually require some kind of editing or modeling or both. Comparison of seismic results to log data may indicate that further processing of the seismic is needed, and the calibration cycle is repeated - often several iterations are needed

23.01 Seismic Data Acquisition
Seismic data originates with a source of impulsive acoustic energy near the surface of the ground. Dynamite, air guns, or vibrating plates are used on land or marine surveys. In addition, a falling weight is sometimes used for shallow land surveys and a sparker survey is sometimes used for shallow marine work. The acoustic energy passes through the earth in all directions. Some is reflected back to the surface by acoustic impedance barriers, as shown in Figure 23.02 (top). Here, the returned energy is sensed by geophones planted on the ground or floating near the surface of the water.

FIGURE 23.02: Geometry of seismic acquisition

Some reflected energy bounces back and forth more than once. These events are called ghosts if they occur in the near surface, and multiples if they come from deeper reflectors. Multiples and ghosts are a form of interference which is usually eliminated by suitable data processing.

Each geophone signal is recorded on digital magnetic tape or disc and presented as a wiggly trace of energy amplitude versus arrival time. Raw traces are seldom delivered as the final product. Considerable data processing is performed to correct for geometry, the filtering effect of the earth, and amplitude decay with depth.

In 2-D seismic, the source and geophones are located in a straight line, resulting in a seismic cross section. If the line cannot be straight due to topography, the data is processed to collect data in short approximations to straight lines.

For 3-D seismic, receivers and sources are set up in a pattern which allows simultaneous recording of many intersecting lines of data. These can be processed to provide a volumetric view of the subsurface. An illustration of 2-D and 3-D geometry is shown in Figure 23.02.

4-D seismic is a term used to describe surveys taken on the same grid several years apart and are used to show changes in reservoir properties over time. These can only be due to changes in fluid content from production or injection. The results are used to evaluate production efficiency, the effectiveness of waterfloods, or monitor aquifer influx.

4-C seismic is a relatively new form of marine survey and refers to four component recording of the seismic signal. The components are the usual compressional or P-wave from a geophone, plus in-line and cross-line shear arrivals, as well as a compressional wave recorded on a hydrophone. The different response of the geophone and hydrophone to reverberations in the water allow specialized processing to remove interference.

Anywhere from 24 to several hundred surface points are recorded for each surface shot. The energy source is then moved a multiple or sub-multiple of the geophone spacing and the signal recorded again. If the movement is less than half the geophone spread length, then the same subsurface points will be recorded more than once, resulting in multi-fold coverage. The seismic traces from different surface layouts that fall at common depth points are collected and stacked together to improve signal to noise ratio. The improvement in data quality with increased coverage is shown in Figure 23.03.

FIGURE 23.03: Multi-fold (stacked) seismic sections

The presentation of seismic sections has evolved over the years, from plain wiggly traces to variable intensity black and white or color displays, as in Figure 23.04. The color can represent signal amplitude, frequency content, or any other desired (and determinable) property of the seismic signal. 3-D seismic is often presented in color as vertical or horizontal slices, as isometric views, or as contour maps, as in Figure 23.05.

FIGURE 23.04: Seismic presentation options

FIGURE 23.05: Mapping, slicing, and 3-d isometric views of seismic data

Most log analysts are unfamiliar with these displays, making it difficult for them to communicate well with geophysical interpreters.

23.02 Seismic Velocity Analysis Techniques
Seismic ray paths, as shown in Figure 23.02, are rarely straight lines, but rather follow the laws of refraction, resulting in curved ray paths. The exact shape of the path is determined by the rock velocity distribution versus depth and the spread geometry. In Figure 23.06, the apparent curvature of the seismic reflections is caused by these factors.

Sound velocity in the subsurface can be found directly from the seismic data by observing the difference in arrival time of reflectors versus distance from the source. The time difference is called normal moveout or NMO.

FIGURE 23.06: Raw 24 trace seismic record showing normal moveout

A large number of semi-automatic techniques have been developed to derive seismic velocity from seismic data independently of well log data, Figure 23.07 and 23.08 show two types of velocity analysis display. Many others have been developed over the years. The principle of such techniques is based on computerized approximations to normal moveout analysis which used to be done by hand calculation. The normal moveout is often called delta-T, just as the sonic log interval transit time is called delta-T by log analysts (abbreviated DELT in this book), so we will always use the term normal moveout (abbreviated NMO) to avoid confusion. My first job as a geophysical engineer in 1967 was to plot NMO versus depth, pick the slopes of the graph and enter the data into a TIAC computer to calculate velocity.

FIGURE 23.07: Velocity analysis display -velocity panels

FIGURE 23.08: Velocity analysis display cross correlation panels

23.03 Calculating Seismic Velocity From Logs
If the interval thicknesses and interval velocities are given, for example by a digitized sonic log or a seismic model of a hypothetical rock sequence, we can calculate what the seismic times would be at various reflectors, recorded at detectors spaced along a geophone spread. When this model is plotted to scale with all major reflectors and the resulting ray paths, it is called ray tracing. The equations work for both shear and compressional waves when the respective interval velocities are used.

The near trace time to a reflector is:
_____1: To = 2 * Sum (Hi / Vi)

The interval travel time in a layer is:
_____2: Ti = 2 * Hi / Vi

The average velocity is defined as:
_____3: Vavg = 2 * Sum (Hi) / To

The RMS velocity is defined as:
_____4: Vrms = ((Sum ((Vi ^ 2) * Ti)) / To) ^ 0.5

The far trace time, from ray path geometry is:
_____5: Tx = (2 * ((X / 2) ^ 2 + (Sum (Hi)) ^ 2) ^ 0.5) / Vrms

Normal moveout is:
_____6: NMOc = Tx - To

These equations are used to evaluate various layered models (ray tracing), and the reverse equations (shown below) are used to calculate interval velocity from seismic data. The equations apply to both compressional and shear data, when appropriate inputs are used.

Note that seismic times are always two way times, and that integrated sonic log times are usually one way times.

23.04 Calculating Interval Velocity From Seismic Data
If we are interested in finding interval velocity from seismic data, instead of from a sonic log, we must use equations that represent the physics of the seismic recording process. Since the geometry of the reflecting horizons is unknown, we must make some assumptions that may later turn out to be untrue. The equations work for both shear and compressional waves when the respective two way times are used.

The far trace time is:
_____1: Tx = To + NMO

Stacking Velocity is:
_____2: Vstk = (X ^ 2 / (Tx ^ 2 - To ^ 2)) ^ 0.5

Based on the geometry shown in Figure 23.09, interval velocity between any two points in the section is:
_____3: Vint = (((Vstk2 ^ 2) * T2 - (Vstk1 ^ 2) * T1) / (T2 - T1)) ^ 0.5

FIGURE 23.09: Geometry for interval velocity calculation

The apparent velocity, or stacking velocity (Vstk), is the velocity which yields the exact NMO from the NMO equation. Therefore, it could be derived from real seismic data with a T^2 - X^2 plot, an NMO analysis, or a computer velocity analysis (CVA). It is the only velocity which will yield the best stacked seismic section.

The average velocity from a log analysis will yield very poor stacking results unless the interval velocities are quite constant and there is no dip on the reflections. The RMS velocity from log analysis will yield adequate results, as long as dip is not extreme. In horizontal bedding, we usually assume that Vstk from seismic equals Vrms from log analysis. Thus NMO could be predicted from well log data by using the appropriate arithmetic. Conversely, stacking velocity can be transformed into interval velocity, which will lead to correct time to depth conversion.

23.05 Time to Depth Conversions
If we assume that we have a list of seismic times and corresponding stacking velocities (Vstk) from a constant velocity stack, or similarly derived seismic velocity analysis, we can calculate the interval velocity between each velocity pick, provided we also assume Vrms = Vstk. Average velocity and true depth are then computed from the equations previously given, and presented in tabular or in graphical form. The interval velocities will compare with sonic log data, if the seismic processing, seismic data, and sonic log data are of good quality.

Straight ray depth is:
_____0: Dstr = Vstk1 * T1 / 2

Straight ray depth is easy to calculate but is NOT the correct way to calculate depth. Unfortunately, this conversion is still in common use. Depths predicted by this method will be too deep by several hundred feet in a 5000 foot well. You can compute it if you want to, but DON'T USE IT.

Interval velocity (as derived previously) is:
_____1: Vint = (((Vstk2 ^ 2) * T2 / 2 - (Vstk1 ^ 2) * T1 / 2) / (T2 - T1)) ^ 0.5

Interval thickness is:
_____2: Hi = Vint * (T2 - T1) / 2

True depth with curved ray path is:
_____3: Dtrue = Sum (Hi)

And the average velocity is:
_____5: Vavg = 2 * Sum (Hi) / T1

The straight ray depth is always deeper (larger value) than the true depth, so this simplification should not be used for time to depth conversions. A typical graphical display of this kind of data is shown in Figure 23.10.

FIGURE 23.10: Velocity analysis results

23.06 Seismic Velocity and Dipping Beds
Another typical problem in common depth point shooting is shown in Figure 23.11. Ray paths are shown only from reflector 3 for clarity. Normal incident points are those points on reflecting horizons which would be recorded from a common ground point, that is from a source and a receiver in the same position.

FIGURE 23.11: Stacking common depth point problem with dipping beds

When we do a velocity analysis on a computer, we compare the arrival times of the incident points from the reflector. For example, from reflector 3, we measure the change in NMO from the reflection which originates at C with that which originates at B and at A. The computer derived velocity is the velocity which will take the three reflecting points, C, B, and A, and stack the reflections from them in phase.

Points A, B, and C are not common depth points, but that doesn't matter for stacking purposes. The computer derived velocity analysis forces the time of the stacked trace to be that of the zero offset or normal incident point C. The presence of dip, as on reflection 3, will obviously reduce the observed NMO between C and A. Therefore we can conclude in this example that as dip increases, we will have less NMO and the apparent velocity will be higher.

It is clear, when using seismic velocity analyses in the absence of migration, that the computed interval velocities may differ very greatly from the true interval velocities. It is just as important to use these velocities for stacking, because a true velocity or an average velocity derived from interval velocities will not move the reflections from points A, B, and C to a common point. Only an apparent, or stacking, or computer derived velocity will do the proper stacking job when beds dip appreciably.

For structural or stratigraphic interpretation, the computer derived seismic velocity functions must be reworked to compensate for dip and the effect of ray path geometry before being used. The correction procedure is called seismic migration. In Figure 23.11, Reflector 3 appears to be at a position given by line G-H because the vertical times Y-H and X-G are equal to the normal incidence times Y-K and X-J respectively. The steeper the dip, the worse the discrepancy becomes. Since we see line G-H on the seismic section, we underestimate dip and depth to the actual reflecting horizon.

To correct this, we note that:
_____1: Sin (DIPact) = Tan (DIPapp)

Where:
__DIPact = dip of actual reflector (line G-H, degrees)
__DIPapp = dip of apparent reflector (line J-K, degrees)

By using the power of the computer, we can identify coherent reflectors, compute their apparent dips, find their true dips, and migrate all reflections from their apparent to their true positions. This is not a trivial task. When dip direction is not the same as the line of the seismic section, migration must also involve dip rotation into the line of section, referred to as 3-D migration. If seismic interval velocity does not match the sonic log, formation dip may be the reason.

23.07 Seismic Modeling Concepts
Seismic modeling is a loosely defined term. It has been taken to include any or all of the following:

1. compute seismic response from a postulated rock sequence, using assumed velocity and density values for successive layers

2. compute seismic response from unedited well logs (sonic or density or both)

3. compute seismic response from modeled and edited log data values to reflect real or hypothetical fluid, porosity, shale, and matrix rock quantities or types

The latter form of modeling is by far the most successful, but it requires an extra step - quantitative log analysis and log reconstruction. Synthetic seismograms will NOT be adequate unless this extra work is done.

There are five major reasons why log editing or log modeling may be necessary:

1. large or rough boreholes that prevent accurate logs from being recorded (Eg: cycle skips on sonic logs)

2. invasion of drilling fluid into gas or light hydrocarbon zones

3. vuggy or isolated porosity types

4. rock alteration by the drilling process

5. missing density or sonic information (not recorded or tool not working properly)

There are a number of different techniques commonly used to overcome bad or missing data. These are discussed below.

23.08 Step By Step Procedure for Seismic Modeling
We will start with four basic definitions:

1. forward seismic modeling - making a synthetic seismic trace from EDITED sonic and density log data.

2. inverse seismic modeling - making a synthetic acoustic impedance log from a seismic trace.

3. seismic interpretation - making correlations, picking horizons, and mapping seismic data, with geological and well log data for control.

4. modeling a log - calculating what a log should read in a given rock and fluid mixture, including the creation of synthetic sonic and density logs from other logs or geological/seismic data.

Modeling of logs to repair recording problems and to create logs for hypothetical rock and fluid mixtures is covered in Chapter Twenty-One.

Just as log analysts use core data as ground truth to calibrate their log analysis results, geophysicists use synthetic seismograms created from edited log data as their ground truth.

Structural interpretation of seismic sections involves identification of reflective horizons and picking of seismic travel times for each trace on each horizon. From these data, time, and subsequently depth, maps are made for each horizon. Choosing the exact horizon time is sometimes difficult, since the bed boundary may not be on a peak or a valley of the seismic trace. A synthetic seismogram is often used to calibrate these time picks.

Stratigraphic interpretation is based on the ability to correlate seismic character to subsurface geology, instead of just seismic two-way times. The correlation is established at a well with a synthetic seismogram computed from edited well log data. The shape and amplitude of the reflection on each trace is important in its interpretation.

The steps in making a synthetic seismogram are:

1. edit the sonic and density log for borehole and recording problems, based on regional trends, offset logs, and mathematical models of log response

2. model sonic and density logs in formations which have been affected by invasion or rock alteration, based on a comprehensive quantitative log analysis

3. model effects in carbonates caused by porosity type or density contrast, based on log analysis and geologic data

4. integrate the modeled and edited sonic log

5. interpolate equal time increment values for sonic and density (and other log) values from depth data

6. calculate acoustic impedance and reflection coefficients from modeled and edited logs

7. generate an appropriate wavelet

8. convolve wavelet with reflection coefficients

9. plot synthetic seismogram on an appropriate time scale

10. check results against real seismic data

11. revise edits or log models over intervals that do not match real seismic, OR improve seismic processing, OR change wavelet characteristics

12. make "What-if" models to test alternate interpretations and sensitivity to fluid, porosity, and lithology assumptions, as well as wavelet shape and frequency

13. remodel zones which do not tie as to time, amplitude, or character, and check again

If possible, the wavelet to be used is based on frequency and phase measured from the seismic data we are trying to match. A seismic wavelet can be extracted from real seismic data by autocorrelation of a seismic trace. If this is not done, a wavelet can be generated from various mathematical expressions, or by filtering a unit impulse with a band pass filter.

Synthetic seismograms derived from unedited and un-modeled data are common commercial products from all vendors. Two typical examples are shown in Figure 23.12. They are useful for gross correlation and major reflector identification. Most vendors allow the purchaser to customize the product to include the type of editing and modeling described in the following sections. However, the responsibility for the model parameters rests with the client as few suppliers have the log analysis skills to prepare adequate log models.

FIGURE 23.12: Commercial synthetic seismogram with integrated sonic log and interval velocity table

The unedited synthetic of Figure 23.13 demonstrates the serious mismatch that can occur. The strong events at the top of the Cretaceous and in the Lower Jurassic are artifacts of bad hole condition affecting the sonic log. Small mis-ties are caused by these problems and other rock alteration. Since I personally made this synthetic in 1969, I feel comfortable criticizing it now, but I thought it was a good match at the time. Today, it would be inadequate and could be improved.

FIGURE 23.13: Synthetic seismogram with corresponding marine seismic section

Synthetics and seismic models can also include multiple reflections (Figure 23.14) or shear wave models from shear wave sonic logs (Figure 23.15). Normal moveout is rarely modeled on a synthetic, but it might be useful to identify reflectors on raw seismic records (Figure 23.16).

FIGURE 23.14: Synthetic seismograms with multiple reflections

FIGURE 23.15: Synthetic seismograms from shear sonic data

FIGURE 23.16: Synthetic model of normal moveout

Using the log response equation, discussed in Chapter Twenty-One, we can inject new layers, delete existing layers, replace existing layers with new ones, change the fluid content from water to gas or oil (or vice versa). This is called "What-if" modeling. Models of this type will show why some water zones have bright spots, why some coal beds have no reflections, and why some carbonate porosity is visible on seismic (and some is not). What-if models allow the interpreter to test different structural and stratigraphic solutions against the actual data before committing to one interpretation.

Another variant of this modeling method is to replace data, or add data to the bottom of one well, using log data from another well which more closely represents the interpreted seismic section. This data must also be modeled and edited to reflect the unaltered rock/fluid mixture. This is often done to aid in design of a sidetrack or whipstock well to find the fault, reef, or salt dome structure that was missed by the original hole.

As far as can be determined, FEW existing seismic modeling workstation performs the necessary log analysis and log recalculation. Therefore YOU must do it first, offline from the modeling package, and enter the edited/modeled logs into the system.

23.09 Integrating the Sonic Log
The reflection coefficient set that is calculated from sonic and density log data could be a string of numbers versus depth. This is because logs are recorded versus depth. However, we need a string of reflection coefficients versus two way seismic time in order to make a synthetic seismogram. To translate depth to two way time, we use a time versus depth plot made from an edited, integrated sonic log, or a computer representation of such a plot.

Integration is a summation of the sonic log readings taken at equal depth increments. This is often adjusted to a datum depth or time horizon, not necessarily the surface. Because the sonic log depth is measured relative to the surface but cannot often be recorded all the way to the surface, we also have to estimate or tie the sonic integrated time to a known horizon below the surface casing. The checkshot survey plays an important role in tying the sonic to surface or some other datum.

The formula is:
_____1: T2way = Tsurf - Tdatum + 2 * Sum (DELTcor * INCR)

Where:
__Tsurf = Two way time from surface to start of sonic log (ms)
__Tdatum = Two way time from surface to desired datum (ms)
__DELTcor = Edited sonic log reading adjusted to SRS or VSP (us/m or us/ft)
__INCR = Digitizing increment (meters or feet)

Once the integration is finished, a depth and acoustic impedance can be calculated for each two way time, by interpolating linearly between two digitized data samples. For more precision, a spline interpolation can be used. The reflection coefficient string is then calculated from the interpolated data. The sonic and density data should always be edited before integration.

The depth digitizing increment should be fine enough to allow reasonable interpolation between the time sample points. The time sample rate is chosen to include the highest frequency (thinnest) events desired on the synthetic. Sample rates from 0.25 to 8 ms are used, with 1 ms the most common. Depth sample rate to match the time sample rate depends on the rock velocity. Too fine a depth sample rate wastes computer time and storage space; too coarse gives poor results.

Some sonic logs have integrated travel time indicated by tick marks along one margin of the log. These are ONE WAY travel times. Hole volume integration may also be indicated by tick marks, and can be confused with sonic integration ticks. These ticks were made before editing, so they are seldom used directly in making a synthetic seismogram. Some integration ticks are recorded incorrectly, and each case should be checked against the actual log data to verify the validity of the integration.

23.10 Acoustic Impedance and Reflection Coefficients
Sound is reflected back toward the source of energy whenever an acoustic impedance boundary occurs or Poisson's ratio changes. Acoustic impedance is the product of velocity and density. Energy is also lost due to reflection and spherical divergence. These are calculated from:

For vertical incidence :
_____1: AcImp1 = KD4 * DENS1 / (DELTc1 * KS3)
_____2: AcImp2 = KD4 * DENS2 / (DELTc2 * KS3)
_____3: Refl = (AcImp2 - AcImp1) / (AcImp2 + AcImp1)
_____4: Atten = Prod (1 - Refl ^ 2)

Where:
__KD4 = 1000 for Metric units
__KD4 = 10^6 for English units
__KS3 = 1.00 for Metric units
__KS3 = 3.281 for English units

For non-vertical incidence:
_____1: K = (Vavg - Vo) / DEPTH
_____2: ANGLE = Arctan ((DEPTH * X + Vo * X/K)/(DEPTH^2 + 2*Vo*DEPTH / K - X^2 / 4))
_____3: Vrat = Vc2 / Vc1
__OR 3a: Vrat = DELTc1 / DELTc2
_____4: Drat = DENS2 / DENS1
_____5: C = (Vrat^2 + (1 - Vrat^2) / (Cos(ANGLE))^2) ^ 0.5
_____6: Refl = (1 - Vrat * Drat * C) / (1 + Vrat * Drat * C)

The reflection coefficient will vary with incidence angle, equivalent to a variation with offset distance. Attenuation is seldom applied to reflection coefficient data, as synthetics are often compared to gain equalized data, in which attenuation has been compensated.

If density log data is missing or cannot be used due to bad hole conditions, an appropriate constant value or a value derived from the empirical chart in Figure 23.17 can be used. A complete reconstruction of the density log can be made if lithology is known, by using the modeling equations given in Chapter Twenty-One.

FIGURE 23.17: Acoustic impedance from velocity

23.11 Calculating Seismic Wavelets
If a wavelet can be extracted by autocorrelation of a real seismic trace, it should be used to make the synthetic. If this cannot be done, wavelets are generated from equations which describe the frequency content of the wavelet. There are two types - Ricker wavelets are generated directly in the time domain and Klauder wavelets are generated in the frequency domain. Both are called zero phase wavelets and are a minimum delay. That is, the maximum energy is at the beginning of the wavelet. Ricker wavelets are shown in Figure 23.18.

FIGURE 23.18: Standard computed wavelet shapes

There are four different Ricker wavelets, called Ricker Type 1, Type 2, Type 3, and Type 4 respectively. The higher types have more peaks and valleys, representing more "ringing", and are less symmetrical than Type

Ricker wavelets are used when the synthetic seismogram is to match dynamite or airgun data. They are computed from the following:

Ricker Type 1:
_____1: G = 3.4641 * CFREQ / 1000
_____2: FOR I = 1 TO WVLTLENGTH / SAMPRATE
_____3: X = (I * SAMPRATE - 1) * G - 3.75
_____4: Wvlt [I] = (X * X - 1) / EXP (X * X / 2)
_____5: NEXT I

Ricker Type 2:
_____1: G = 3.65 * CFREQ / 1000
_____2: FOR I = 1 TO WVLTLENGTH / SAMPRATE
_____3: X = (I * SAMPRATE - 1) * G - 3.85
_____4: Wvlt [I] = (X * X - X - 1) / EXP (X * X / 2)
_____5: NEXT I

Ricker Type 3:
_____1: G = 2.9 * CFREQ / 1000
_____2: FOR I = 1 TO WVLTLENGTH / SAMPRATE
_____3: X = (I * SAMPRATE - 1) * G - 3.85
_____4: Wvlt [I] = (-X^3 + 3 * X) / EXP (X * X / 2)
_____5: NEXT I

Ricker Type 4:
_____1: G = 2.35 * CFREQ / 1000
_____2: FOR I = 1 TO WVLTLENGTH / SAMPRATE
_____3: X = (I * SAMPRATE - 1) * G - 3.2
_____4: Wvlt [I] = 0.1*(-X^6+X^5+15*X^4-10*X^3-45*X*X+15*X+15)/EXP(X*X/2)
_____5: NEXT I

When comparing synthetic seismograms with different wavelets, it is important that each wavelet have the same initial energy. The energy can be found by squaring and summing each of the Wvlt[I] terms. The sum is then compared to a standard energy of 1.0 and a scale factor derived, which is multiplied against each Wvlt[I] term to obtain normalized wavelet values.

The Klauder wavelet is generated in the frequency domain, which gives the amplitude and phase spectra of the wavelet. The autocorrelation of the spectra produces the wavelet shape in the time domain. This wavelet is used when the synthetic is to match Vibroseis data.

Klauder Wavelet:

STEP 1: Generate chirp signal amplitude and phase
_____1: CFREQ = 0.5 * (LFREQ + HFREQ)
_____2: KLFREQ = (HFREQ - LFREQ) / SWEEPLENGTH
_____3: K = SAMPRATE / 1000
_____4: FOR I = 1 TO SWEEPLENGTH / K

Calculate phase
_____5: I[I] = 2 * PI * (LFREQ * I * K + KLFREQ * I * I * K * K / 2) + INITPHASE

Calculate amplitude
_____6: R[I] = 1250 * COS(I[I])

Put phase in all 4 quadrants
_____7: I[I] = TAN (I[I])
_____8: I[I] = ATN (I[I])
_____9: I[I] = I[I] * 180 / PI + 180 * (SIN (I[I]) < 0) + 90
_____10: NEXT I

STEP 2: Calculate autocorrelation to get wavelet in time domain. This double nested loop takes the first 100 data points and cross correlates them with the last 100 data points of the sweep. Since the number of data points is small, there is little advantage to using the FFT for the correlation.
_____1: I = SWEEPLENGTH / SAMPRATE - 1
_____2: K = 1
_____3: FOR J = SWEEPLENGTH/SAMPRATE TO SWEEPLENGTH/SAMPRATE - 100 BY -1
_____4: Auto[K] = 0
_____5: FOR L = 1 TO J
_____6: Auto[K] = R[L] * R[I - J + L] + Auto[K]
_____7: NEXT L
_____8: K = K + 1
_____9: NEXT J

To see the wavelet shape, the correlation string is reversed and paired with itself:
_____10: FOR I = 1 to 100
_____11: Wvlt[I] = Auto[100 - I + 1]
_____12: NEXT I
_____13: FOR I = 101 TO 199
_____14: Wvlt[I] = Auto[I - 99]
_____15: NEXT I

The convolution process uses the phase and amplitude obtained earlier, so if you don't want to look at the sweep wavelet, the above steps are not necessary.

23.12 Convolution of a Wavelet with the Reflection Coefficients
With the wavelet and reflection coefficient strings both in the time domain, we are now able to create a synthetic seismogram. The process is called convolution. In the time domain, the following steps are needed:
1. reverse the order of the coefficients of the wavelet.
2. multiply each of the W values of the reversed wavelet with the first W values of the reflection coefficient string and sum the results; this is the first value of the output trace. (W is the number of samples in the wavelet).
3. shift the wavelet down the coefficient string by one sample and multiply/sum as in step 2; this is the second output sample.
4. repeat until the wavelet has been shifted to the end of the reflection coefficient string.
5. to prevent unwanted end effects, the reflection coefficient string must be padded with zeros at both ends before the convolution is started; the length of the pad is equal to the wavelet length.

NUMERICAL EXAMPLE:

Notice how the tail end of the first reflection interferes with the beginning of the second reflection at data point 6.

This technique, while simple to program, is computer intensive. The number of multiplications equals N squared, where N is the length of the padded reflection coefficient string.

A more efficient technique uses the Fast Fourier Transform, which converts the time series wavelet and reflection coefficient strings to their equivalent frequency domain amplitude and phase spectra. Convolution of the two involves the following:
1. calculate phase and amplitude of wavelet, using FFT algorithm.
2. calculate phase and amplitude of reflection coefficients, using FFT algorithm.
3. multiply the two amplitude spectra and add the two phase spectra.
4. convert these spectra to the time domain using the inverse FFT algorithm.

This method uses N * log2(N) multiplications instead of N squared. For example, a 5000 point transform takes 61,438 multiplications, while the direct approach takes 25,000,000, a speed improvement of over 600 times. Thus a computer which did the FFT in one minute would take 6.7 hours to compute longhand. For a million data points, the long method would take four weeks compared to a one minute FFT.

A sample program used for performing the FFT is shown below. As written, the program uses single precision numbers, and results may be a bit noisy. Real programs use double precision for the important array variables.

1 REM **********************************************************************
3 REM FAST FOURIER TRANSFORM PROGRAM - FORWARD
5 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
23 REM **** CAUTION: GWBASIC sets true = NEGATIVE!!!!!!!!!!!!!!!!
24 REM only on arithmetic operators - line 4070 needs fixes if true = POSITIVE
25 REM **** CAUTION: GWBASIC expects angles in radians
40 OPTION BASE 1 'array indexes start at 1
41 PI = 3.141593 'single precision
42 WVLTLENGTH = 256 'data points
43 SAMPRATE = 1 'milliseconds
44 SWEEPLENGTH = 1 'seconds
45 INITPHASE = 0 'radians
80 FFTLENGTH = WVLTLENGTH/2
83 L = LOG(FFTLENGTH) / LOG(2) 'number of cycles in FFT
84 M = FFTLENGTH 'abbreviated name for use in FFT
85 FREQINT = 1000 / (2^(L+1)) / SAMPRATE 'output frequency sample rate from FFT
86 DIM R [M], I[M], WVLT[WVLTLENGTH] 'Real and Imag arrays for FFT
4430 P1 = 0: IF FFTFLAG$ = "IFT" THEN P1 = 1: GOTO 4800 'Inverse Transform
4440 PRINT "FFT Loop 1": K = 0: FOR J = 1 TO M-1: I = 2 'FFT Starts Here
4455 IF K < M/I THEN GOTO 4470
4460 K = K - M / I: I = I + I: GOTO 4455
4470 K = K + M / I: IF K <= J THEN GOTO 4510
4480 A = R[J+1]: R[J+1] = R[K+1]: R[K+1] = A 'Binary Sort
4490 A = I[J+1]: I[J+1] = I[K+1]: I[K+1] = A
4510 NEXT J: G = .5: P = 1
4511 PRINT "FFT Loop 2": FOR I = 1 TO L: G = G + G: C = 1: E = 0
4520 Q = ((1-P)/2)^.5 * (1-2*P1)
4522 IF I = 1 THEN P = -1 ELSE P = ((1+P)/2)^.5
4525 FOR R = 1 TO G
4530 FOR J = R TO M STEP G+G: K = J + G
4540 A = C*R[K] + E*I[K]
4550 B = E*R[K] - C*I[K]
4560 R[K] = R[J] - A
4570 I[K] = I[J] + B
4580 R[J] = R[J] + A
4590 I[J] = I[J] - B
4600 NEXT J: A = E*P + C*Q: C = C*P - E*Q: E = A: NEXT R: NEXT I
4610 IF FFTFLAG$ = "IFT" THEN RETURN 'End of IFT
4801 REM **********************************************************************
4803 REM FAST FOURIER TRANSFORM PROGRAM - INVERSE
4805 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4810 PRINT "FFT Loop 3" 'FFT Continues
4820 A = 4 * ATN(1) / M: P = COS(A) 'Start IFT
4830 Q = (1 - 2*P1) * SIN(A): A = R[1]: R[1] = A+I[1]: I[1] = A-I[1]
4835 IF FFTFLAG$ = "IFT" THEN GOTO 4850
4840 R[1] = R[1]/2: I[1] = I[1]/2
4850 C = 1 - 2*P1: E = 0: FOR J = 2 TO M/2: A = E*P + C*Q: C = C*P - E*Q:
E = A: K = M - J + 2
4855 A = R[J] + R[K]
4860 B = (I[J] + I[K])*C - (R[J] - R[K])*E: U = I[J] - I[K]
4865 V = (I[J] + I[K])*E + (R[J] - R[K])*C: R[J] = (A+B)/2: I[J] = (U-V)/2
4870 R[K] = (A-B)/2: I[K] = -(U+V)/2: NEXT J
4875 I[M/2+1] = -I[M/2+1]: IF FFTFLAG$ = "IFT" THEN GOTO 4440 'IFT Continues^
4880 FOR I = 1 TO M: R[I] = R[I]/128: I[I] = I[I]/128: NEXT I
4885 RETURN 'End of FFT
4901 REM **********************************************************************
4903 REM RECREATE WAVELET FROM PHASE AND AMPLITUDE
4905 REM ********* Author: E.R. Crain, P.Eng. 29 Oct 1987 *******************
4907 PRINT "Recreating Wavelet"
4910 FOR I = 1 TO FFTLENGTH: R[I] = A[I]*COS (P[I]*PI/180): I[I] = P[I]*SIN(P[I]*PI/180): NEXT I
4915 FFTFLAG$ = "IFT": GOSUB 4400 'Do IFT on Wavelet
4920 PRINT "FFT Loop 4R": RECSCALE = 5000*SAMPRATE
4925 FOR I = 1 TO FFTLENGTH: WVLT[2*I-1] = -I[I]/RECSCALE: NEXT I
4926 FOR I = 1 TO FFTLENGTH: WVLT[2*I] = -R[I]/RECSCALE: NEXT I
4927 WVLT[1] = 0
4930 PRINT "Plotting Recreated Wavelet"
4935 VIEW (300,10)-(500,110)
4940 WINDOW (-20,-2)-(128,2)
4945 LINE (0,0)-(0,0),0
4950 FOR I = 1 TO LENGTH: LINE -((I-1) * SAMPRATE, WVLT[I]),7:NEXT I
5000 REM **********************************************************************
5002 REM END FOURIER TRANSFORM SUBROUTINES
5007 REM *********************************************************************

23.13 Case Histories: Synthetic Seismograms
Figures 23.19 through 23.22: Wabamun Salt
This example is taken from "Geophysical Aspects of Wabamun Salt Distribution in Southern Alberta", by N.L. Anderson, R.J. Brown, and R.C.Hinds, CSEG Journal, Dec 1988. Figure 23.19 shows sonic logs for two wells, one with Wabamun and Cairn salt beds, the other without. Both sonic logs have noise and spikes which should have been edited.

FIGURE 23.19: Wabamun example - raw logs

The synthetic in Figure 23.20 (top) shows a model with two salt layers based on sonic data only. Adding the density data increases the amplitude of the reflections (Figure 23.20, bottom). Figure 23.21 illustrates the salt free case. Comparison of the models make it relatively easy to itemize the characteristics of salt and salt free cases. Notice that the wavelet frequency is critical to identification of the formation tops.

FIGURE 23.20: Synthetic with salt layer - sonic only (top), density added (bottom)

FIGURE 23.21: Wabamun example - no salt case

Figure 23.22 portrays a synthetic seismic section based on interpolation of the sonic logs between the salt and salt free cases. The model would be improved by use of density data augmenting the sonic.

FIGURE 23.22: Wabamun example - synthetic seismic section

Figures 23.23 through 23.29: Swan Hills/Beaverhill Lake (D-3) Reef
This model shows the capability of a particular commercially available seismic modeling package (GMA). The system takes edited and modeled well log data as its basic input. The stratigraphic cross section is digitized, along with the log data, or is derived from well history data in the data base (Figure 23.23).

FIGURE 23.23: Correlated stratigraphic depth cross section

In this example, the reef was edited to a geologically believable shape (Figure 23.24). The acoustic impedance cross section is then computed and displayed (Figure 23.25). This cross section uses raw sonic log data that has been interpolated between well control. The sonic logs may need editing for bad hole, casing, rock alteration, or gas effects. If gas or carbonates are present, the model should incorporate density information as well.

FIGURE 23.24: Edited correlations and hypothetical reef

FIGURE 23.25: Acoustic impedance cross section

A synthetic seismic section is calculated by creating reflection coefficients from impedance and convolving a wavelet with each modeled trace (Figure 23.26).

FIGURE 23.26: Synthetic seismic section

This is the usual form of synthetic modeling, but this example used logs which received no editing or correction for invasion. If the reef contained gas or had higher porosity than the platform rock, the synthetic section would be wrong because the modeled velocity and density would be wrong.

Thus a number of different formation models could be created. The seismic response of the models can be compared to each other and to field data using discriminant analysis or frequency analysis techniques or, more commonly, by visual inspection. If sufficient discrimination between various models can be detected, then the comparison of these models with field data should provide a valuable interpretation aid.

Another form of modeling is ray path tracing (Figure 23.27), a time consuming computer solution. The ray paths to each key horizon are computed from the laws of refraction, and a synthetic section is prepared using the noise free reflection coefficients at each layer modeled. Both normal (perpendicular) and vertical incidence sections are displayed in this example (Figures 23.28 and 23.29). These are used to judge the quality of the reflections to be expected from the steeply dipping beds.

FIGURE 23.27: Ray path tracing on reef model

FIGURE 23.28: Synthetic seismic section from ray path tracing - normal incidence

FIGURE 23.29: Synthetic seismic section from ray path tracing - vertical incidence

Figures 23.30 through 23.33: Keg River (Granite Wash) Sand
This example is taken from "Seismic Perspective on Panny and Trout Fields of North Central Alberta", by N.L. Anderson, R.J.Brown, and R.C.Hinds, CSEG Journal, Dec 1989. Figures 23.30 and 23.31 show sonic log models for a thinning lower anhydrite contrasted with a thinning lower clastic. These models help interpreters determine which case might be more likely on a seismic section.

FIGURE 23.30: Keg River model - thinning anhydrite

FIGURE 23.31: Keg River example - thinning sand

A typical synthetic seismogram is given in Figure 23.32 as well as a geologic cross section from logs. The seismic model, synthetic section, and real seismic section are given in Figure 23.33. The basement high and surrounding basal sand are the key elements in this model. Note that there is no significant reflection from the bald high because the evaporites above mask the event. The model would be more realistic if density data were incorporated.

FIGURE 23.32: Keg River synthetic seismograms and log cross section

FIGURE 23.33: Keg River synthetic section and real seismic

Figures 23.34 through 23.38: More Sophisticated Models
Seismic modeling programs can generate the diffraction patterns from steep events, so that their interfering effects can be assessed. Figure 23.34 shows such an example. Here the diffractions might be mistaken for reflection on an anticline instead of the structural low that really exists.

FIGURE 23.34: Modeling refractions from steeply dipping beds

Stratigraphic traps, such as channels, beaches, and bars can also be modeled by using adequately modeled log data, as in Figure 23.35.

FIGURE 23.35: Stratigraphic trap model

23.14 In Conclusion
There are many connecting links between the seismic and well logging domains. Both develop velocity, density, and lithologic relationships from their measured data. A synthetic view of seismic response can be made from well logs, as can the inverse process create a filtered sonic log from seismic data. A clear understanding of the sources and definitions of acoustic velocity information, and the ability to communicate these differences, will go a long way toward integrating exploration and evaluation techniques.

23.15 Exercises for Chapter Twenty-Three
1. Assume an interval velocity of 10000 ft. per second and a single layer 5000 ft. thick. The total spread distance (X) is 4400 ft. Find average velocity and normal moveout. (10 marks)

2. Calculate stacking velocity for the above example. (10 marks)

3. Assume Vrms = Vstk and calculate average velocity and true depth from the table of data given below:

4. Under what conditions will these interval velocities compare with sonic log data? How does straight ray depth compare with true depth? (10 marks)

5. How does dip affect the above conclusions? (10 marks)

6. Why are synthetic seismograms needed and what major effects on the logs need to be accounted for? (10 marks)

7. List the steps required to make an accurate synthetic seismogram. (20 marks)

8. Explain the repercussions of using unedited or un-modeled data. (10 marks)

9. Choose one of the Case Histories in Section 23.13. Write a brief report covering the available input data, the steps taken to create the various outputs, the purpose for producing each output, and an interpretation of what these results show us. (20 marks)

23.16: Bibliography for Chapter Twenty-Three
SEISMIC TECHNIQUES

1: Procedure and technique in refraction profile interpretation; Source unknown; Course notes, 1968

2: Seismic energy sources 1968 handbook; United Geophysical Corp.; 38th Annual Meeting of Society of Exploration Geophysicists, 54 p., 1968

3: Automatic migration; Lindseth,R.O.; Course notes, 20 p., 1971