Theoretical and practical aspects of the Vibroseis method

(1)

Theoretical and practical aspects of the

Vibroseis method

(2)

Theoretical and practical aspects of the

Vibroseis method

(3)

ISBN 90-9002817-X Druk: ALEVO - Delft April 1989

Cover page: l.V.I. "birdwagen MARK IV off-road vehicle with 27 ton vibrator. (courtesy: I.V.I.)

(4)

Theoretical and practical aspects of the

Vibroseis method

Proefschrift

ter verkrijging van de graad van doctor aan de

Technische Universiteit Delft, op gezag van de Rector

Magnificus, prof. drs. P.A. Schenck, in het openbaar te

verdedigen ten overstaan van een commissie

aangewezen door het College van Dekanen op

donderdag 11 mei 1989 te 14.00 uur

door

Guido Jozef Maria Baeten

mijningenieur

(5)

Dit proefschrift is goedgekeurd door de promotor

Prof. A.M. Ziolkowski M.A., Ph.D., M.Sc. (Econ.)

(6)

The research reported in this thesis has been performed at the Delft University of Technology, Department of Mining and Petroleum Engineering, Section Applied Geophysics. Part of this research has been financially supported by Industrial Vehicles International, Inc., Tulsa, Oklahoma, U.S.A.

(7)

Notations and conventions

In this thesis, vectors and matrices are typed bold-faced. To locate a point in space, orthogonal Cartesian coordinates X\,x2 andx3 are employed with respect to a given orthogonal Cartesian

reference frame with origin O and the three mutually perpendicular unit vectors i}, i2 and i3; in the given order, ij, \2 and i3 form a right-handed system. The unit vector i3 points downwards.

The vector x = X\ ij + x2 i2 + *3 «3 denotes the three-dimensional position vector in the xhx2, x-i space.

Occasionally, the subscript notation for (Cartesian) vectors and tensors is used. Latin subscripts are to be assigned the values 1, 2 and 3, while Greek subscripts are to be assigned the values 1 and 2; for repeated subscripts, the summation convention holds.

Use is made of the temporal and spatial Fourier transform pairs. Quantities in the time domain are denoted by lowercase symbols, quantities in the frequency domain by the corresponding capitals, and quantities in the wavenumber domain are denoted by a circumflex (A) on top. In most cases, the parameters "co" and "/" are omitted for brevity. The temporal Fourier transform pair is defined as

f(t)=(2n) I F(ti) exp(-im) doo ,

1

F(a>) = I ƒ « exp(icor) dr .

The spatial Fourier transformation pair is given by -2 r°° r°°^

F(x}, x2, (0) = (In) I I F(kx, k2, a>)txp[i(k-ix-l+k2x2)]dk-[dk2 ,

F(k} , k2, a)- I I F(x: ,x2, co) c%p[-i(kix1 +k2x2)]dxidx:

(11)

List of major symbols

Latin symbols

03, ar accelerations of baseplate and reaction mass, respectively (ms- 2) cp compressional (P) wave velocity ( m s- 1)

cR Rayleigh wave velocity (ms- 1)

cs shear (S) wave velocity (ms- 1)

Cjjpq stiffness tensor ( N m- 2)

D flexural rigidity (Nm) etj strain tensor

j applied force delivered by the drive mechanism of the vibrator (N) Gfcx Green's function for flexural rigidity method (s2kg"')

Hn , Hn Bessel functions of the third kind (also called Hankel functions), of

order n

/„, Kn Modified Bessel functions of order n

ï j , ï2> b Cartesian unit vectors

./„ Bessel function of the first kind of order n

Mp, Mr, Mh masses of baseplate, reaction mass and holddown mass, respectively

(kg)

kp, ks, kR P - S - and Rayleigh wavenumber (m_ 1)

t time coordinate (s)

f,- traction vector (Nm- 2)

t;j stress tensor ( N m- 2)

Hi particle displacement vector (m)

»3, ur, uh displacements of baseplate, reaction mass and holddown mass,

respectively (m)

v (- particle velocity vector (ms~]) x position vector

.V], x2: v3 Cartesian coordinates

