Determination of Reflectivity and Water Depth for Firing Marine Sources as an Inverse Distortion Operator

This invention relates to marine seismic exploration and more particularly to the use of parameters measured from conventionally produced seismograms to control the firing of an array of sources to produce seismograms with enhanced reflections.

In marine seismic exploration, seismic energy can be generated from a towed array of sources. The resultant acoustic pressure waves are reflected from the subsurface formations and recorded as seismograms. Unfortunately, the seismograms include distortion which obscures the characteristics of the subsurface formations. In marine exploration, this distortion is often severe. It includes water bottom multiples as depicted in FIG. 1 and reverberations as depicted in FIGS. 2 and 3.

As the initial shot pulse travels to the bottom of the water, it is partially transmitted and partially reflected by this boundary. This transmitted part is reflected by the discontinuity beneath the water to become a primary reflection. The reflected part bounces back to the surface where it is again reflected to the water bottom. Once more, part of this pulse passes through the boundary and part is reflected. The process continues indefinitely. Each time a pulse bounces between the surface and the bottom of the water, a pulse emerges from the boundary at the bottom to follow the initial pulse into the section. These multiple pulses are delayed equal time intervals relative to each other, and their amplitudes are successively reduced by the reflection coefficient at the bottom of the water. The net result of this is to produce a seismogram consisting of a single primary followed by a train of equally spaced multiples or reverberations.

In U.S. Pat. No. 3,689,874 to Foster et al there is described a method and apparatus for removing the reverberation distortion present in seismic data by an operation which separates the reverberation distortion component of the data from the component representing the characteristics of the subsurface formation. This operation is known as dereverberation and more particularly involves the obtaining of a dereverberation operator in the form of an inverse filter that is applied to the seismic data to produce a signal which approximates the signal which would have been produced by a desired shot pulse interacting with the reflectivity of the subsurface formation in the absence of the reverberation distortion.

While dereverberation has been successfully used on many marine seismograms, it is not effective in all areas. Particularly where the water bottom is hard, the reverberations are so severe that the complete dynamic range of the recording system is required just to faithfully record the reverberations. Accordingly, it is desirable to suppress the reverberations before recording.

This invention relates to marine seismic exploration wherein the characteristics of the medium through which the seismic energy travels are measured from conventionally produced seismograms. These measurements are converted into control parameters. When the sources are fired under control of these parameters, seismograms having enhanced reflections are produced.

More particularly, this invention relates to measuring the reflectivity of the water bottom and the two-way travel time T.sub.w directly from seismograms produced by conventional shooting of at least one seismic source in a linear array of sources. Then, the measured reflectivity and travel time T.sub.w are converted to a three-point inverse distortion operator. The array of sources is then fired in three groups at a time sequence which produces this three-point operator. Firing of the second group of sources is delayed from the firing of the first group by the measured two-way travel time, and the firing of the third group of sources is delayed by twice the two-way travel time. The number of sources in each group is controlled by the measured reflectivity.

The foregoing and other objects, features, and advantages of the invention will be better understood from the drawings and the following more detailed description and appended claims.

FIG. 1 depicts water bottom multiples;

FIGS. 2 and 3 depict water reverberations;

FIGS. 4 and 5 show a marine seismic exploration system;

FIGS. 6-9 illustrate two- and three-point time domain inverse operators;

FIG. 10 depicts a conventionally produced seismogram combined with assumed values of R and T.sub.w ;

FIG. 11 depicts an example of a system for carrying out this invention;

FIG. 12 shows the correlation function of a selected window; and

FIG. 13 depicts conversion from an estimated reflectivity to actual reflectivity.

FIGS. 4 and 5 show a marine seismic exploration system. A vessel 11 traverses a seismic exploration path in surveying the subsurface formation 12 below the water layer 13. A linear horizontal array of seismic sources 14, 15, 16 and others is towed behind the boat by the cable 17. Surface floats help to maintain the equal spacing between the sources.

A hydrophone streamer 18 is also towed behind the boat. A surface support buoy 19 helps to maintain the proper relationship between the hydrophone streamer and the source array.

A seismic trace which is generated during seismic marine exploration by this system may be described mathematically as follows:

where,

S(t) = seismic trace,

B(t) = shot pulse,

H(t) = reverberation distortion,

R(t) = reflectivity function, and

N(t) = noise.

