GB2048481A - Seismic exploration system - Google Patents

Abstract

A method of determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth comprises employing one or more first point sound sources and one or more second point sound sources to produce first and second waves, the energy of the elastic radiation of the or each first source differing by a known factor from the energy of the elastic radiation of the or each second source, detecting reflections of said first and second sound waves from within the earth and generating therefrom respective first and second seismic signals and subjecting these two seismic signals to analysis and comparison. The invention also provides apparatus for performing the method.

Description

SPECIFICATION Seismographic method and apparatus This invention relates to a method of determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface layers in the earth and to apparatus for this purpose.

There is well-known a seismic reflection technique which employs a sound source at or near the earth's surface to emit an impulsive sound wave at a known time. As this sound wave passes through the earth it encounters boundaries between the different sub-surface layers. At each boundary some of the sound is transmitted and some is reflected. A receiver at or near the surface close to the source detects the reflected waves which arrive at later and later times.

A record (or seismogram) made of the receiver response is then processed to determine the amplitudes and arrival times of the individual reflections. These may then be used to determine the locations of the rock boundaries within the earth and/or the acoustic properties of the rock either side of each boundary.

The accuracy of such analysis depends on the ability of the processing technique to separate the individual reflections from one another. One of the reasons that this is a non-trivial task is that it is extremely difficult to generate a purely impulsive sound wave. The sound wave generated by most seismic sound sources has a duration which is longer than the smallest separation time interval of which the recording system is capable. In other words, the series of reflections which arrives at the receiver is not the series of impulses one desires (the reflectivity series); it is a series of over-lapping wavelets. The processing step which is used to remove the effect of the source from the recorded signal in an attempt to recover the reflectivity series is usually known as deconvolution.

The conventional description of a seismic signal regards the propagation of seismic waves as a linear elastic process in which the signal x,(t) is obtained as the convolution of the impulse response of the earth g(t) with a far field source wavelet s(t). Usually some additive noise is also present so that x1(t) = s(t)"g(t) + n1(t) (1) where the asterisk * denotes convolution. One wishes to extract g(t) uncontaminated by either s(t) or n,(t). However, n,(t) is not normally known, and often s(t) cannot be measured or predicted and must also be regarded as unknown.

Since s(t), g(t) and n,(t) are all unknowns, the problem of finding g(t) from the measurable quantity x,(t) is basically that of solving one equation containing three unknowns. It cannot be done, of course. Even when the noise can be ignored the essential difficult remains; that of deconvolving s(t) and g(t). Unless s(t) is known g(t) cannot be found without a lot of assumptions.

The noise term from n,(t) will be small compared with the signal term s(t) * g(t) provided there is enough signal energy. In order to achieve an adequate signal-to-noise ratio it is sometimes necessary to repeat the experiment a number of times in the same place, up to a total of, say, p times, using the same sound source, or identical sound sources. The series of received seismic signals x,(t), x2(t), . . xp(t) is summed to form a composite signal x(t), where

Equation (1) then becomes x(t) = p.s(t)*g(t) + n(t), (1 a) where n(t) is a composite noise signal, given by

This summation is known as a "vertical stack" and p is the "fold of stack" and is an integer greater than or equal to 1.The same, or similar, result may sometimes be obtained by generating p identical impulsive sound waves simultaneously. If the sources do not interact, each one will generate an identical far field wavelet s(t), and the received signal x(t) will be described by equation (lea).

For more than twenty years much ingenuity has been devoted to devising methods for solving equation (1) or (la) using assumptions which are as realistic as possible. But the fact remains that these assumptions are made purely for mathematical convenience. They are not substitutes for hard information.

The best known example of such a method is the least-squares time-domain inverse filtering method used throughout the industry. For this method to be valid it is required that (1) g(t) be a stationary, white, random sequence of impulses; (2) s(t) be minimum-phase and have the same shape throughout the seismigram; (3) there be no absorption.

All these assumptions are very strong, and they must all be correct simultaneously if the method is to work. This condition is very difficult to satisfy, especially since the assumptions are not mutually reinforcing. For example in attempting to satisfy the stationarity assumption, some sort of spherical divergence correction must first be applied. This has the effect of distorting s(t) unevenly down the seismogram which immediately invalidates the assumption that the shape of s(t) remains constant; it also introduces a tendency for s(t) to be non-minimum-phase in the early part of the seismogram.