Yn Bessel function of the second kind (also called Weber's function) of

(12)

Greek symbols

A Rayleigh denominator (s4m~4)

A , \x Lamé elastic parameters ( N n r2)

7 , ys vertical P - and S wave slownesses ( s n r ' )

p volume density of mass (kgirr3)

o surface density of mass of the baseplate (kgni-2)

v Poisson's ratio

(13)

Chapter 1 An introduction to Vibroseis

In seismic exploration, the use of a vibrator as a seismic source has become widespread ever since its introduction as a commercial technique in 1961. This chapter explains the principles on which the Vibroseis™ method is based. The mechanism which allows the seismic vibrator to exert a pressure on the earth is explained. The two basic features of the force function generated by the seismic vibrator are discussed: first, the signal generated by a seismic vibrator is not impulsive, but has a duration of several seconds; second, the monitoring system of the seismic vibrator allows control over the outgoing vibrator signal. The choice of the signal that is to be monitored on the vibrator in order to control the outgoing signal is a controversial issue. A brief review of existing methods is given. Finally, the advantages and disadvantages of Vibroseis over most impulsive sources are discussed, and the objectives of the present investigations are explained.

1.1 The mechanics of the seismic vibrator

An example of a seismic vibrator is shown in Figure 1.1.

(14)

The vibrator is a surface source, and emits seismic waves by vibrating a plate which is kept in contact with the earth.

Figure 1.2 The force-generating mechanism of the seismic vibrator source.

The driving force applied to the plate is supplied either by a hydraulic system, which is the most common system in use, or an electrodynamic system, or by magnetic levitation, the latter being a new development in the field of seismic vibrator technology. The direction in which the plate vibrates can also vary: P wave vibrators (where the motion of the plate is in the vertical direction) as well as S wave vibrators (vibrating in the horizontal direction) are used. Finally, a marine version of the seismic vibrator has been developed.

For all these vibrator types, the general principle which governs the generation of the driving force applied to the plate (usually referred to as the baseplate) can be described by the configuration shown in Figure 1.2. A force F is generated by a hydraulic, electrodynamic or magnetic-levitation system. A reaction mass supplies the system with the reaction force necessary to apply a force on the ground. The means by which this force is actually generated is illustrated in Figure 1.3, in which the principle of the hydraulic drive method is shown. By pumping oil alternately into the lower and upper chamber of the piston, the baseplate is moved up and down. The fluid flow is controlled by a servo valve. The driving force acting upon the baseplate is equal and opposite to the force acting upon the reaction mass, as can easily be

(15)

Reaction mass

Piston

/Fluid pressure

from

servo valve

Baseplate

Figure 13 Schematic view of the generation of the driving force for a hydraulic vibrator. inferred from Figure 1.3. In general, the peak force is such that the accelerations are in the order of several g's, so that an additional weight has to be applied to keep the baseplate in contact with the ground.

For the hydraulic and electrodynamic vibrators, the weight of the truck is used for this purpose. This weight, commonly referred to as the holddown mass, is vibrationally isolated from the system shown in Figure 1.3 by an air spring system with a low spring rate (shown in Figure 1.4), and its influence on the actual output of the system is usually neglected. The resonance frequency of the holddown mass is in the order of 2 Hz, the lowest frequency of operation in Vibroseis seismic surveys for exploration purposes being usually not less than 5 Hz.

1.2 The shape of the signal emitted by a seismic vibrator

The signal emitted by the seismic vibrator is not impulsive, but typically has a duration of 10 sec. This seems to be in contradiction with the fact that seismic exploration methods aim to detect the impulse response of the earth. This apparent contradiction can be clarified by taking a closer look at the properties of an impulse and its earth response.

A perfect impulse at time f=0 contains all frequency components with equal amplitude and zero phase. This is illustrated in Figure 1.5.

(16)

STILTS VIBRATOR ACTUATOR

AIR SPRINGS

(a) (ft)

Figure 1.4 (a) schematic view of the Vibroseis truck with the air springs, the baseplate and the vibrator actuator (reaction mass), and (b) detailed view of the middle pan of the truck.

'lime t = 0

(fl)

Amplitude Phase

-► frequency frequency

TIME DOMAIN FREQUENCY DOMAIN

(b)

Figure 1.5 The notion of a perfect impulse, (a) in the time domain, and (b) its corresponding frequency domain

(17)

In practice, one cannot generate a perfect impulse because this would require an infinite amount of energy; the best one can achieve is to emit a bandlimited impulse, resulting in a finite-amplitude wavelet whose time duration is small compared with any time interval of interest.

The Vibroseis source emits a bandlimited, expanded impulse; The band limitation has two aspects: at the low frequency end, it is dictated by the mechanical limitations of the system (the amount of travel or "stroke" of the reaction mass, and the amount of hydraulic flow that can be supplied by the pump) and the size of the baseplate. The high frequency limit is determined by the mass and stiffness of the baseplate, the compliance of.the trapped oil volume for a hydraulic vibrator and mechanical limitations of the drive system.

The notion of an "expanded" impulse can be explained in terms of the amount of energy per unit time, further referred to as energy density. In an impulsive signal, all energy is concentrated in a very short time period, leading to a very high energy density. In the Vibroseis method, the same amount of energy can be transmitted over a longer time, so that the energy density of the signal is reduced considerably. This reduction in energy density can be achieved by delaying each frequency component with a different time delay, while keeping the total energy contained in the signal constant. Thus, instead of emitting a signal with a flat amplitude spectrum and a zero phase spectrum, a signal is created which has the same flat amplitude spectrum in the bandwidth of interest, but has a non-zero phase spectrum. The frequency-dependent phase shifts cause time delays which enlarge the duration of the signal. However, the total energy of the signal is determined only by its amplitude spectrum. Since the amplitude spectrum of the expanded impulse is equal to the amplitude spectrum of the true, non-expanded impulse, in this way a signal is created which has the same amount of energy as is contained in the zero phase impulse.

The effect which the increased time duration of the emission has on the recorded seismogram needs to be eliminated. This can be achieved by controlling the phase function of the emitted signal. Then, the signal received at the geophone can be corrected for the non-zero phase spectrum of the source wavelet by performing a cross-correlation process. To clarify this point, let the source wavelet be denoted by s(t). If the convolutional model is adopted to describe the response at the geophone, x(t), the following expression is obtained in the absence of noise:

*(») = s(0 * S(0 , (1-D

where g{t) denotes the impulse response of the earth, and "*" denotes a convolution. The validity of the convolutional model for the Vibroseis configuration is discussed in Chapter 2. Transforming equation (1.1) to the frequency domain yields

(18)

If the received signal x(t) is cross-correlated with the source signal s(t), the signal c(t) is obtained which, in the frequency domain, is given by

2

C(o))=X(fl))S*(©) =|S(©)| G(co) , (i.3) since cross-correlation of x(t) with s(t) in the time domain corresponds to a multiplication in the

frequency domain of X(a) with the complex conjugate of S(co). In this equation, the complex conjugate is denoted by the superscript "*". Since the amplitude spectrum of S(co) is flat over the frequency band, and zero outside this frequency band, it follows that by cross-correlating the measured seismogram x(t) with the source function s(t) the (scaled) bandlimited impulse response of the earth is obtained. The scaling factor is simply the (constant) amplitude spectrum of the source signal squared. In the time domain, the cross-correlated seismogram is a convolution of the earth's impulse response with the autocorrelation of the sweep. In Vibroseis applications, the measured seismogram x(t) is ususally referred to as "vibrogram", whereas the seismogram after cross-correlation, c(t), is denoted as "eorrelogram".

Therefore it is possible to emit a signal with a flat amplitude spectrum and an arbitrary phase spectrum, and still obtain the exact (bandlimited) impulse response of the earth, provided the phase spectrum of this source wavelet is known. It may be observed that this cross-correlation is merely a special deconvolution process, in which we exploit the fact that the amplitude spectrum is constant, thus avoiding the design of deconvolution filters.

This completes the general framework into which the vibrator's source signal is placed. The detailed shape chosen for the vibrator signal in practice is now described. A signal q(t) (usually referred to as the "sweep") is generated by the sweep generator of the seismic vibrator. The signal s(t) introduced in equation (1.1) is different from the sweep q(t), for reasons that are explained when the control óver the vibrator source is discussed in Section 1.3. The generated sweep q(t) has the general form

q(i) = a{t)sin [ 2 * 0 ( r ) ] . (1.4) a(t) is a taper function which is usually chosen to be a linear or cosine roll-off taper, and

typically has a length of 250 msec at each end of the sweep. The taper is applied only at the beginning and at the end of the sweep to reduce Gibbs' oscillations. The sweep typically has a duration T of 10 sec. The function 0(0 determines the frequency of the sine wave as a function of time; more specifically, the derivative of 0(0 is equal to the instantaneous frequency ƒ ""'(0 of the emitted sine wave:

(19)

dm

dl

= e-(o =/"'(o .

(1.5)

As an example, a sweep is considered where the instantaneous frequency varies linearly with time, as illustrated in Figure 1.6:

r\

t)

=

h+

hzf±

t

,

(1.6)

in which f0 is the starting frequency of the sweep, ƒ, its end frequency and T the sweep

duration.

Figure 1.6 The instantaneous frequency of a linear sweep as a function of time.

This sweep with a linearly-varying instantaneous frequency is commonly referred to as a linear sweep, and is given by

qit) =a(f)sin[2ff(/o +jf-±j/°t )t) (1.7)

If the first sweep frequency fa is smaller than the last frequency ƒ,, the sweep is called an

upsweep; if/0 is larger than/,, it is called a downsweep.

Although the relation between the sweep and its instantaneous frequency as a function of time is readily established, the Fourier transform of the sweep cannot be detennined analytically. At this point, a clear distinction must be made between the concept of "instantaneous frequency", which is the single frequency contained in the signal at a single time instant, and the concept of

(20)

a frequency component as used in conjunction with a Fourier-transformed time function. In the latter case, the spectrum of q(t) can in general not be obtained analytically; however, Rietsch (1977a, 1977b) has shown using the method of stationary phase that the spectrum Q( ƒ) of q{t) can be approximated by

QV) = A(f) e x p [ - 2 a i ( 0 ( T ) - / T ) + i J ] , (1.8)

in which the different sign convention for the temporal Fourier transform adopted by Rietsch and the different sweep definitions (Rietsch uses a cosine function instead of a sine function to define the sweep) are taken into account. In this equation, T is the stationary point, given by the solution of the equation

0 •(*)-ƒ=() , (1.9)

and/if/) is given by

*(/>=

a(T)

,

7

- o.io)

2 [ 0 " ( T ) ] 2

In equation (1.8), an upsweep is assumed; this, however, is not essential. The meaning of this equation can be illustrated using the previous example of the linear upsweep. For this signal, the stationary point T is given by

*=——lf-fo] . (1.11) / l 'JO

and it can easily be verified that the phase spectrum of the linear upsweep is a second degree polynomial in/. Since the second derivative of 6(t),

e r

(

T > - ^ . a.n)

is a constant, the amplitude spectrum of the sweep is a scaled version of the taper function

A(f), which is flat over most of the bandwidth. Sometimes a nonlinear sweep is taken as a

source signal, for instance to account for frequency-dependent absorption (Hargrove et al, 1984).

(21)

50 100 Frequency (Hz) (b) S 150000 -50 100 Frequency (Hz) (c)

Figure 1.7 An 8 sec, 10-100 Hz upsweep with a taper length of 250 msec, (a) the sweep in the time domain; the frequency range for this Figure is 1-5 Hz for display purposes, (b) the amplitude spectrum of the sweep, (c) the phase spectrum of the sweep, in degrees, and (d) the autocorrelation of the sweep.

In Figure 1.7 the concepts outlined above are illustrated for the example of a linear upsweep. An 8 sec, 10-100 Hz linear upsweep is used with a taper length of 250 msec. Figure 1.7a shows the sweep q(t). Because the oscillations in the sweep are too rapid to yield a clear picture, the frequency limits for this figure are 1-5 Hz. Figures Mb and 1.7c show the amplitude and phase, respectively, of Q(f). It can be observed from these figures that the phase indeed is a quadratic function of frequency, and that the amplitude spectrum of the sweep is constant over the bandwidth, apart from some Gibbs' oscillations at the beginning and end frequencies. Finally, Figure Md shows the autocorrelation of the sweep.

(22)

1.3 The control over the Vibroseis source signal 1.3.1 Engineering aspects

The Vibroseis source is designed in such a way that the operator has, at least in principle, perfect control over the emitted wavelet. It follows from the previous discussion that knowledge of the source wavelet emitted by the vibrator is essential. The Vibroseis method presents the opportunity for this kind of signal control. This is achieved by the system shown in Figure 1.8.

feedback signal

f(0

compensator inputforce modification

vibrator

Pre-determined

sweep input q(t)

Figure 1.8. The feedback system of a seismic vibrator.

The pre-determined sweep input q(t) (also called the pilot sweep) is compared with a feedback signal f(t), which is a signal measured on the vibrator. In principle, any vibrator motion that can be measured can serve as a feedback signal. The main function of the feedback system is to produce a feedback signal which resembles the pre-determined sweep as close as possible. Amplitude and phase differences between these two signals are corrected for by a change in the time-varying input force. This correction for differences between the feedback signal and the pre-determined sweep is necessary because the vibrator electronics can, for the example of a

(23)

hydraulic vibrator, only control the oil flow in the hydraulic piston, which in general is not equal to the feedback signal chosen on the vibrator.

Until recently, the control system of the vibrator was designed to correct only for phase differences between die pre-determined sweep and the vibrator feedback signal. Since the early 1980's, amplitude control techniques have been in use which allow the vibrator to operate as close to the maximum holddown weight as possible without decoupling.

1.32 Geophysical aspects

In the previous section, the control of the vibrator source by means of a feedback system was discussed. It was shown that the problem of ensuring that the pre-determined sweep q(t) and the feedback signal f(t) are equal is an engineering problem.

The geophysical problem of the vibrator is the determination of the far field wavelet from measurements on the vibrator. For conventional geophysical applications, the wavelet s(t) that appears in the convolutional model is of interest. In the convolutional model, it is assumed that the wavelet s(t) does not change its shape while propagating, and is the same in all directions. The first demand, that 5(0 remains constant, implies that s(t) is the far field wavelet; the second demand, the source being non-directional, implies that the source dimensions must be small compared with a wavelength. The latter condition is usually satisfied when a single vibrator is considered. It should be noted that the term "far field wavelet" usually refers to the far field

particle velocity, since in normal land seismic applications geophones are used for signal

detection.

We intend to show that it is theoretically possible to calculate the far field wavelet from measurements on the vibrator. This calculated wavelet is called svib(t). The band-limited

impulse response of the earth is obtained by cross-correlating the received seismogram x(t) with this far field vibrator signal svih(t). The frequency domain expression for the seismogram after

cross-correlation, c(t), is thus given by

C{co)=X(co)Svib*(.co) = 5(w) Svib' G(co) . (1.13)

If s(t) and svib(t) are the same, the desired zero-phase impulse response g(t) is obtained.

Whether they are the same or not depends on the validity of the theory.

1.3.3 Combined engineering and geophysical aspects

The two aspects of vibrator control that are discussed in the two previous sections can be summarized as follows:

(24)

1. Engineering problem: make the pre-determined sweep q(t) and the feedback signal f(t) equal. 2. Geophysical problem: determine svib(t) from measurements on the vibrator such that svib(t)

is equal to the far field wavelet s(t).

In all present applications of the Vibroseis technique, the feedback signal f(i) is different from the far field vibrator signal svib(t). One may wonder why it is not normal practice to make f(t)

equal to svib(t). There are several reasons for this discrepancy between feedback signal and far

field vibrator signal:

1. There is still discussion about the actual shape of the far field wavelet s(t), and, consequently, about the choice of the far field vibrator signal svib(t).

2. The far field vibrator signal svib(t) can only be calculated approximately from measurements.

3. Often, only the phase information of the feedback signal is considered for source control purposes.

4. Even if the far field vibrator signal is known, and can be determined exactly from measurements on the vibrator, practical reasons may prevent the feedback signal f(t) from being chosen equal to the far field vibrator signal svlb(t). For example, it was stated in the

previous section that vibrators are usually operated as close to the maximum holddown-weight as possible without decoupling, to maximize the energy output. This condition need not necessarily have any consequences for the control of the far field wavelet shape. The possible discrepancy that may exist for practical reasons between the feedback signal f(i) and the far field vibrator signal svib(t) (point 4 above) is not essential; in the remainder of this

thesis it is assumed that the feedback signal on the vibrator is chosen such that it attempts to predict the far field vibrator signal, and thus the far field wavelet.

In practice, the measured seismogram x(t) is cross-correlated in real time with the pre determined sweep q(t). Only the resulting correlogram c(t), in the frequency domain given by

C(a>) = X(co) Q\co) , (1.14)

is recorded. It is thus assumed in practice that the feedback system of the vibrator ensures that the feedback signal f(t), which is a measure of the actual vibrator motion, and the pre determined sweep q(i), which is the output of a sweep generator, are equal. One of the reasons for this assumption is that the feedback signal contains the effect of harmonic emissions due to the non-linear behaviour of the vibrator hydraulics.

(25)

The question of which feedback signal should be used to predict the far field wavelet has been a controversial issue since it was first raised by Lerwill in 1981. To establish the correct feedback signal, two questions have to be addressed:

1. What is the relation between the shapes of the near field and the far field wavelet?

2. Given this relation, how can the far field wavelet shape be related to measurements on the vibrator?

So far, three different feedback signals have been in use in Vibroseis surveys: baseplate acceleration, reaction mass acceleration, and a weighted sum of baseplate acceleration and reaction mass acceleration. All three methods produce different answers to the two questions stated above.

Using baseplate acceleration as a feedback signal, it is assumed that

1. The far field particle velocity is in phase with the near field particle velocity. 2. The near field particle velocity is equal to the baseplate velocity.

Using the reaction mass acceleration as a feedback signal, it is assumed that 1. The far field pressure or panicle velocity is in phase with the near field pressure. 2. The near field pressure is proportional to the reaction mass acceleration.

Using a weighted sum of baseplate acceleration and reaction mass acceleration as a feedback signal, it is assumed that

1. The far field displacement is in phase with the integrated surface traction.

2. The integrated surface traction is equal to a weighted sum of baseplate and reaction mass acceleration.

1.4 The use of Vibroseis in seismic exploration

One may wonder why it is not normal practice in seismic exploration to use an impulsive source, since, after all, it is the earth's impulse response we are after. As can be seen in Figure 1.9, the most well-known impulsive seismic source, dynamite, is indeed used very often in land seismic surveys. There are, however, some distinct disadvantages related to the use of an impulsive source like dynamite.

(26)

1962 1967 1972 1977 1982 1987 Year

&0Ö3 Vibroseis I I dynamite

Figure 1.9 The contribution of Vibroseis and dynamite to the total number of crew months spent in land petroleum exploration, in %, for the years 1962-1987*. First of all, due to the high energy density of the dynamite explosion, severe harm can be done to the environment. In any case, the destructive nature of the dynamite source prohibits its use in densely populated areas. Second, a hole has to be drilled for every shotpoint in which the dynamite charge is placed. Third, the high energy-density of the explosion results in a non linear zone surrounding the explosion. Although the ignition time of the dynamite itself is short compared with any time duration of interest in seismic exploration, this nonlinear zone results in a distorted wavelet. It has been shown by Ziolkowski and Lerwill (1979) that the high frequency content of the signal decreases when the charge size is increased (the low frequency content increases). This yields a trade-off between penetration and resolution: a large charge size has better penetration, but lacks high frequencies. Another disadvantage of the creation of a nonlinear zone around the dynamite explosion is that effectively a wavelet is transmitted into the earth that is not an impulse, and that has a shape which is unknown and cannot be measured directly.

Foi 1962-1971, the figures for Vibroseis must be interpreted as maximum numbers, since only a division between dynamite and non-dynamite was made for these years. Also, marine surveys were included in the 1962-1973 figures, but the contribution of marine surveys never exceeded 10 % in these years. No figures were available for 1969. Data arc from the Geophysical Aclivity Reports in Geophysics (1963-1981) and Geophysics: The Leading Edge of Exploration (1982-1988).

(27)

The Vibroseis source has some distinct advantages over the dynamite source. First, the emitted signal contains an amount of energy that is (roughly) comparable to the energy contained in a dynamite signal: Janak's (1982) results from a land seismic source study indicate that a 2.5 kg dynamite charge buried at 5 m depth delivers an amount of energy of 9.5 x 106 Joules, whereas a single Litton model 311 P wave vibrator emitting a 10-79 Hz sweep (the drive level is not mentioned in the report) delivers 1.2 x 106 Joules. Because of the use of an expanded impulse, the energy density of the source wavelet in the Vibroseis technique is much less than the energy density of the dynamite wavelet. Therefore, destructive effects are much less severe. Second, Vibroseis provides us with a direct means to measure and control the outgoing wavelet. Third, there is no need to drill holes when using Vibroseis.

There are, however, also some disadvantages connected with the use of Vibroseis as a source. First, a single vibrator in general does not deliver a sufficient amount of energy required for seismic exploration purposes, so that arrays of vibrators have to be used. Typically, 4 vibrators vibrate at each vibration point simultaneously (Figure 1.10).

Figure 1.10 A seismic survey using an array of 4 vibrators.

Second, as vibrators are surface sources, large amounts of Rayleigh waves are generated. The generation of Rayleigh waves can be suppressed in a dynamite survey by placing the charge at the bottom of the weathered layer. In Vibroseis surveys, the Rayleigh waves have a very high amplitude and are an undesired feature on the seismogram. Third, the Vibroseis method can be

(28)

employed only in areas which are accessible to the seismic vibrator trucks, whose weight may exceed 20 tons. Fourth, correlation noise limits the ratio between the largest and smallest detectable reflections.

In spite of all its disadvantages, the Vibroseis method is now a standard method in the seismic exploration for hydrocarbons. In 1987, Vibroseis was used more often in land seismics than dynamite (the contribution of Vibroseis to the total number of crew months spent in land seismic petroleum exploration was 49 %, whereas the contribution of dynamite was 48.3 %). The operational advantages of the Vibroseis method over the conventional dynamite survey result in an average cost per kilometre of Vibroseis which is only two-thirds of the cost per kilometre for a dynamite survey (figures for 1987). Also, the average number of kilometres that can be covered per crew month is 30 % higher for Vibroseis surveys than it is for dynamite surveys (figures for 1987). This cost-effectiveness and efficiency, together with the increasing importance of signal control in the search for higher resolution of seismic data and the non destructive character of the method explains the increasing popularity of Vibroseis.

7.5 Objectives of the research in this thesis

In this section, the objectives of the research described in this thesis are explained.

The cross-correlation of the received seismogram with the sweep (which can be either the pre determined sweep input q(t) or the feedback signal/(/J) does not yield a seismogram in which the earth's impulse response is convolved with a zero phase wavelet. The reason is that the sweep and the far field wavelet s(t) are different. This thesis is devoted to

(a) explaining why this is so, (b) what can be done about it, (c) tests on real data.

The contents of this thesis is now discussed in more detail.

The first issue that is investigated is the propagation of the wavefield emitted by the seismic vibrator into the earth. Two questions are of particular interest:

1. For forward modelling purposes, how can the wavefield at any point in the earth be described, given the vibrator motion?

2. What vibrator motion should the received seismogram be deconvolved with in order to obtain the desired impulse response of the earth?

(29)

These two closely related questions are investigated in Chapter 2, where expressions are derived for the wavefield in a plane-layered elastic medium, and for the wavefield in the far field of an elastic halfspace. It is shown that the far field displacement is essentially equal to the reaction force exerted by the earth on the baseplate, or, in other words, to the ground force.

Once this choice of the far field vibrator signal has been established, it may seem straightforward to define the appropriate feedback signal in terms of measurements that can be made on the seismic vibrator. However, the complexity of the vibrator behaviour prohibits a simple connection to be made between the far field vibrator signal (the ground force, or, more precisely, the time derivative of the ground force), and vibrator signals that can be measured in practice (the accelerations of the baseplate and the reaction mass). Measurements show that the traction directly underneath the baseplate is not distributed uniformly over the plate. The ground force is equal to the integration of the traction over the entire baseplate area. Since the ground force has to be determined from measurements on the vibrator, a thorough understanding of the nature of this non-uniform traction distribution is required, as well as a good model for the relation between the traction underneath the baseplate on the one hand, and the accelerations of baseplate and reaction mass on the other hand. First, a mechanical model of the seismic vibrator is developed, which is described in Chapter 3. Since the motion of the vibrator is a result of the interaction between the vibrator and the earth, the coupled earth-vibrator model is investigated in Chapter 4. In this chapter, the distributions of traction and displacement on the baseplate are investigated using different baseplate models. Also, the power output of a vibrator and the radiation impedance of the baseplate is discussed. Chapter 5 addresses the question of which vibrator motion yields the best estimate for the ground force. In this chapter, theoretical results are compared with measurements made on a seismic vibrator, and in a borehole containing a downhole geophone array.

The question of which feedback signal must be chosen on the vibrator is also important in the marine environment. Since the marine vibrator differs in many aspects from the land vibrator, a special chapter (Chapter 6) is devoted to the marine Vibroseis configuration.

(30)

Chapter 2 The wavefield of a surface source in a horizontally

layered elastic medium

2.1 Introduction

A seismic vibrator is a surface seismic source which transmits waves into the earth by applying a force at the surface of the earth. In this chapter, the wavefield generated by an arbitrary distribution of tractions on the surface of an horizontally layered elastic earth is derived. It is not necessary to specify any characteristics of the vibrator; the theory can be applied to P and S wave vibrators, and to single vibrators as well as to arrays. In Chapter 3, the origin and nature of these surface tractions are specified.

Expressions are derived for the components of particle velocity and stress within or at the surface of a horizontally layered, elastic medium on top of which the surface source is acting. Each medium is assumed to be homogeneous, isotropic and perfectly elastic. The source is described by the distribution of traction components it generates at the surface, hereafter referred to as the surface traction components.

The well-known reflectivity method is used to derive the expressions for the wavefield in the medium. In this method, the propagation of the wavefield through the different layers is described mathematically by means of reflection matrices, containing the reflectivity properties of the individual layers as well as the propagation effects within each layer. The reflection matrix has to be calculated for each layer separately; this can be achieved by means of a recursive scheme. The reflectivity method is modified to include the surface tractions, resulting in an expression for the components of stress and displacement at any point in the medium in terms of the surface tractions. All derivations are made in the frequency - wavenumber domain.

The relations for an elastic half space are obtained from the expressions derived for the layered medium. These lead to expressions for the shape of the far field wavelet.

2.2 Configuration

A stack of A' plane elastic homogeneous isotropic layers is sandwiched between a lower halfspace and an upper surface on which the surface traction components are acting. Each layer

j (J=\,...., N) has a thickness /( and Lamé elastic parameters A ■, u ■ and density p associated

with it. The bottom of layer y' coincides with depth level .r3=.v, . The configuration is shown in Figure 2.1.

(31)

Bull Su hl rfare / Surface TraciionComponents / / Layer 1 A, . /i, . p , *l .1 Layer 2 A2, / i2, p2 h2 Layer N lK.fs,pN *v Layer A'+ 1 AA. , . pM , p ^ . , L w f r halfspace — ' S . ! —"3.2 —xy A-i - » 3 . A '

Figure 2.1 The eanh model, consisting of/V plane elastic homogeneous isotropic layers, bounded by a halfspace at the bottom and a surface at the top. On the surface, a distribution of tractions is present, as illustrated for points A and B.

2.3 Basic relations

The integral equations that relate surface traction components and the components of stress and particle velocity in each layer are derived from the equations that govern the wave motion. The relevant layer indices are, for the moment, omitted for simplicity; they are attached to the resulting equations in Section 2.5 when the coupling of the wave fields at the layer boundaries is included.

In the absence of body forces, the linearized equation of motion in the frequency domain is given bv

(32)

in which 7^ denotes the stress tensor of rank two, p the mass density of the medium, co the angular frequency and Ui the displacement vector. The components of stress and strain are

connected by means of the constitutive relation. For a perfectly elastic solid, the constitutive relation is

* ij — Cijpq kpq < (2.2)

where c^ is the stiffness tensor, and E is the strain tensor. For isotropic, homogeneous them

be written as

media, the number of independent elements in the stiffness tensor reduces to two, and c- can

CÜP? = A 5ij 5Pq + »<>Sip 5k + 8iq SjP> • (2-3)

In the absence of volume sources, the strain tensor E is given by

E„ = ±(BpUq+dfJp). (2.4)

To solve these equations, the procedure is followed as outlined by Kennett (1974, 1979, 1983), Ursin (1983) and du Cloux (1986). First, the linear wave equations are subjected to a two-dimensional Fourier transformation with respect to the space coordinates x1 and x2- Then,

the stress components which are not continuous across a layer interface (7"n , T12. ^21 anc*

r2 2) are eliminated. This yields the relation

33F = -ico E F , (2.5)

where F is a (6x1) column vector containing the field components:

F = [ F1, F2]T, (2.6)

F ^ l V j , ? , , ^ ]7, (2.7)

F2= r 7i3 , P1, V2]T, (2.8)

in which the superscript "T" denotes the transpose, E is the system matrix, \7,are the components of the particle velocity and 7",- are the components of the traction:

(33)

Ti = Tij rij = Ta (2.9)

The system matrix E has the following structure:

E = _(2.10)

in which the submatrices E, and E2 are given by

E , = A+2/i

k>\

Xp2 A+2/i Ap, X+2/u Ap2 X+2n 4/j(A+/i)p] 2 p W92 A+2jz H(3X+2n)p]p2 P X+2/i H(3X+2n)plp2 A+2/J 4p(X+n)p2 2 W i X+2p. X+2n (2.11) p P i Pi P\ 1 0 P2 0 1 (2.12) where Pa=-co (2.13)