Equation (1) states that the seismic trace is produced by the shot pulse convolved with the reverberation distortion convolved with the reflectivity, plus coherent and random noise. [Convolution is denoted by an asterisk (*)]. Convolution is the process of filtering. Consequently, the seismic trace is the result of the shot pulse being filtered by the reflectivity of the earth and by the reverberation distortion.

Basically, the reverberation distortion is a function of both the magnitudes and signs of the reflection coefficients, R, at the bottom of the water and at the water surface. The reverberations of the shot pulse are delayed equal time intervals relative to each other, and their amplitudes are successively reduced by these reflection coefficients. The reflection coefficient at the water bottom can be either positive (+) or negative (-) in sign depending on the nature of the water bottom, while the reflection coefficient at the water surface is always negative (-) and is assumed to be unity (-1). Further, each multiple pulse is reduced in amplitude from the previous one by the magnitude of the water bottom reflectivity R. Accordingly, a dereverberation operator can be described by the two-point operator i(t):

This operator, consisting of a first impulse of unity magnitude at t = 0 and a second impulse of magnitude R at t = T.sub.w, can be illustrated in two forms. The first is shown in FIG. 6 for a positive water bottom reflection coefficient, +R. The second is shown in FIG. 7 for a negative water bottom coefficient, -R.

In those circumstances wherein the water bottom is relatively flat, the downward-going reverberation distortion at the shotpoint, and the upward-going reverberation distortion at the detector are approximately the same. Consequently, for a flat water bottom, a total dereverberation operator I(t) can be represented by the convolution of the two-point operator at the shotpoint, i.sub.s (t), with the two-point operator, i.sub.d (t), at the detector.

For the case of a positive water bottom reflectivity, +R, the dereverberation operator can be illustrated as in FIG. 8; while for a negative water bottom reflectivity, -R, the dereverberation operator can be illustrated as in FIG. 9. It can therefore be seen that the total dereverberation operator I(t) for a flat water bottom is a three-point operator consisting of three impulses -- the first of unity magnitude at t = 0, the second of magnitude 2R at t = T.sub.w, and the third of magnitude R.sup.2 at t = 2T.sub.w.

By firing the sources in various groups, with different firing times for each group, the resultant acoustic pressure wave has the time domain characteristics depicted in FIGS. 8 and 9. As an example, the operator of FIG. 8 is generated by firing a first group of pulses at time t = 0. The amplitude of the pulse generated by this group is arbitrarily designated unity. A second group of pulses is fired at time T.sub.w. The amplitude of the pulse produced by the second group is related to that of the first group by a factor 2R, where R is the reflectivity of the water bottom. A third group of pulses is fired at time 2T.sub.w. The amplitude of the pulse produced by the third group is related to the amplitude of the first pulse by a factor R.sup.2. It can be shown that firing the sources in a sequence such as the foregoing produces suppressed reverberations.

In order to fire the array of sources in such a sequence, it is necessary to determine the reflection coefficient R and the two-way travel time T.sub.w of acoustic energy through the water depth. In accordance with this invention, R and T.sub.w are obtained from a conventional seismogram produced by the firing of at least one of the sources 14-16 and others simultaneously. The sources are then fired in groups which produce a three-point inverse distortion operator. The seismic energy of this inverse operator produces seismograms which are substantially free of reverberations.

As one example of the manner in which R and T.sub.w can be determined from a conventional seismogram, consider FIG. 10. FIG. 10 depicts a conventional seismogram which has been combined, or stacked, in accordance with:

for different assumed values of R and T.sub.w. The seismic sections under the lines 21 . . . 22 represent reflections detected by the geophones 18 as a function of time (ordinate) after the shot. The seismic sections under the lines 23-25 represent the same seismic sections stacked with an assumed reflectivity R = 0.1 and for various assumed values of T.sub.w. The seismic section under the line 24 is the sum of the conventionally produced section plus the conventionally produced section weighted by a factor of 2(0.1) delayed by 96 msecs., plus the conventionally produced trace weighted by a factor of (0.1).sup.2 delayed by 192 msecs. The section under the line 25 is the sum of the conventionally produced trace weighted by the same factors but delayed by 98 msecs. and twice 98 msecs. The last section in this set under the line 26 is the sum of the conventionally produced trace weighted in the same manner but delayed by 120 msecs. and twice 120 msecs.