According to a first aspect of the present invention there is provided a method of determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth which method comprises employing one or more first and one or more second point sound sources to produce respectively first and second sound waves, the energy of the elastic radiation of the or each first source differing by a known factor from the energy of the elastic radiation of the or each- second source, detecting reflections of said first and second sound waves from within the earth and generating therefrom respective first and second seismic signals and subjecting these two seismic signals to analysis and comparison.

According to a second aspect of the present invention there is provided apparatus for determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth which apparatus comprises one or more first point sound sources and one or more second point sound sources adapted respectively to produce first and second sound waves in the earth, the energy of the elastic radiation of the or each first source differing by a known factor from the energy of the elastic radiation of the or each second source, receiver means for detecting reflections of said first and second sound waves from within the earth and generating therefrom respective first and second seismic signals, and means for analysing and comparing said first and second seismic signals.

In the present invention the first and second sound sources may be individual point sound sources or there may be employed a plurality of identical non-interacting point sound sources which will produce a seismic signal having a greater signal-to-noise ratio.

Alternatively the receiver means may be adapted to sum a series of identical seismic signals obtained by repeated production of identical sound waves by one or more identical sound sources.

The term "point source" is employed throughout the specification to denote a source whose maximum dimension is small compared with the shortest wavelength of the useful radiation it generates. If this source is buried in an homogenous isotropic elastic medium it will generate spherically symmetric radiation at distances greater than about a wavelength. This is the far-field region in which any aspherical distortions of the wavefield from this point source will occur only at high frequencies outside the useful bandwidth.

The method of the present invention is suitable for buried point sources, on land and at sea. It requires none of the assumptions demanded by known methods. In particular, nothing is assumed about differences in the amplitude or phase spectra of s(t) and g(t). The present invention is based upon the fact that the wavelet obeys a scaling law of the type:

In this equation, , and T2 are both very nearly equal to r = t - r/c where t is time measured from the shot instant, r is the distance from the sound source to a point in the far field, and c is the speed of sound in the medium; s'(T2) is the far field wavelet of a source similar to that which generates 5(to), but which contains a3 times as much energy.

Fig. 1 shows diagrammatically how this source scaling law affects the far field wavelet.

There is excellent experimental evidence for the existence of such a scaling law for a variety of point sources and this law can readily be derived for explosives, for example, if the following assumptions are made: i. that the elastic radiation from the source possesses spherical symmetry; thus it will be applicable to most marine sources such as a single air gun, a single water gun, a marine explosive such as that available under the Trade Mark "Maxipulse", a marine source employing high pressure steam to cause an implosion such as that available under the Trade Mark "Vaporchoc", or a sparker and to explosives buried on land- but probably not to surface sources because their radiation is not spherically symmetric; ii. that the fraction of the total available energy stored in the explosive which is converted into elastic radiation is a constant for a given type of explosive and a given medium; iii. that the volume of the explosive may be neglected relative to the volume of the sphere of anelastic deformation produced by the explosion;; iv. that the elastic radiation produced by the explosion could be obtained by replacing the sphere of anelastic deformation by a cavity at the interior of which there is applied a timedependent pressure function P(t) and that P(t) is independent of the mass of the explosive and is constant for explosives of the same chemical composition in the same medium; and v. that " for an explosion of a first mass and T2 for an explosion of a second mass can be taken to be approximately equal and equal to T. This is sufficiently accurate if the time interval Ar between 1 and r2 iS unobservable within the frequency band of interest, i.e. Ar should be less than about one sample interval. This approximation will suffice for values of a up to about 5 or so.

To exploit the scaling law a seismic signal x(t) as described by equation (1a) is generated. The experiment is then repeated in the same place using a source of the same type but containing a3 times as much energy. This will generate a seismogram: x1(t) = qs'(t)*g(t) + n'(t) (3) where s'(t) is the far field wavelet of the source and is defined in equation (2); g(t) is the same as in equation (1) because it is the response of the earth to an impulse in the same place; the noise n'(t) may be different from n(t) in equation (lea); q is a known integer greater than or equal to 1 and which may be different from p in the equation (lea).