2.4 Decomposition into upgoing and downgoing waveficlds

To solve equation (2.5), the wavefield is decomposed into upgoing and downgoing waves by diagonalizing the system matrix E. The diagonalization can be achieved by calculating the eigenvalues of the system matrix, and subsequently determining an independent set of eigenvectors. The decomposition is given by the equation

(34)

E = P Q£P - _(2.14)

In equation (2.14), Q£ is a diagonal matrix with the following structure:

Q£=

Q o o - Q

Q is a diagonal matrix containing the eigenvalues of the matrix E: Q = diag^, y^y,), with rP.s=\—-p % ( R e , l m ) yp s> 0 (2.15) (2.16) (2.17) and 2 2 2 P =P\+P2- (2.18)

cp and cs denote the velocity of the compressional and shear waves, respectively, and are given

by

\X + 2H %

. c, =

-%

(2.19)

The composition matrix P contains the eigenvectors of the system matrix. There is some freedom in the scaling of these eigenvectors, but the choice of a set of independent eigenvectors is by no means essential. Ursin's (1983) scaling is adopted because of the resulting simplicity of the inverse of the composition matrix. Ursin obtains

P =

4ï

Pi

P2 - P i

(2.20)

(35)

p , =

M

,PI

Mï)

:A

&

-2p