The sections under the lines 27-29 are combined, or stacked, sections with the same time delays, 94, 96 . . . 120 msecs. but with an assumed reflectivity of R = 0.2. The sections under the lines 30, 31 and 32 are stacked sections with an assumed reflectivity of R = 0.5, and with the same time delays. The values of R and T.sub.w can be determined from such an analysis by selecting the values of R and T.sub.w which produce the minimum energy in a section. Visual analysis of FIG. 10 shows that the section 33 has the minimum energy therein. This indicates that a correct value of R is 0.5, and the correct value of T.sub.w is 114 msecs.

The values of R and T.sub.w can be automatically determined in this manner but a very large number of computer operations is required, and generally it is not practical to make such an analysis on shipboard.

FIG. 10 illustrates an important feature which makes possible a simplified process for determining R and T.sub.w on board the ship. This process, which is subsequently described, is the invention disclosed and claimed in copending application Ser. No. 844,549, filed Oct. 25, 1977, Allen et al. FIG. 10 shows that the correct value of T.sub.w can be determined for any assumed value of R. That is, the section 34 has the minimum energy in its set and this indicates a value of T.sub.w of 114 msecs. even though an incorrect value of R = 0.2 was used in the stacking. Similarly, the section 35 has the minimum energy in its set and this indicates a correct value of T.sub.w = 114 msecs. even though an incorrect value of R = 0.1 was used in the stacking. This phenomena is utilized by assuming a value of reflectivity in determining T.sub.w. Then, reflectivity can be determined from the known value of T.sub.w.

FIG. 11 depicts one system for determining R and T.sub.w on board ship. A seismogram from a conventionally produced shot is recorded on magnetic tape 40 as is standard. A portion, or window of this seismogram is played back as indicated at 41. This window is auto-correlated by the auto-correlator 42.

The auto-correlation function is shown in FIG. 12. As is well known, auto-correlation requires the point-by-point multiplication of the samples in the window by themselves to produce a set of auto-correlation coefficients which are then summed to produce one point on the auto-correlation function of FIG. 12. Then, the window is time shifted, the samples are again multiplied one by the other and summed to form another point on the auto-correlation function. The process is repeated for various time shifts.

As an example, with zero time shift, the samples in the window are multiplied by themselves to produce the zero lag auto-correlation coefficients. These are summed to produce the value S.sub.0 on the curve of FIG. 12. Then, the samples of the window are shifted by one time sample, multiplied one by the other to produce the coefficients at one time sample lag. These are summed to form another point on the correlation function of FIG. 12. In general, the window is shifted with respect to itself by a number of time samples designated n.

The zero lag auto-correlation coefficients S.sub.0 are combined with the auto-correlation coefficients for a lag of n and for a lag of 2n for various values of n. The change in the value of n is indicated at 43 in FIG. 11. The combination takes place in the multipliers 44. The zero lag value S.sub.0 is multiplied by 2.0625. The value of the correlation function at a lag of n is multiplied by 2.5 and the value at a lag of 2n is multiplied by 0.5. These products are summed at 45. The minimum power, that is, the least sum, in the combined auto-correlation functions is selected at 46. The minimum power sum specifies the correct n which is designated M.

Since M is in sample times, this can be directly converted to T.sub.w because the time per sample is known. For example, the sections of FIG. 10 have a sample interval of 2 msecs. Assuming that a minimum power in the combined auto-correlation coefficients is detected when n = 57, then the output of 47 is a T.sub.w of 114 msecs.

The determination of reflectivity is an improvement on the technique described in Pflueger, "SPECTRA OF WATER REVERBERATIONS FOR PRIMARY AND MULTIPLE REFLECTIONS," GEOPHYSICS, Vol. 37, No. 5 (October 1972), pp. 788-796. In that prior art technique, reflectivity is related to the ratio of the amplitude of the side lobes of the auto-correlation function occurring at multiples of T.sub.w. This procedure provides a good estimate of reflectivity when the auto-correlation window is quite long. However, the use of a long auto-correlation window introduces other problems, notably the introduction of noise. It is for this reason that the window selected at 41 is a relatively short multiple of T.sub.w.

The determination of reflectivity as illustrated in FIG. 11 makes use of the Pflueger technique to determine a first estimate R. The divider 48 forms the ratio:

where S.sub.1M is the determined value of the correlation function at T.sub.w and S.sub.2M is the value of the auto-correlation function at 2T.sub.w.

It can be shown that the estimate R is related to true reflectivity R in the manner depicted in FIG. 13. F is the length of the correlation window in terms of the number of T.sub.w 's which are included in the window. As an example, suppose the auto-correlation window extends from 3T.sub.w to 7T.sub.w and R is determined to be 0.4 from Pflueger's technique. Then, F = 4 and FIG. 13 shows that R is approximately 0.425. In actual practice, reference is not made to a graph like FIG. 13. Rather coefficients are stored in memory which are applied to the estimate R to convert it to the true value of R. The coefficients are A.sub.0, A.sub.1 and A.sub.2 which are applied to R as indicated at 49.

The coefficients are a least squares fit to a surface described by plotting the true value of R as a function of N and the first estimate of R. Specifically, the coefficients are determined as indicated at 50. The value of F for this procedure is found by dividing the length of the window by T.sub.w as indicated at 51.

The length of the auto-correlation window, selected at 41, is a function of the two-way travel time as determined from the fathometer 52. The auto-correlation lag M is determined by dividing twice the fathometer reading by the velocity of sound in water, assumed to be 4800 feet per second. This ratio is determined at 52A. The length of the auto-correlation window is a multiple of the auto-correlation lag M. Generally, a window extending from 5M to 15M for each trace will be suitable.

Having determined T.sub.w and R, the sources are fired in groups which produce the desired three-point operator. The sources are fired in three groups designated a, b and c. The a group is fired at time t = 0. Group b is fired at t = T.sub.w, and group c is fired at t = 2T.sub.w. The firing circuits 53 for accomplishing this may be of the type shown in U.S. Pat. No. 3,687,218 -- Ritter.

The number of sources in each group is determined by the value of R. If a, b and c denote the number of sources in each group, these numbers are related to the determined reflectivity by the following:

a typical array which has been used in practice includes forty sources which are designated by the numerals "1" through "40". In this example, four sources are used as spares and the sources are fired in the following groups:

if R = 0.3

a = 40,39,37,36,33,32,30,29,27, ,23,22,19,18,13,12,10,9,7,4,3,1

b = 38,34,31,28,24,21,20,16,14,11,8,5,2

c = 25,17

if R = 0.4

a = 40,38,36,33,31,29,27,23,22,19,18,13,11,9,7,4,3,1

b = 39,37,34,32,30,28,24,20,16,14,12,10,8,5,2

c = 25,21,17

if R = 0.5

a = 40,38,36,33,31,29,25,22,18,14,12,10,8,5,3,1

b = 39,37,34,32,30,27,23,21,19,17,13,11,9,7,4,2

c = 28,24,20,16

if R = 0.6

a = 39,37,34,32,28,24,23,19,12,11,9,7,4,2

b = 40,38,36,33,31,30,27,22,20,18,16,14,10,8,5,3,1

c = 29,25,21,17,13

if R = 0.7

a = 39,37,34,30,27,22,18,14,11,8,4,2

b = 40,38,36,33,31,29,25,23,21,19,17,13,10,9,7,5,3,1

c = 32,28,24,20,16,12

if R = 0.8

a = 39,37,32,31,28,23,19,16,11,7,3

b = 40,38,36,34,30,27,24,22,20,18,14,12,10,8,5,4,2,1

c = 33,29,25,21,17,13,9

Because R changes slowly as the vessel moves, R is not recomputed for each shot. Rather, an optimum value of R occurring over the subsurface length of the cable is used to control the number of sources in each group for the three-point operator shots.

Conventional analog or hard wired digital circuits may be used in implementing the invention. However, the preferred embodiment of the invention makes use of a small general purpose digital computer which is carried on board the vessel. One example of such a computer which has been successfully used is the Texas Instruments 980A computer.

In the specific example just described, the control parameters are the reflectivity and water depth which are used to control the firing of the sources in time domain operators. Another example is where the directivity of the array of sources is controlled from measurements of dip made from conventionally produced seismograms. My copending application Ser. No. 793,680, filed May 4, 1977, describes an example of this. The disclosure of that application is specifically incorporated herein.

While a particular embodiment has been shown and described, other modifications are within the true spirit and scope of the invention. The appended claims are, therefore, intended to cover such modifications.