Let us consider these equations together for the case where the noise is negligibly small: x(t) = ps(t)*g(t) (4) x1(t) = qs(t)*g(t) (5) sl(t)= as(t/a) (2) In these three independent equations there are three unknowns: s(t), s'(t) and g(t). Therefore, in principle when the noise is negligibly small, we can solve for all three exactly without making further assumptions.

By taking Fourier transforms and by manipulation we obtain the equation p qS(af)= -S(f). R(f) (6) a2 Equation (6) suggests a recursive algorithm of the form:

where N is dicated by the highest frequency of interest and the process must be initiated with a guess at fO. If a > 1 equation (7) enables us to work up the spectrum calculating values at afO, a2 fO, . . . aNfO starting with a guess at fO.

To compute values at frequencies less than fO, equation (6) can be rearranged:

such that we obtain the recursion:

where M is dictated by the lowest frequency of interest. This now enables the values at frequencies f0/a, fO/a2, . . . , f0/aM to be computed.

Thus from the recursion scheme of equations (7) and (8) we can obtain values at frequencies fO/aM, f0/aM-l, . . , fO/a, fg, aft, . . . . . , αNf0.

We can now use an interpolation routine to find a value at another specified frequency, say f,.

and use the recursion to calculate values at af1, a2f, etc. This procedure is repeated until sufficient values have been computed. Once S(f) has been calculated, s(t) is obtained by taking the inverse Fourier transform.

It should be noted that the quantities involved in the algorithm are complex. One can operate either with the modulus (amplitude) and argument (phase), or with the real and imaginary parts.

The real and imaginary parts have been used in the example as these are considered to be the most 'basic' components of the complex numbers in a computer, whereas amplitude and phase are admixtures of these quantities, The Initial Guess The algorithm is initiated with a guess. If this guess is wrong, the final result will be wrong.

The guess at f0 is a complex number which, in all probability, will not be the true value at fO. In fact, the guess SG(fo) is related to the true value S(fo) in the following way:

wherein r e'O is the unknown complex error factor. If this error is not taken into account there will be generated the values:

which, with sufficient interpolation yield the function SG (fro) for fO/amKfsfO. the range can be extended to the origin by defining SG(O) = 0 which is compatible with a time series SG (t) with zero mean.

The effect of the initial error can be seen by substituting equation 9 into equation 10; thus:

It is evident that the error factor is constant for all the values deduced from the algorithm.

Thus far the algorithm has allowed computation of the function.

(12) where f0 has been assumed to be positive.

Two problems now exist. First the transform must be completed by generating values of SG(f) at negative frequencies.

Secondly the error factor must be found to obtain S(f) from equation (12). Both these problems can be solved by consideration of the physical properties of s(t), which impose constraints on the properties of S(f).

It is known that s(t) is real, and therefore the estimated wavelet should be real. This constraint imposes Hermitian symmetry on S(f). That is, the real and imaginary parts of S(f) must be even and odd functions, respectively. Thus if S(f) is known for positive frequencies, S(f) can easily be computed for negative frequencies using this condition.

However, only SG(f) is known, which is in error by a phase shift O and a scale factor r. The scale factor is unimportant because it has no effect on the shape of s(t), and consequently cannot affect our estimate of the shape of g(t). It can therefore be ignored.

However, the phase error 8 cannot be ignored, because this will make SG(t) non-causal, and it is known that s(t) is causal. That is, s(t) is zero for times t less than zero. In the frequency domain causality imposes the condition that the odd and even parts of the Fourier transform are a Hilbert transformer pair. It can be shown that this causal relationship is destroyed unless the phase error 9 is zero.

This consideration suggests a trail-and-error procedure for improving the estimate of s(t). This is as follows: 1. Compute SG(f) from an initial guess at f0 as described above, noting that SG(f) and S(f) are related as in (12).

2. Multiply SG(f) by a correction factor e-'BG where 0G is a guess.

3. Impose Hermitian symmetry.

4. Check for causality. If the recovered wavelet is non-causal, return to step 2 and repeat, using a different 8,. This procedure is repeated until the causality condition is met.

Thus the equations may be solved in the frequency domain using the algorithm described above and applying the constraints which allow from two physical properties of s(t): it is real and causal. The final estimate of s(t) will be in error only by a scaling factor r, which is trivial.