2

nm

and P , =

PM

2

'Ah

P) I r / 2 | J _ :

JP-W

Pi

V*

2

\c* I

p{pyf

A

\^-

2

P

-pApyf

A

-pi{pyf

/2 m

1

2W —

/ \A

\PI

*f

/ \A

AP)

%

P

iM

n

p\i \A P

4M'

%

-

P

Mn

M

V

\-

X

A

(2.21) (2.22)

Ursin (1983) derived similar expressions for the components of the eigenvector matrices, but his corresponding expression for Pj contains an error in the second and third element of the

- j - 2 p I is l

As a result of this choice of the eigenvectors, the following symmetry property is obtained: second column; the term| 2p is missing

, - i 1

Pi

(2.23)

(36)

F = P \V . (2.24)

Substituting equation (2.24) into equation (2.5), and using that from equation (2.14) it follows that

Q£= P E P , (2.25)

the relation

a3W = -iü> Q£W (2.26)

is obtained and, since Q£ is a diagonal matrix, the desired decomposition has been achieved. The vector W can be expressed as

(2.27) ^-^

w =

UP D where UP = [£//>,, UP2, UPj] (2.28)

is a vector containing the three components of the upgoing wavefield, and

D = [D1,D2,53]T (2.29)

is a vector containing the three components of the downgoing wavefield. Ü P and D do not correspond to physical quantities such as displacement and traction. Rather, they represent the upgoing and downgoing wavefield potentials, which makes a physical interpretation of these components difficult. When applying the boundary conditions at layer interfaces in a later section, these potentials must be transformed into the physical quantities panicle velocity and traction by applying the composition matrix P.

However, the introduction of these wave potentials leads to a simple solution of equation (2.26), given by

(37)

and

^ ^DOWN

D = A e\p(icoQx3). (2.31)

^ UP ^ DOWN

Here, A and A are vectors containing the amplitude factors for the upgoing and downgoing wavefields, respectively, and Q is the diagonal matrix containing the vertical P and S wavenumbers (see equation (2.16)).

The wavefield in the layer can now be constructed from these upgoing and downgoing waves by applying the composition matrix P (equation (2.24)):

F = P W , (2.32) ^ ^ UP DOWN

where the vector \V contains the unknown amplitude factors A and A . These amplitude factors can be determined using the boundary conditions at the interfaces of each layer. In the next section it is shown that application of the boundary conditions at the layer interfaces relates the upgoing and downgoing waves at the layer boundaries by means of a scattering matrix containing local reflection and transmission matrices.

2.5 The local reflection of the wavefield at a layer interface

In the previous section, relations for the wavefield in each individual layer were derived using the procedure of Kennett. As the coupling of these wavefields across layer boundaries is now discussed, the relevant subscripts are attached to the field quantities and medium parameters to indicate to which layer they refer. This notation does not interfere with the use of the summation convention: in the remainder of this chapter, the summation convention is only used in conjunction with Greek subscripts, whereas the layer indices are denoted by Latin subscripts.

Assuming firm contact at each interface, the components of particle velocity and traction are continuous across each interface. Since the components of particle velocity and traction are contained in the vector F , this boundary condition can be expressed as

(38)

in which the index "bottom" indicates that the relation holds at the bottom of the relevant layer, and the index "top" refers to the top of the relevant layer. Using equation (2.24), the formulation of this boundary condition in terms of upgoing and downgoing waves is obtained:

PjWj (bottom) = P;+1W;+1(t0p) . (2.34)

Using this equation, the scattering matrix containing the local reflection and transmission matrices for the upward and downward going waves can be introduced, yielding

U P ; (bottom) D;+i(top) DOWN UP DOWN UP v.y+i rJ.j+i Dj (bottom) UPj + 1(top) (2.35)

This equation expresses the concept of upgoing and downgoing waves being reflected and transmitted at each layer interface, as illustrated in Figures 2.2 and 2.3.

Dj.(bottom)

UP,-+1 (top)

Laf-erj +1

r r S „ N DOWN* ,, . .UP rrp; U P ; (bottom) = Tj J + 1 Dj (bottom) + t;- ;+ 1 U P7 +, (top)

(39)

\ ~ \ _{\ .DOWN} \ hj+i

f r

up

^«v

/ r;.y+l \ / UP,+1 (top) R , . .DOWNfï D>+, (top) = t;- ;+1 D7 Layer] Layerj +1 D7+1 (top)

(bottom) + ry,.,-+iUPj+i(top)

Figure 2.3 The composites of ihe downgoing wavefield at the top of layer j +1.

The reflection matrices in equation (2.35) are denoted "local" because they describe only the reflection and transmission response at a single layer interface, without accounting for the contribution of the other layers to form the total wavefield. These contributions are considered in Section 2.7, where "global" rather than "local" reflection matrices are introduced.

The elements of the local reflection and transmission matrices can be identified with the reflection and transmission of compressional (P), vertically polarized shear (SV) and horizontally polarized shear (SH) waves, as was shown by du Cloux (1986). They are given by

DOWN DOWN _ rPPJ.j+\ rPSV,j,j+\ u DOWN DOWN _ rSVP,j,j+l rSVSV,j,j+] u DOWN u u rSHSH.j,j+\ = r P L+i P i . ; + P l . ,+i P 2 . , r1[ P l . /+i P 2 .J- P l y+i F i , , ] . (2.36) which denotes the local reflection matrix for the downgoing waves incident at the interface from

layer/, DOWN

(40)

UP r;.>+i " UP UP „ rPP,j,j+) rPSV,j,j+l U UP UP „ rSVP,j,j+l rSVSV.j,j+\ u „ „ U P u u rSHSH,j,j+\ . P L - P I . M + PT.>P2.M -1 T ] [ P ; , P T (2.37)

which denotes the local reflection matrix for the upgoing waves incident at the interface from layer y+1, DOWN 7.7+1 0 DOWN DOWN tpPJJ+l lPSV,j,j+\ DOWN DOWN „ tsvpj.j+i tsvsv.jj+i u 0 0 _{'silSll.j.j+i}DOWN 2 [ P l y P].J + 1 + P1',yP2,;+ 1]"1. (2.38)

which denotes the local transmission matrix for the downgoing waves incident at the interface from layer), and

UP V. 7+1 UP fpp. 7. UP 'SVP, j 0 7+1 j'+i UP 'PSV, j , UP lSVSV, j 0 7+1 .7+1 0 0 UP hnsii. j 7+1 2 [ P 2 , y+l P l . ; + P l , ;+l P2.7] (2.39) which denotes the local transmission matrix for the upgoing waves incident at the interface from layer y+1. In these equations, use has been made of the symmetry properties of the inverse of the decomposition matrix (equation (2.23)).