Having obtained a satisfactory estimate of s(t) , g(t) can be obtained using equation (1), by standard methods.

The algorithm described depends on a complex division in the frequency domain. There are two problems associated with this. First, the ratio will become unstable at any frequency at which the amplitude of the denominator is too small. Secondly, if the denominator contains nonminimum-phase components which are not contained in the numerator then the quotient becomes unstable in the sense that it is non-realisable.

To solve the first problem it is usual to add a small threshold of white noise to the denominator to negate the possibility of zero or near zero division. An alternative but more timeconsuming method is to search for low values in the denominator and to replace them with small positive values.

Finding the inverse of non-minimum phase wavelets is a well-known problem. However, the problem can be avoided simply by applying an exponential taper of the form e-to both x(t) and x'(t). By choosing y large enough the quotient R(f) can be forced to be stable, but then the estimates of s(t), s'(t) and g(t) will be distorted. In practice the distortion may be removed simply by applying the inverse taper el4 to these functions.

In the presence of noise the problem is to obtain a reliable estimate of the ratio spectrum R(f), for then the scaling law and recursive algorithm can be used to find S(f) as described above.

From equation (6) we define R(f) in the absence of noise as: X1(f) qS'(f) R(f) = ~~~~ = (13) X(f) pS(f) It follows that s'(')= r(t)*s(t) (14) where r(t) is the inverse Fourier transform of R(f) and, since s(t) and s'(t) are both real and causal, r(t) must also be real. However, r(t) will not be causal unless s(t) is minimum-phase. Both s(t) and s'(t) must be forced to be minimum-phase by applying the exponential taper to x(t) and x'(t) as described above. Under these conditions r(t) will be teal and causal.

In the noise-free case it is also true that: x'(t) = r(t)"x(t) (15) and it will be seen that r(t) is simply a one-sided filter which shapes x(t) into x'(t), provided the correct exponential taper has been applied. When noise is present the estimate of r(t) must be stabilized and this can easily be done using a least-squares approach (N. Levinson, in N. Wiener, 1947; Extrapolation, Interpoiation and smoothing of Stationary Time Series, Wiley, New York).

That is, a filter r'(t) is found which, for an output x(t) will give an output which is the best fit in a least-squares sense to x'(t). This filter r'(t) will be the best estimate of r(t).

In other words, in the presence of noise r(t) can be calculated in the time domain using standard programs, and then its Fourier transform taken, whence s(t) etc, can be found as described above.

It will be understood that although the problem has been discussed in terms of the scaled energies of the elastic radiation of the sources, normally the particle velocities or the sound pressures generated by the source may be detected and recorded using respectively a geophone or a hydrophone as conventionally employed.

It will further be understood that the individual elements of the apparatus of this invention may be chosen at will to be suitable for the particular purpose for which they are required thus air guns, water guns, "Maxipulse", "Vaporchoc", sparkers etc. may be employed as the sources. Similarly any suitable analysers, receivers etc may be employed as necessary.

It is believed that a may have a value of from 1.1 to 5, more preferably from 1.5 to 3.

EXAMPLE Applying the above wavelet deconvolution scheme to a synthetic example.

Two independent synthetic far field source wavelets as shown in Fig. 2 were generated. Each wavelet was calculated using a model as described in the Geophysical Journal of the Royal Astronomical Society 21, 137-161, for the signal generated by an air gun in water. The model is based on the nonlinear oscillations of a spherical bubble in water and takes into account nonlinear elastic effects close to the bubble. The model predicts waveforms which very closely match measurements.

The top wavelet s(t) of Fig. 2 was computed for a 10 cubic inch gun at a depth of 30 feet, a firing pressure of 2000 p.s.i. and a range of 500 feet from the gun. No sea surface reflection has been included. The bottom wavelet s'(t) was computed using the same computer program for a 80 cubic inch gun, at the same depth, firing pressure, and range. In other words, only the volume was changed.

Secondly each of these wavelets was convolved with the synthetic reflectivity series g(t) shown in Fig. 3. The result of performing these convolutions is shown in Fig. 4. The top trace x(t) represents the convolution of g(t) with the upper wavelet s(t) of Fig. 1; the bottom trace x'(t) represents the convolution of g(t) with the lower wavelet s'(t) of Fig. 1. Thus these two traces, x(t) and x1(t), were constructed entirely independently without any use of the scaling law.

It was then assumed that these two traces had been obtained knowing only that they were from the same place and that the top one was made using a 10 cubic inch gun, while the bottom was made using an 80 cubic inch gun at the same depth and pressure.

Since only the gun volume was changed, the scaling law can be invoked. In this case a3 = 8; therefore a = 2.

Solving for s(t) and g(t) as described above using the set of simultaneous equations (4), (5) and (2) with p = q = 1, substituting a = 2, the recovered wavelet and reflectivity series are shown in Fig. 5; they compare very well with the top wavelet of Fig. 2 and the original reflectivity series of Fig. 3. The small difference between the recovered series and the original are attributed primarily to computer round-off error.

This example shows that the method is valid in principle.

Thus by means of the present invention in-the absence of noise the impulse response of the earth can be obtained exactly. In.the pressure of noise a stable approximaton to this impulse response can be obtained, the accuracy of which approximation is dependent on the noise level present.

Claims (20)

1. A method of determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth, which method comprises employing one or more first point sound sources and one or more second point sound sources to produce first and second sound waves, the energy of the elastic radiation of the or each first source differing by a known factor from the energy of the elastic radiation of the or each second source, detecting reflections of said first and second sound waves from within the earth and generating therefrom respective first and second seismic signals and subjecting these two seismic signals to analysis and comparison.

2. A method according to claim 1, wherein a plurality of identical signals is summed to generate said first seismic signal.

3. A method according to claim 2, wherein said plurality of identical seismic signals is obtained by producing a series of identical sound waves by the use of one or more identical sound sources.

4. A method according to claim 1, 2 or 3, wherein a plurality of identical seismic signals is summed to generate said second seismic signal.

5. A method according to claim 4, wherein said plurality of identical seismic signals is obtained by producing a series of identical sound waves by the use of one or more identical sound sources.

6. A method according to claim 1, wherein a plurality of identical non-interacting sound sources is employed simultaneously to produce said first sound wave.

7. A method according to claim 1 or 6, wherein a plurality of identical non-interacting sound sources is employed simultaneously, to produce said second sound wave.

8. A method according to any of claims 1 to 7, wherein each sound source comprises an air gun, a water gun, a marine or a sub-surface land explosion generator, an implosive marine source or a sparker.

9. A method according to any preceding claim, wherein said known factor is from 1.33 to 125.

10. A method according to any preceding claim, wherein said known factor is from 3.375 to 27.

11. A method of determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth, substantially as hereinbefore described.

1 2. Apparatus for determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth, which apparatus comprises one or more first point sound sources and one or more second point sound sources adapted respectively to produce first and second sound waves in the earth, the energy of the elastic radiation emitted by the or each first source differing by a known factor from the energy of the elastic radiation emitted by the or each second source, receiver means for detecting reflections of said first and second sound waves from within the earth and generating therefrom respective first and second seismic signals, and means for analysing and comparing said first and second seismic signals.

1 3. Apparatus according to claim 12, wherein said first sound source comprises one or more identical sound sources arranged to produce a series of identical sound waves and the receiver includes means for summing said series of identical sound waves to produce said first seismic signal.

14. Apparatus according to claim 1 2 or 13, wherein said second sound source comprises one or more identical sound sources arranged to produce a series of identical sound waves and the receiver includes means for summing said series of identical sound waves to produce second seismic signal.

1 5. Apparatus according to claim 12, wherein said first sound source comprises a plurality of identical non-interacting point sound sources arranged to produce simultaneously said first sound wave.

1 6. Apparatus according to claim 12 or 15, wherein said second sound source comprises a plurality of identical non-interacting point sound sources arranged to produce simultaneously said second sound wave.

1 7. Apparatus according to any one of claims 12 to 16, wherein said known factor is arranged to be from 1.33 to 125.

1 8. Apparatus according to any one of claims 1 2 to 17, wherein said known factor is arranged to be from 3.375 to 27.

1 9. Apparatus for determining the location in the earth of sub-surface boundaries and/or the acoustic properties of sub-surface features in the earth, substantially as hereinbefore described.

20. The features hereinbefore disclosed, or their equivalents, in any novel selection.