2.6 The propagation of the wavefield within a single layer

The above equations describe the behaviour of the wavefield when it strikes an interface. Before the wavefield in the medium can be determined, propagation of the wavefield within a layer also has to be described. This relation, which is a simple phase-shift operation, can be easily inferred from equations (2.30) and (2.31), yielding

(41)

e\p[-i(oQj(x3-x3)]

t\p[-iü)Qj(xj-x3)]

_UPj(4)

(2.40)

where Q is defined in equation (2.16), and in which x3 and x3 denote the x3 coordinates of two points that are located within layer). With the aid of the equations derived in this section, the recursive formulae can be calculated, to facilitate the determination of the wavefield at any point in any layer.

2.7 The global reflection response 2.7.1 Backward recursion

Now a backward recursion procedure for the problem is performed according to the standard reflectivity method. A global reflection matrix is defined that relates upgoing and downgoing waves at the bottom of a layer j+1:

— up ^

UPy+1(bottom) = R;-+1 Dy+](bottom) (2.41) and thus (using equation (2.40)),

UP

UPy+](top) = exp(io>Q;+,/iy+j) R;+, exp(ifi>Qy+1/i;+1) Dy+1(top). (2.42)

The boundary conditions at the interface of layers; andj+1 are the continuity of the three components of particle velocity and traction. These are expressed by equation (2.34), which can be written as (see equation (2.35))

-—- DOWN ~ UP -~~^ UPy (bottom) = r, y+] D; (bottom) +ty y+1 UP;+1(top) ^v DOWN *"* UP "—""

D,+ )(top) = ty,y+1 Dy(bottom) + ry y+1 UPy+,(top) .

(2.43)

(2.44)

(42)

. T S „ x f DOWN „UP .. „ . VT>UP . . ry , „

U P; (bottom) = ( r ; .J + 1 +tyii+1expOö)Q;+,^+i)Ry+1exp(ift)Qy+i/i>+1) x

r UP UP ~\ DOWN 1 ""*

x [I-rJiy+iexp(jwQJ+]/!J+i)R;+1exp(iüjQJ+1/i;+1)J ty ; + I j D, (bottom), (2.45) up

so that the global reflection matrix of thej'th layer, R, , can be determined using the recursive formula:

„ U P DOWN .UP ,. „ , . „ U P .. „ , . Rj = Th /+i + lJ. y'+i expaö)QJ+1AJ+1) Rj+1 exp(iCt>Q,+1/i;+1) x

r* UP ,■ ^ . V „ U P ,• /-, , a „DOWN

X [ ^ • " y . y + l ^ P ^ Q M ^ + ^ ^ + l ^ P f ' ^ Q y + l V ^ J Vw+l ' (2-46) The recursion is initialized by using the knowledge that there are no upgoing waves in the

lower halfspace, and thus

R ^ ! = 0 • (2.47)

2.7.2 Forward recursion

In the forward recursion, the source is included by accounting for the surface traction components at the top of the configuration. This is why a new formulation of the forward recursion scheme is needed, which differs from the recursion that can be found in the literature, where either a halfspace or a free surface is present at the top of the configuration (e.g. du Cloux 1986).

From equation (2.24), it follows that

Tj^BjVPj+A^Dj , (2.48)

in which T ■ is a vector containing the three components of the traction,

Ti= ( r3 J, r1 <, - , T2 ) /)T. (2.49) the matrix B ,• is given by

(43)

B

'-7?

, 1

1.7 (2.50)

where P ■ denotes the q row of matrix P „ ,, and

1 V2 ~r2 J P2 p3 (2.51)

Equation (2.48) is applied at the surface of the configuration to obtain a relation between the surface tractions and the upgoine and downgoing wavefields. First the vector S T is defined, containing the traction components at the surface of the earth:

ST = ( 5 7 3 , 5 7 - , , ST2)7 = ( T ^ T , , f2) * „• (2.52)

Since the vector S T contains all characteristics of the surface source, the wavefield needs to be expressed in terms of this vector. For this purpose, equation (2.48) is applied at the top of the configuration, and subsequently rewritten as

Dj(top)= - A i B1U P1( t o p ) + A , S T , _(2.53) in which A , = -

V2

/i,4 _{] " i}

y

-!__•> 2 %

--pyP,Mr,^ '-**

-ip \CsJ I Pi -2P2Ys,i(Pi7P,J

Vl

ViPjJL--**

(p^y^-pH-j—zp \(pjsj/2

%

■yC/'.r,,/^.

4

.

(2.54) and

(44)

A , B , =

A,

-(-i— 2p

2

**\ + *p\**

A

7

sA -4P(-L-2P2\YPI1YJ'

Vs.l ) \Cs,l 1

,]/2

**-4p|-i— *py**

p

.xr,.y

2

i-r-V| - V r „ 7 , ,

*.1 / \Cs . l

i , A % / i

1/ ; , .»2 2\ . 2

* 1

(2.55) where A1 is the so-called Rayleigh denominator, and is given by

A

\ = {-T~

2p2

\ +

A

p\V'.i- (

2

-

56

)

The components of the slowness vector that appear in these expressions are defined in Section 2.4. In view of the structure of equation (2.53), the global reflection matrices relating the downgoing and upgoing waves at the top of layer j are defined as:

^ nnwN -—- ST -—•

D, (top) = R; U P; (top) + R f A, S T . (2.57) ST

The matrix A] is not included in the reflection matrix Ry because Aj is dependent on both px

and p2, and not simply on the slowness p. In the transformation to the space domain, the

DOWN ST

separation between A] and the reflection matrices R,- and R, (which depend only onp) is essential when cylindrical coordinates are introduced.

Using the same procedure as applied in the backward recursion ( applying the boundary conditions at each layer interface using equations (2.43) and (2.44), and successively eliminating D • (top) and U Py (top)), the following recursive scheme is obtained:

„DOWN .DOWN , . , . . . „ D O W N .. „ , . R>+, =«j,y+i expOo)Qjhj)Rj expQcoQ/ij) x

rT DOWN .. „ , .„DOWN .. , . , . -i~ IT UP

(45)

DS T .DOWN , . _ , . f , „DOWN .. „ , .

RJ + 1= tyy+, txp(icoQjhp [ I +R, expOwQjhp

. ST

x [ l - r^avQuQjhpR^ctXiwQjhp] r™» exP(/a> Q, h]) } Rf (2.59) The recursion can be initialized using equations (2.53) and (2.57):

„DOWN

Ri = - AiB] , (2.60)

and

ST

Ri = I (2.61)

2.8 The particle velocity and traction components at an arbitrary point in the medium

To obtain the wavefield at any given point in the medium, the forward and backward recursion have to be coupled to each other at the given depth level. As the expressions derived so far are all in the wavenumber domain, the horizontal spatial coordinates of the observation point only enter the calculations when the wavefield is transformed back to the space domain. An expression is wanted for the wavefield in the k layer, at depth X3 * . The vertical distance of the top of layer k to the observation point is denoted by hk :

, obs obs

hk = x2.k - * 3 , * - i - ( 2 . 6 2 )

as illustrated in Figure 2.4.

From the backward recursion, the following expression is obtained for the upgoing wavefield at depth x? J (equations (2.40) and (2.41)):

T T S , obs . „ , ofc.s UP , . „ ., .obs■ * obs. , , , „

(46)

Observation point , , obs depth x3k , obs hk , , obs hk~hk Layer k - 1 Layer k Layer k + 1 c3.*-l c3 . *

Figure 2.4 The observation point in layer k, and the notation associated with its depth level.

obs

From the forward recursion, the corresponding downgoing wavefield at depth level x3 k is

obtained (equations (2.40) and (2.57)):

£i , obs ^ r. .~ , obs, „DCWN . . „ ,obs,rTt^ , obs.

Dk(x3,k) = txp[icoQkhk ]Rk txp[icoQkhk ]UPk(xXk) +

+ exp[icoQkhlbs]R!?A1ST. _(2.64)

obs

From these two equations, the upgoing and downgoing wavefields at depth level x3<k,

--*^^ obs ^s> obs "^^**'

U P* (x3j(k) and Dk (x3 k ) , can be expressed in terms of the surface traction vector S T :

T>k(xl*ks)=[l-exp[i<oQkh°kbs]Rk ™ exp\io>Qkhk]Rk txp[icoQk(hk- h°kb*)] ] x

x txp\ia>Qkh0k s]Rk A] S T , _(2.65)

(47)

u pi (*3.t) = expL/6)QA Vik-lt°k")] R)ïFzxp[i(üQk (hk- hf")] x

x [ i - e x p l / ü J Q i / ï f ' l R ^ ^ ' e x p l / C ü Q ^ J R ^ e x p t / ö Q ^ / ^ - ^ ^ ) ] ] x

x exp[iü)QA/if*]Rf A , S T . (2.66)

Symbolically, these equations can be written as

Dt( ^ ) = R DA0 « : ^ ) A1S T , (2.67)

and

Ö P ^ U ) = RL't Cx-3^) A, ST , (2.68)

in which the global downgoing and upgoing reflection matrices RD^and RUj. are introduced, and are given by

- l _ _ obs. r . r- r\ i °bs, _ DOWN r. .-. , Tr-lfP ,. „ ., , r f i . , l

RD*U3.t)= [ I - e x p [ i w Qt/ it ]R* expDcoCMJRt expliü)Qt(/it-/it )J J x

x e x p [ / a ) Qt/ ! fs] R f , (2.69) and

RUt (*3?/) = expfüuQ* (/(*- / / f ■')] R^expfiwQi ( / it- /if*)] x

x [l-exp[/üjQA/j^ 5] Rt ' txp[iwQkhk]Rk txp[icoQk (hk- h°k ")] ]

-ï x

x exp[/ö)Qt/jf-']Rf , (2.70) respectively. The global downgoing and upgoing reflection matrices RDt and K Ut directly

(48)

Now the wavefield in terms of its particle velocity and traction components can be determined by applying the composition matrix P to equations (2.67) and (2.68). This composition, expressed by equation (2.24), results in

F~ / °*J\ r O < obs. if, , obs. -2, , obs ^-7

i(*3.*) =iV3(x3,k),T1(x^k),T2(x3k)]

= ^ = Pl t 4[ R UtC x # ) + R DtC x # ) ] A , ST , (2.71)

F"^ , obs s r £, , obs. f) . obs. rj , obs. , T

2 C*3.t ) = [ r3 Cx3.t) > vi C*3.*). v2 C*3l*) ]

= ^ P2.t[ R U4( x ^ ) - R DJ l( x ^ ) ] A1g T . (2.72)

From equations (2.71) and (2.72), it follows that the components of particle velocity and traction at any depth level, in the wavenumber-frequency domain, are given by the multiplication of the surface tractions, and a matrix that will be referred to as the Green's matrix G „ ( a = l , 2 ) . Thus,

r i obs \ p! i obs \ cr

Fa(*3,*) = Ga( x3 i it ) S T , (2.73)

in which the Green's matrices are given by

G ,(*36;) = ^ PU[ R Uk( x f j ) + RDk(x°3bks)] A , , (2.74)

and

G2(x?Sk) = y= P2>* t R U t O r ^ ) - RBk(xf;) ] A , . (2.75)

Expressions are now derived for the components of particle velocity and traction in the spacc-frequency domain in terms of the tractions acting at the surface of the medium. The reflection matrices appearing in these expressions are computed by means of the recursive scheme described above.

(49)

2.9 The transformation to the space domain

In the previous section, it is shown that the components of particle velocity and traction in layer k at depth level .r3 k' , contained in the vectors Fj and F2, can be written in the frequency-wavenumber domain as a multiplication of a Green's matrix and the traction distribution at the surface of the earth:

FB(*3Ü') = Ga( x # ) S T , (2.76)

in which

p . obs. ~ obs C; obs C; . obs. , T

F lfr3,* ) = f V3C*3.* ) . Flfr3.it ) - Tl(x3,k ) J . (2,77) ifr3.*) = ^ P i , J R U * f r3, * ) + R D* f r 3 , i ) ]Ai . (2.78) F2fr3,t ) = [ F3fr3,* ) , Vi(x3 i t) , V2fr'3,* ) ] , (2.79) G jfrf*) = yy P 2.* [ R U*fr3ft ) - RDt(i3J;) ] A j , (2.80) and ST = [ 73, 7 ] , 72]a l^3 = 0 • (2.81)

The global reflection matrices relating the downgoing and upgoing waves to the surface tractions, RD^.and RU^, respectively, are given by equations (2.69) and (2.70), in which the global reflection matrices can be computed by a recursive scheme given by equations (2.46), (2.58) and (2.59). The remaining two matrices, P„ k and A,, depend only on the medium

parameters of layer k and layer 1, respectively. They are given by equations (2.21), (2.22) and (2.54).

Since a multiplication in the wavenumber domain corresponds to a convolution in the space domain, the components of traction and displacement at position x{ are given by

ob ? f°° [°° r, , obs obs obs

(50)

where the Green's matrices in the space domain follow from

n . obs obs obs

« <x(*3,*> x \,k ~ xh x2.k~x2 > ~

= (2K)~ J J Ga{x03_sk,ki,k2)cxp[ikp(x°p^-xp)]dkldk2. (2.83)

To appreciate the symmetry present in the configuration, and to exploit this symmetry to arrive at a more elegant wavenumber-space transformation, the following substitutions are made. First, since the Green's matrices are derived in terms of slownesses pa, the substitution

ka=capa (2.84)

is made, and the analysis is restricted to positive frequencies. The time domain expressions can be obtained from Appendix A, equation (A.47). Then, the cylindrical symmetry in the earth model is accounted for by substituting

\p\=P cos(£) \ . - , • (2.85) | p2 = psinCr) 2 2 ' / wherep = (p} + p2) » an^L x°,k~xi = rcos(0) obs . tA (2.86) x2,k~x2= ''Sin(0)

in which r= [ (x°ks-X\) + (x2,k~x2) ] T h e notation employed for the space

coordinates is illustrated in Figure 2.5.

As the desired integration over the angle % cannot be recognized if the recursive scheme depends on X' it is first noticed that, as could be expected for a cylindrically-symmetric earth model, the reflection matrices are indeed independent of X- Furthermore, it can easily be verified that the symmetry present in the local reflection and transmission matrices is also present in the reflection matrices RUt and RDt:

(51)

RU, = | Ru pp Ru psy O RU SVP RUSvsv O O O RU SHSH (2.87) and RD _PP RD, _{RD svp} RDpSV O RDsvsv O D r , 'SHSH (2.88) Earïh Surface Layeri Layer k ■. 1

Theoretical and practical aspects of the Vibroseis method

Theoretical and practical aspects of the

Vibroseis method

Theoretical and practical aspects of the

Vibroseis method

Theoretical and practical aspects of the

Vibroseis method

Proefschrift

ter verkrijging van de graad van doctor aan de

Technische Universiteit Delft, op gezag van de Rector

Magnificus, prof. drs. P.A. Schenck, in het openbaar te

verdedigen ten overstaan van een commissie

aangewezen door het College van Dekanen op

donderdag 11 mei 1989 te 14.00 uur

door

Guido Jozef Maria Baeten

mijningenieur

Dit proefschrift is goedgekeurd door de promotor

Prof. A.M. Ziolkowski M.A., Ph.D., M.Sc. (Econ.)

Table of contents

Notations and conventions

1

Chapter 1 An introduction to Vibroseis

Reaction mass

Piston

/Fluid pressure

from

servo valve

Baseplate

TIME DOMAIN FREQUENCY DOMAIN

dm

= e-(o =/"'(o .

r\

=

hzf±

,

*(/>=

,

- o.io)

e r

T > - ^ . a.n)

feedback signal

f(0

vibrator

Pre-determined

sweep input q(t)

Chapter 2 The wavefield of a surface source in a horizontally

layered elastic medium

k>\

4ï

M

Mï)

:A

-2p

nm

PM

'Ah

JP-W

V*

2

p{pyf

\^-

P

-pApyf

-pi{pyf

1

/ \A

*f

/ \A

AP)

P

iM

n

4M'

-

Mn

M

\-

A

Pi

_UPj(4)

**\ + *p\**

**-4p|-i— *py**