1 Geophysical Prospecting, 1999, 47, Marine seismic wavefield measurement to remove sea-surface multiples 1 A.M. Ziolkowski, 2 D.B. Taylor 2 and R.G.K...

Author:
Roy Chambers

0 downloads 2 Views 4MB Size

Marine seismic wavefield measurement to remove sea-surface multiples1 A.M. Ziolkowski,2 D.B. Taylor2 and R.G.K. Johnston2

Abstract We propose a new method for removing sea-surface multiples from marine seismic reflection data in which, in essence, the reflection response of the earth, referred to a plane just above the sea-floor, is computed as the ratio of the plane-wave components of the upgoing wave and the downgoing wave. Using source measurements of the wavefield made during data acquisition, three problems associated with earlier work are solved: (i) the method accommodates source arrays, rather than point sources; (ii) the incident field is removed without simultaneously removing part of the scattered field; and (iii) the minimum-energy criterion to find a wavelet is eliminated. Pressure measurements are made in a horizontal plane in the water. The source can be a conventional array of airguns, but must have both in-line and cross-line symmetry, and its wavefield must be measured and be repeatable from shot to shot. The problem is formulated for multiple shots in a two-dimensional configuration for each receiver, and for multiple receivers in a two-dimensional configuration for each shot. The scattered field is obtained from the measurements by subtracting the incident field, known from measurements at the source. The scattered field response to a single incident plane wave at a single receiver is obtained by transforming the commonreceiver gather to the frequency–wavenumber domain, and a single component of this response is obtained by Fourier transforming over all receiver coordinates. Each scattered field component is separated into an upgoing wave and a downgoing wave using the zero-pressure condition at the water-surface. The upgoing wave may then be expressed as a reflection coefficient multiplied by the incident downgoing wave plus a sum of scattered downgoing plane waves, each multiplied by the corresponding reflection coefficient. Keeping the upgoing scattered wave fixed, and using all possible incident plane waves for a given frequency, yields a set of linear simultaneous equations for the reflection coefficients which are solved for each plane wave and for each frequency. To create the shot records that would have been measured if the sea-surface had been absent, each reflection coefficient is multiplied by complex amplitude and phase factors, for source and receiver terms, before the five-dimensional Fourier transformation back to the space–time domain. 1

Paper presented at the 59th EAGE Conference – Geophysical Division, Geneva, Switzerland, May 1997. Received May 1997, revision accepted April 1999. 2 Department of Geology & Geophysics, University of Edinburgh, Grant Institute, West Mains Road, Edinburgh EH9 3JW, UK. q 1999 European Association of Geoscientists & Engineers

841

842

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Introduction The sea-surface causes all upgoing seismic waves to be reflected down again and makes the processing and interpretation of seismic data difficult. The theory of migration assumes that the data consist of primary reflections only. If multiples are present, they are treated in the migration process as primary reflections, thus confusing the interpretation. Multiples can also confuse amplitude-versus-offset (AVO) analysis, which uses subtle variations in the amplitude of primary reflections to deduce contrasts in the elastic properties at a reflecting interface. Thus a solution of the free-surface multiple problem is not only critical for the interpretation of the seismic data for basic geological structure, but it is also essential to the interpretation of more subtle effects caused by variations in fluid content and stratigraphy. Many different solutions to the problem of suppressing multiples in marine seismic data have been used in the last four decades, including two of the real workhorses of seismic data processing: common-midpoint stacking (Mayne 1950, 1962) and predictive deconvolution in the time domain (Peacock and Treitel 1969). Recently Taner, O’Doherty and Koehler (1995) extended the least-squares predictive deconvolution method of Peacock and Treitel (1969) to two dimensions: time and offset. Based on traveltime arguments, the theoretically correct range of application is thus extended from normal incidence to non-normal incidence over plane horizontal layers. This is equivalent to predictive deconvolution in the t–p domain (Tatham, Keeney and Noponen 1983; Tatham 1989). The first rigorous wave-theoretical approach to the suppression of multiples was formulated for 1D and 2D by Riley and Claerbout (1976), who first identified the seasurface as the source of difficulty, as illustrated in Fig. 1. Riley and Claerbout (1976)

Figure 1. Subsea reflections (a) with the effect of the sea-surface, (b) without the sea-surface.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 843

Figure 2. (a) Configuration of a line source and line of hydrophone receivers in a water layer overlying a stratified half-space; (b) as (a), except that there is no surface to the water layer.

stated: ‘The solution to eliminating this class of multiples amounts to finding a transformation which converts data recorded with the free surface present to the equivalent data which would have resulted if the free surface were absent’. They also identified the requirement that the source wavefield be known and, in the absence of source measurements, proposed a least-squares minimization of the energy of the multiple wavetrain to estimate the inverse of the far-field source waveform. We illustrate the effect of the free surface on synthetic data. Figure 2 shows a schematic of the model set-up with a line source and an array of hydrophones in water above a stratified elastic half-space: in Fig. 2a the water is a layer with a surface; in Fig. 2b the water is a half-space. The layered earth model is given in Table 1. Figure 3 shows the Ricker wavelet used as the source time function. Figure 4a shows the full synthetic seismogram for the configuration of Fig. 1a, calculated using the reflectivity method (Fuchs and Mu¨ller 1971); Fig. 4b shows the scattered field component of Fig. 4a, which is the total seismogram minus the incident field (the direct wave from the source and its sea-surface reflection). Figure 4c shows the full synthetic seismogram for the configuration of Fig. 2b and Fig. 4d shows the scattered field component of Fig. 4c, which is the total seismogram minus the incident field (the direct wave from the source). Removal of the incident field and the sea-surface renders the data far more interpretable: the primary reflections and P–S converted reflections can be seen clearly in Fig. 4d, whereas it is impossible to distinguish them in Figs 4a or b. Other previous work on this problem using a wave-theoretical approach and point sources includes that of Kennett (1979), Berkhout (1982), Fokkema and van den Berg

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

844

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Table 1. Layered earth parameters for Fig. 2 synthetic seismograms. Layer number

Thickness (m)

Depth (m)

P-wave velocity (m/s)

S-wave velocity (m/s)

Density (kg/m3)

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

70 100 150 200 150 200 150 200 150 200 200 100 200 250 200 150 200 Half-space

70 170 320 520 670 870 1020 1220 1370 1570 1770 1870 2070 2320 2520 2670 2870 ∞

1500 1900 2000 2200 2500 3000 2800 4500 3800 4500 6000 3500 3800 4600 4900 2800 2500 5750

0 950 1000 1100 1250 1500 1400 2250 1900 2250 3000 1750 1900 2300 2450 1400 1250 2875

1000 2040 2100 2200 2260 2600 2300 2300 2800 2300 2800 2600 2800 2550 2500 2400 2460 2610

(1990, 1993), Carvalho, Weglein and Stolt (1991, 1992), Verschuur, Berkhout and Wapenaar (1989, 1992), Verschuur and Berkhout (1992), van Borselen, Fokkema and van den Berg (1996), and Weglein et al. (1997). These approaches are discussed briefly below. Kennett (1979) solved the problem for a point source and a line source at the surface of a horizontally stratified elastic or fluid earth. Berkhout (1982) showed that the effect of the free surface is essentially a feedback problem, and the reflection response can be recovered by undoing the feedback effect of the free surface by a recursive process of prediction and subtraction. For the general case of a three-dimensional earth, a two-dimensional array of sources is required for each receiver and a two-dimensional array of receivers for each source. He considered the case of a dipole source at the water-surface for each shot. One problem not discussed by Berkhout (1982) (or by Riley and Claerbout 1976) is the removal of the incident field. Their analyses apply to the scattered field only. To obtain the scattered field from the measured data, the incident field must be removed. This is a non-trivial problem since the incident field and scattered field may overlap, as we have illustrated above. Berkhout’s (1982) scheme has been developed by Verschuur et al. (1989, 1992), Verschuur (1991) and Verschuur and Berkhout (1992). There are three elements of this development that should be mentioned. First, following Berkhout (1982), the incident field is ignored. In practice it is muted out (Verschuur 1991), and part of the

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 845

Figure 3. (a) Ricker wavelet used in synthetic seismograms; (b) its amplitude spectrum.

scattered field is removed in the process. What remains is thus an incomplete representation of the scattered field. Second, the theory assumes the source is a dipole at the surface, whereas, in practice, it is an array of airguns a few metres below the surface. Third, although the true source time function is in the formulation, the authors claim it does not need to be known. Instead, following Verschuur et al. (1989) (and the suggestion of Riley and Claerbout 1976) they assume the energy in the upgoing wave is a minimum when the multiples have been removed. The source wavelet is then

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

846

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Figure 4. (a) Synthetic seismogram corresponding to the configuration in Fig. 2a; as Fig. 4a but with the incident field removed; (c) synthetic seismogram corresponding to the configuration of Fig. 2b; (d) as Fig. 4c, but with the incident field removed.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 847

essentially redefined as the function that satisfies this minimum-energy assumption. In this formulation the unknowns are the coefficients of this wavelet, whose length is arbitrary. The problem of source directivity, arising from the reality that the source is an array below the surface rather than a dipole at the surface, is tackled in a similar way: additional source parameters are introduced and are found using the minimum-energy assumption. An exact wave-theoretical formulation for a monopole source below the sea-surface was presented by Fokkema and van den Berg (1990, 1993). For the monopole source a complete two-dimensional array of receivers is required and, because reciprocity is an essential part of the formulation, each receiver must also be a monopole, and for each receiver a complete two-dimensional array of monopole sources is required. The approach of Fokkema and van den Berg (1990, 1993) begins with Rayleigh’s representation theorem and arrives at an integral equation that needs to be solved. The formulation explicitly includes both the incident and scattered fields. The theory of Fokkema and van den Berg (1990, 1993) has been developed by van Borselen et al. (1996), who applied it to two-dimensional synthetic data. The integral equation is solved in the double Radon domain. The Fokkema–van den Berg–van Borselen scheme can be applied in two modes: either a direct solution to the integral equation, or by expanding in a series, which is essentially a prediction-and-subtraction approach. In both schemes the time function of the monopole source must be known and the incident field must be removed. Carvalho et al. (1991, 1992) presented a wave-theoretical multiple removal scheme, using an inverse scattering approach and assuming the source to be a monopole whose signature needs to be known. In a follow-up paper, based on the same theory, Carvalho and Weglein (1994) suggested that the signature of the monopole may be estimated using the energy-minimization criterion using simulated annealing to find the wavelet that meets the criterion. Ikelle, Roberts and Weglein (1997), using the same inverse scattering theory, described an iterative approach to find the wavelet to minimize the energy of first-order multiples. The more complete formulation of the inverse scattering approach by Weglein et al. (1997) states that the incident field must be subtracted from the total field to recover the scattered field and that therefore knowledge of the source is a prerequisite for the solution, but the formulation is still limited to a monopole point source. The objective of this paper is to present an alternative exact formulation for the removal of the incident field and the water-surface for a three-dimensional earth and seismic source arrays. Three problems associated with earlier work are solved: we can accommodate source arrays, rather than point sources; we can remove the incident field without simultaneously removing part of the scattered field; and we avoid the minimum-energy criterion to find a wavelet. The additional information required to solve these problems comes from measurements of the source wavefield made during data acquisition. In common with the previous work based on a rigorous wave-theory approach, this paper formulates the problem for multiple shots in a two-dimensional configuration

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

848

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

and for multiple receivers in a two-dimensional split-spread configuration for each source position. That is, the data are five-dimensional and must be adequately sampled in each dimension. In contrast to previous work, we regard the source as knowable, if appropriate measurements are made during seismic data acquisition (Ziolkowski et al. 1982). The unknown in our formulation is the earth and we formulate the problem for a realistic array of airguns, rather than for a monopole or a dipole. We first present the idea that the reflection response of a three-dimensional earth can be formulated as the ratio of the upgoing to the downgoing fields. This is not new. For one dimension, it is in Riley and Claerbout (1976, equation (1)). The complexity introduced by the water-surface is then analysed, first for a one-dimensional earth, and then for a three-dimensional earth. The solution is first derived heuristically, but a full wave-theoretical analysis is then presented. The analysis of the wavefield from a single shot is divided into two parts: the first part treats the wavefield that would exist in the absence of the sea-floor, which is here called the incident field, and describes how this wavefield can be computed from pressure measurements made close to the source. The second part describes the total measured field as the sum of this incident field plus the scattered field caused by the introduction of the three-dimensional earth forming the sea-floor. The incident field is known, therefore the scattered field may be isolated from the total field. The remainder of the paper is devoted to a wave-theoretical formulation of the steps needed to separate the plane waves and recover the response of the earth without free-surface multiples. The proposed scheme imposes severe constraints on the acquisition of the data: the data are five-dimensional and all five dimensions must be adequately sampled. This should come as no surprise, following previous work in this subject and the excellent book by Vermeer (1990), which demonstrates the symmetric sampling technique, required here, as theoretically optimum. Plane-wave reflection response: no free surface The problem is approached in terms of the plane-wave reflection response. Consider a horizontally stratified earth, as illustrated in Fig. 5a. The reflection response depends on the angle of incidence Θ and on the angular frequency q. Using the horizontal wavenumber kx ¼ (q sinΘ)/c1, where c1 is the speed of sound in medium 1, the reflection response may be written as ˆ ˆ 12 ðkx ; qÞ ¼ U1 ðkx ; qÞ ; R ˆ 1 ðkx ; qÞ D

ð1Þ

ˆ 1(kx, q) is the downgoing plane wave in medium 1 at frequency q and angle in which, D ˆ 1(kx, q) is the corresponding upgoing plane wave at the same of incidence v, U frequency, and the angle of reflection equals the angle of incidence because all the interfaces in the model are horizontal.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 849

(a)

D1 (k x , ω )

U1 (kx , ω)

θ θ

(b) D D D1 (k x , k y , zw , ω )

U U D D U1 (kx , k y ; kx , k y , zw , ω)

z = zw Figure 5. (a) Plane-wave reflection response of a stratified half-space: D1(kx, q) is the downgoing plane wave in medium 1 with angular frequency q and angle of incidence v, corresponding to horizontal wavenumber kx; U1(kx, q) is the reflected upgoing plane wave with the same horizontal wavenumber because the angle of reflection equals the angle of incidence. (b) Plane-wave reflection response of a three-dimensional earth, referred to a plane z ¼ zw above D the sea-floor: D(kD x , ky , zw, q) is the downgoing plane wave in the water with angular frequency q D and horizontal wavenumbers kD x and ky ; the upgoing wave may be decomposed into plane-wave U D D components, one of which is U(kU , k x y , kx , ky , zw, q).

Now consider a three-dimensional earth covered by water with no free surface, ˜ˆ (kD, kD, z , q) incident on the as illustrated in Fig. 5b, and a downgoing plane wave D w w y D horizontal plane z ¼ zw, taken to be in the water, above the earth. kD x and ky are the two horizontal wavenumbers of this downgoing wave. This downgoing wave is reflected by the three-dimensional earth. The upgoing reflected wave is, in general, not plane. However, it can be decomposed into plane-wave components, and at the plane z ¼ zw U D D U U ˆ˜ (kU one of these is U x , ky ; kx , ky , zw, q), in which the parameters kx and ky are its horizontal wavenumbers. We can now define a reflection coefficient at the plane z ¼ zw as the ratio of this component of the reflected upgoing plane wave to the incident downgoing plane wave: ˜ U U D D ˜ˆ U ; kU ; kD ; kD ; z ; qÞ ¼ Uˆ ðkx ; ky ; kx ; ky ; zw ; qÞ Rðk w x y x y ˜ˆ D ; kD ; z ; qÞ Dðk w x y

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

ð2Þ

850

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

and we can define similar coefficients for all the other plane-wave components of the reflected wave. The reflection coefficient defined by (2) is used as a boundary condition in our wave-theoretical formulation. We note that there are five independent coordinates U D D (kU x , ky , kx , ky and q) required to define this coefficient at depth level z ¼ zw. Plane-wave reflection response: free surface included We now introduce the sea-surface for the one-dimensional problem, as shown in Fig. 6. We assume that a finite number of discrete plane waves is available to define the response. Consider the ith plane pressure wave Pi incident on the sea-floor, where the subscript i denotes a unique frequency–wavenumber combination. The reflection response of the layered earth is then Ri, as defined in (1). The upgoing primary reflection just above the sea-floor is a plane wave RiPi and the angle of reflection equals the angle of incidence. This upgoing wave travels to the sea-surface, where it is totally reflected, and then travels down to the sea-floor again, with, in general, a different phase, which is accounted for by the factor si which includes the p phase shift at the seasurface. The phase factor si corresponds in the time domain to the time delay caused by the two-way propagation in the water layer plus the –1 reflection coefficient at the water-surface. This downgoing scattered wave is thus siRiPi with the same angle of incidence as the incident wave. This is therefore reflected by the sea-floor with the same factor Ri. The upgoing wave is reflected back down again, and so on, creating an infinite sequence of upgoing and downgoing waves.

π

π

si Ri Pi

Pi

si Ri si Ri P

Ri P Ri

π

si Ri si Ri si Ri Pi

Ri si Ri P Ri

Ri si Ri si Ri P Ri

Ri

Figure 6. A water layer over a stratified half-space, showing multiple reflections of a single planewave component: Pi is one incident plane-wave component, Ri is the reflection coefficient of the stratified half-space for this component, and si is the phase lag associated with the two-way propagation in the water layer plus the p phase change at the water-surface (corresponding to a reflection coefficient of –1).

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 851

The wavefield in the water therefore consists of three parts: an incident field Pi, an upgoing scattered field Ui which contains an infinite number of multiple reflections, given by the equation Ui ¼ Ri Pi þ Ri si Ri Pi þ Ri si Ri si Ri Pi þ Ri si Ri si Ri si Ri Pi þ . . . ;

ð3Þ

and a downgoing scattered field Di which also contains an infinite number of multiple reflections, given by the equation Di ¼ si Ri Pi þ si Ri si Ri Pi þ si Ri si Ri si Ri Pi þ . . . :

ð4Þ

From (3) and (4) the upgoing wavefield is related to the incident field and the downgoing scattered field simply by Ui ¼ ðPi þ Di ÞRi :

ð5Þ

We see the reflection response Ri is the ratio of the upgoing wave Ui to the downgoing wave, with the downgoing wave consisting of the downgoing incident wave Pi plus the downgoing scattered wave Di. We now introduce the free surface of the water at z ¼ 0 for the case of the threedimensional earth, as shown in Fig. 7. Again a single wave Pi is considered, incident on a plane z ¼ zw. The resulting scattered wavefield consists of a superposition of upgoing and downgoing plane-wave components, each of which contains an infinite number of multiple reflections. Consider now a single component Uji of the upgoing scattered field in which the subscript j indicates a unique combination of wavenumbers, not necessarily the same as those of the ith downgoing incident wave. This wave is caused by a reflection of the downgoing incident field Pi with reflection coefficient Rji, as defined in (2), and by reflections of each component Dki of the downgoing wavefield, in which the subscript k indicates a unique combination of wavenumbers. This may be expressed as X Rjk Dki ; ð6Þ Uji ¼ Rji Pi þ k

or, using frequency–wavenumber notation, ˜ˆ U U D D U D D ˆ˜ D D U˜ˆ ðkU x ; ky ; kx ; ky ; zw ; qÞ ¼ Rðkx ; ky ; kx ; ky ; zw ; qÞ:Pðkx ; ky ; zw ; qÞ XX ˜ U U ˜ˆ ˆ x ; ky ; kx ; ky ; zw ; qÞ:D Rðk þ ð2pÞ¹2 Dkx Dky kx

ky

D × ðkx ; ky ; kD x ; ky ; zw ; qÞ;

ð7Þ U D D ˆ˜ (kU in which U x , ky ; kx , ky , zw, q) is the upgoing scattered plane-wave component, ˜ D D ˆ D(kx, ky kx , ky , zw, q) is a plane-wave component of the downgoing scattered field, D P˜ˆ (kD x , ky , zw, q) is the downgoing incident plane wave, and Dkx and Dky are the wavenumber sampling intervals.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

852

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

0

y

x

Sea-surface : pressure = 0.

z U ji

Notional boundary z = zw

Uki

Pi Rji

Dki

Rki

Figure 7. A water layer overlying a three-dimensional earth: Pi is one incident plane-wave component, Rki is the reflection coefficient of this downgoing component for the kth upgoing scattered plane-wave component Uki, which, on reflection at the water-surface, becomes a downgoing scattered plane-wave component Dki; Uji is another possible upgoing scattered planewave component; all waves are referred to the plane z ¼ zw.

We have arrived at this result heuristically. Below we provide a wave-theoretical formulation to arrive at this same result. The wave theory also helps in formulating the acquisition requirements and data-processing steps. The incident field and the upgoing and downgoing scattered fields can be found from measurements, as described below. Keeping the upgoing wavenumber U components (kU x , ky ) and frequency q fixed, and varying the wavenumber components D D (kx , ky ) of the incident field, leads to simultaneous equations that can be solved for the ˜ˆ (kU, kU; kD, kD, z , q). reflection coefficients R w x y x y There are three issues to be addressed: (i) determination of the incident field; (ii) isolation of the scattered field from the measured data; and (iii) separation of the upgoing and downgoing components of the scattered field. These issues are discussed below in this order. There are also several practical issues of how to obtain the required multishot twodimensional split-spread seismic data, which are discussed briefly at the end of the paper. Analysis of a single-shot record: decomposition of the wavefield into incident and scattered fields It is convenient to divide the analysis into two parts: first, the wavefield that would be present without the sea-floor and second, the wavefield that exists when the sea-floor is present. In the first part, shown in Fig. 8a, the sea-surface is the only scatterer. There are no upgoing scattered waves because there is no scatterer below the source. The resulting downgoing scattered field is related to the upgoing field from the source by the

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 853

Figure 8. Separation of the total field in the water layer into two parts: the incident field and the scattered field. (a) The incident field exists alone in the absence of the sea-floor. (b) The total field is the sum of the incident field (a) and the scattered field, which are the acoustic waves trapped between the sea-surface and the sea-floor. The scattered field can be found by subtracting the incident field from the measured total field.

zero-pressure condition at the sea-surface. Provided the source wavefield can be measured, as discussed below, the wavefield for the water layer as a half-space can be determined. This wavefield, in the absence of the sea-floor, we now call the incident field. Second, we consider the presence of the sea-floor, as shown in Fig. 8b. The sea-floor and the earth beneath the sea-floor add an upgoing component to the scattered field in the water. The sea-surface reflects all upgoing waves, so this scattered field contains both upgoing and downgoing waves. We assume that the sea-floor is sufficiently far away to have a negligible effect on the way the source behaves. The total field in the water is then the field that would be present in the absence of the sea-floor, which we call the incident field, plus the scattered field introduced by the presence of the sea-floor. The recorded data are the sum of the incident field and the scattered field. These two fields must be separated. The incident field can be computed directly from the known configuration of the source and receivers, and can be subtracted from the total field to leave the scattered field introduced by the sea-floor. At the sea-surface the total pressure field is zero. Since the incident pressure field is zero at the sea-surface, it follows that the scattered pressure field is also zero at the sea-surface. This information allows the upgoing and downgoing components of the scattered field to be separated in the frequency–wavenumber domain.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

854

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Measurement of the source wavefield First, the source wavefield must be measured. Typical source arrays consist of airguns arranged singly or in clusters, each such volume-injection element being small compared with a wavelength and therefore emitting a spherical wave. If there are n elements in the source array there are n sources of spherical waves in the array. Each element emits an air bubble that oscillates non-linearly and is affected by the pressure field generated by the other bubbles. Nevertheless, each element is still a source of spherical waves and the pressure in the water is the superposition of these spherical waves plus their reflections off the sea-surface. Of course there is also a component of the pressure field caused by scattering off subsea interfaces. This component can be neglected, because of the effect of spherical divergence, if the distance from the source element to the point of measurement is very much smaller than the distance from the source to the sea-floor. In principle, and in practice (Ziolkowski et al. 1982; Parkes et al. 1984; Ziolkowski and Johnston 1997; Pieprzak et al. 1998), this condition can normally be met in normal continental shelf conditions if the gun–hydrophone distance is about 1 m. The principle of the method of measurement, using near-field hydrophones, is illustrated in Fig. 9. We assume a right-handed Cartesian coordinate system for the survey with the origin at the sea-surface and with the z-axis pointing downwards. To ease the discussion we assume the vessel towing the source is steaming parallel to the x-axis. The shot-point has coordinates (xs, ys, zs) and is, in fact, a point somewhere in the middle of the source array. We designate a point at the sea-surface vertically above this point on the array as the origin of relative coordinates for positioning the source array elements (airguns) and near-field hydrophones, with the x-axis pointing in the direction of the vessel and the z-axis pointing vertically downwards. In this relative

Figure 9. Principle of near-field measurement of the wavefield generated by an array of sources.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 855

coordinate system (xhj , yhj , zhj ) is the position of the jth hydrophone and (xgk, ygk, zgk) is the position of the kth source element. Ignoring the reflection from the sea-floor, the jth hydrophone measures the following pressure: n X rkj 1 gkj 1 qk t ¹ ¹ qk t ¹ ; ð8Þ pj ðtÞ ¼ c c r gkj k¼1 kj in which qk(t) is the source time function of the kth notional source element, c is the speed of sound in water, q ð9Þ rkj ¼ ðxgk ¹ xhj Þ2 þ ðygk ¹ yhj Þ2 þ ðzgk ¹ zhj Þ2 is the distance from the kth source element to the jth hydrophone, and q gkj ¼ ðxgk ¹ xhj Þ2 þ ðygk ¹ yhj Þ2 þ ðzgk ¹ zhj Þ2

ð10Þ

is the distance from the virtual image of the kth source element to the jth hydrophone. At least n near-field hydrophones are required for an n-element source array. Each hydrophone yields a different measurement defined by (8), because the rkj and gkj are different for each hydrophone. From the n pressure measurements of the source wavefield the n ‘notional’ source signatures sk(t), k ¼ 1, 2, . . . , n, can be found (Ziolkowski et al. 1982; Parkes et al. 1984). In principle, the source signatures may vary from shot to shot, so these measurements need to be made on a shot-by-shot basis. However, as we shall see, the theory requires the incident field to be the same from shot to shot. These source measurements can be used to determine the best mean incident field. Shot-to-shot deviations from the mean introduce noise. It has recently been observed (Henman, paper presented at UKOOA meeting, London, 1999; Strijbos, paper presented at UKOOA meeting, London, 1999) that, due to the focusing effect of the airgun array, the sea-bed reflection in water 90 m deep with a hard bottom can be detected on the near-field hydrophones with amplitudes of the order of a few per cent of the near-field signal. A major improvement in the measurement procedure has been devised by Ziolkowski (1998). This allows the pressure to be measured only a few centimetres from the gun ports – essentially inside the bubble – where the pressures are about 20 times greater than at the near-field hydrophones placed in the linear field. The relative amplitude of the sea-bed reflections is consequently about 20 times smaller and effectively negligible. The incident field: no sea-floor First we consider the part of the wavefield that would exist if the sea-floor were absent and we aim to obtain an expression for this field in the frequency–wavenumber domain in terms of the known notional source signatures. For a shot at source position (xs, ys, zs) the pressure in the water at position (x, y, z) would be p0(x, y, z, t; xs, ys, zs). This is

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

856

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

R the sum of the direct wave from the source pD 0 (x, y, z, t; xs, ys, zs) plus the wave p0 (x, y, z, t; xs, ys, zs) reflected from the sea-surface, thus

R p0 ðx; y; z; t; xs ; ys ; zs Þ ¼ pD 0 ðx; y; z; t; xs ; ys ; zs Þ þ p0 ðx; y; z; t; xs ; ys ; zs Þ:

ð11Þ

At the sea-surface, the pressure p0(x, y, z ¼ 0, t; xs, ys, zs) is zero. Therefore, from (11), R pD 0 ðx; y; z ¼ 0; t; xs ; ys ; zs Þ ¼ ¹p0 ðx; y; z ¼ 0; t; xs ; ys ; zs Þ:

ð12Þ

The direct wave satisfies the following inhomogeneous wave equation: 2 n X ∂ ∂2 ∂2 1 ∂2 D þ þ ¹ ðx; y; z; t; x ; y ; z Þ ¼ ¹ sk ; p s s s 0 ∂x2 ∂y2 ∂z2 c2 ∂t 2 k¼1

ð13Þ

in which sk ¼ qk ðtÞdðx ¹ xs ¹ xgk Þdðy ¹ ys ¹ ygk Þdðz ¹ zgk Þ

ð14Þ

is the kth notional monopole source in the array and c is the velocity of sound in water. Superposition applies for the notional monopoles (Ziolkowski et al. 1982). The reflected wave satisfies the homogeneous wave equation 2 ∂ ∂2 ∂2 1 ∂2 R þ þ ¹ ð15Þ p0 ðx; y; z; t; xs ; ys ; zs Þ ¼ 0: ∂x2 ∂y2 ∂z2 c2 ∂t 2 We now consider the transformation of the field from the time–space domain to the frequency–wavenumber domain, with a temporal Fourier transform followed by a two-dimensional spatial Fourier transform. The temporal transform and its inverse are ∞ po ðx; y; z; t; xs ; ys ; zs Þ expfiqtgdt ð16aÞ Po ðx; y; z; q; xs ; ys ; zs Þ ¼ ¹∞

and po ðx; y; z; t; xs ; ys ; zs Þ ¼ ð2pÞ

¹1

∞ ¹∞

Po ðx; y; z; q; xs ; ys ; zs Þ expfiqtgdq;

ð16bÞ

while the spatial Fourier transform and its inverse are ∞ ∞ ˜ Po ðx; y; z; q; xs ; ys ; zs Þ expf¹iðkx x þ ky yÞgdxdy Po ðkx ; ky ; z; q; xs ; ys ; zs Þ ¼ ¹∞ ¹∞

ð17aÞ and Po ðx; y; z; q; xs ; ys ; zs Þ ¼ ð2pÞ¹2

∞ ∞ ¹∞ ¹∞

P˜ o ðkx ; ky ; z; q; xs ; ys ; zs Þ

× expfiðkx x þ ky yÞgdkx dky :

ð17bÞ

The change from lower case to upper case denotes the transform from time to

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 857

frequency, and the tilde indicates the transform from space to wavenumber. Using (12), the effects of the transform can be considered separately on the direct wave and the reflected wave. Equation (13) transforms to ¹k2x

¹

k2y

n X ∂2 q2 ˜ D þ 2 þ 2 P0 ¼ ¹ Sk ; ∂z c k¼1

ð18Þ

in which Sk ¼ Qk ðqÞ: expf¹i½kx ðxs þ xgk Þ þ ky ðys þ ygk Þÿg:dðz ¹ zgk Þ;

ð19Þ

while (15) transforms to ¹k2x

¹

k2y

∂2 q2 ˜ R þ 2 þ 2 P0 ¼ 0: ∂z c

ð20Þ

The vertical wavenumber kz is now recognized as 2 1 2 q 2 2 ¹ kx ¹ ky ; kz ¼ 2 c

ð21Þ

and the solutions of (18) and (19) are, respectively, n i X S˜ ðzÞ; P˜ 0D ðkx ; ky ; z; q; xs ; ys ; zs Þ ¼ 2kz k¼1 k

ð22Þ

in which S˜ k ðzÞ ¼ Qk ðqÞ: expf¹i½kx ðxs þ xgk Þ þ ky ðys þ ygk Þÿg: expfikz | z ¹ zgk |g

ð23Þ

and P˜ 0R ðkx ; ky ; z; q; xs ; ys ; zs Þ ¼ P˜ 0R ðkx ; ky ; z ¼ 0; q; xs ; ys ; zs Þ expfikz zg:

ð24Þ

The solution (22) is simply the superposition of the plane-wave decomposition of all the monopole sources, the plane-wave decomposition of a single monopole being given by Aki and Richards (1980, p. 197). The i/2kz factor on the right-hand side of (22) is the well-known ‘obliquity factor’ that arises in the plane-wave decomposition of a spherical wave. The solution (24) states that the scattered field can be only the downgoing wave reflected from the sea-surface, since no upgoing scattered waves can exist. At the sea-surface (12) applies. Therefore, combining the direct and reflected waves yields the incident field in the frequency–wavenumber domain: n i X D S˜ ðzÞ:ð1 ¹ expfikD P˜ 0 ðkD x ; ky ; zw ; q; xs ; ys ; zs Þ ¼ z 2zk gÞ 2kz k¼1 k

as required.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

ð25Þ

858

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

The scattered field: the total measured field minus the incident field The sea-floor is now included and the pressure field in the water, called the total field, now also includes the scattered field, i.e. the waves trapped between the sea-floor and the sea-surface. The pressure in the water at (x, y, z), created by the array of monopole sources centred on the shot position at (xs, ys, zs), is p(x, y, z, t; xs, ys, zs), which is the sum of the incident and scattered pressure fields. The scattered field is therefore the total field minus the incident field: pSCAT ðx; y; z; t; xs ; ys ; zs Þ ¼ pðx; y; z; t; xs ; ys ; zs Þ ¹ p0 ðx; y; z; t; xs ; ys ; zs Þ:

ð26Þ

In this equation the incident field p0(x, y, z, t; xs, ys, zs) is computed from independent measurements made near the source, as described above. The scattered field may therefore be obtained from the measured total field by subtracting the known incident field. We have not fully investigated the problems of doing this subtraction. The important point is that the subtraction may be done if the source wavefield is properly measured.

Response to a single downgoing plane wave and its decomposition into plane-wave components The incident field from a single shot is not a plane wave. In the formulation of the problem given by (6) and (7), it is necessary to obtain the scattered field response to a single plane wave. We begin with pSCAT(x, y, z, t; xs, ys, zs), the scattered field at the point (x, y, z), due to the source at position (xs, ys, zs). This is obtained from the shot gathers by subtracting the incident field, as described above. We now re-order the data into common-receiver gathers, that is, keeping the receiver point (x, y, z) fixed, and then transform over time and over horizontal shot coordinates (xs, ys) to the frequency– wavenumber domain. The spatial transform must use the opposite sign in the exponential from that used in (17a), because the waves are travelling from the transform variable coordinates (xs, ys) in this case, whereas they are travelling to the transform variable coordinates (x, y) in (17a). Thus the triple Fourier transform of the common-receiver gather is ∞ ∞ ∞ D ; k ; z Þ ¼ pSCAT ðx; y; z; t; xs ; ys ; zs Þ Pˆ SCAT ðx; y; z; q; kD x y s ¹∞ ¹∞ ¹∞

D × expfiðqt þ kD x xs þ ky ys Þgdt dxs dys ;

ˆ SCAT

kD x,

ð27Þ

kD y,

(x, y, z; q, zs) is the scattered pressure response in the frequency in which P domain at receiver position (x, y, z) to a plane wave with vertical wavenumber r q2 D 2 D 2 ¹ ðkD kz ¼ ð28Þ x Þ ¹ ðky Þ : c2 This transformation is applied to all receiver gathers. Then, for each frequency component, a further transformation over the horizontal receiver coordinates allows

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 859

this response to be decomposed into individual plane-wave components: ∞ ∞ SCAT ˜ D D D ˆ ðkx ; ky ; z; q; kx ; ky ; zs Þ ¼ Pˆ SCAT ðx; y; z; q; kD P x ; ky ; zs Þ ¹∞ ¹∞

× expf¹iðkx x þ ky yÞgdx dy:

ð29Þ

Separation of upgoing and downgoing waves: the boundary condition at the sea-surface The total pressure field is zero at the sea-surface and, as we have seen, the source and its ghost reduce the incident field to zero at the sea-surface. It follows that the scattered field is zero at the sea-surface. The resulting scattered field we have created in (29) is a superposition of scattered fields each of which is zero at the sea-surface and each of D which satisfies a homogeneous wave equation. Therefore P˜ˆ SCAT(kx, ky, z; q, kD x , ky , zs) also satisfies the homogeneous equation ∂2 q2 ˜ˆ SCAT 2 2 D ðkx ; ky ; z; q; kD ð30Þ ¹kx ¹ ky þ 2 þ 2 P x ; ky ; zs Þ ¼ 0; ∂z c which has the solution SCAT ˜ˆ U D D D ðkx ; ky ; z; q; kD P˜ˆ x ; ky ; zs Þ ¼ P ðkx ; ky ; z ¼ 0; q; kx ; ky ; zs Þ expf¹ikz zg D D þ P˜ˆ ðkx ; ky ; z ¼ 0; q; kD x ; ky ; zs Þ expfþikz zg;

ð31Þ

in which the first term on the right-hand side is the upgoing component of the scattered pressure field and the second term is the downgoing component. Since the scattered field is zero at the sea-surface, the upgoing and downgoing wave components are equal and opposite at the sea-surface. Therefore, from (31), the upgoing and downgoing scattered plane waves in response to the single incident plane-wave component with D coordinates (q, kD x , ky ) may be written as SCAT D ðkx ; ky ; z; q; kD P˜ˆ U x ; ky ; zs Þ expfþikz zg D ; k ; z Þ ¼ P˜ˆ ðkx ; ky ; z ¼ 0; q; kD x y s 1 ¹ expfþikz 2zg

ð32Þ

and SCAT D ðkx ; ky ; z; q; kD P˜ˆ D x ; ky ; zs Þ expfþikz zg ˜ D D ˆ : P ðkx ; ky ; z ¼ 0; q; kx ; ky ; zs Þ ¼ ¹ 1 ¹ expfþikz 2zg

ð33Þ

Equations (32) and (33) simply describe deconvolution for the receiver ghost to perform the up/down separation. It is not necessary to do the transform over source coordinates first to achieve the separation. Our motive in performing the transform over source coordinates first is to separate the scattered field responses for different incident plane waves. Whichever order is chosen, it appears to be necessary to perform five Fourier transforms to separate the data into plane-wave components.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

860

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

In (32) and (33) the denominator can become zero. In this and in all other divisions in frequency–wavenumber space we have used in our implementation well-known stabilization methods such as those proposed by Deregowski (1971, 1978). Computation of the reflection response: the boundary condition at the horizontal plane z ¼ zw Now consider the reflection response in the water, below the source and the receivers, at some depth z ¼ zw, above the sea-floor. We choose an upgoing vertical wavenumber kU z . The upgoing scattered wave at this depth is known from (32) and is U ð34Þ U˜ˆ ðkU ; kU ; z ; q; kD ; kD ; z Þ ¼ P˜ˆ ðkU ; kU ; z ¼ 0; q; kD ; kD ; z Þ expf¹ikU z g: x

w

y

x

s

y

x

y

x

s

y

z

w

The downgoing wave at depth level zw consists of many parts: an incident field component and many scattered field components. Each scattered field component of the downgoing wave is known from (33) and is ˜ˆ ; k ; z ; q; kD ; kD ; z Þ ¼ P˜ˆ D ðk ; k ; z ¼ 0; q; kD ; kD ; z Þ expfþik z g: ð35Þ Dðk x

y

w

x

y

s

x

y

x

y

s

z w

We now need to compute the downgoing incident field at depth z ¼ zw. The contribution to this wave from one shot is given by (25), which, for the chosen downgoing wavenumber, may be written as n i X D ; k ; z ; q; x ; y ; z Þ ¼ ð36Þ S˜ ðz Þ:ð1 ¹ expfikD P˜ 0 ðkD w s s s x y z 2zk gÞ; 2kz k¼1 k w with g g g D D S˜ k ðzw Þ ¼ Qk ðqÞ: expf¹i½kD x ðxs þ xk Þ þ ky ðys þ yk Þÿg: expfikz ðzw ¹ zk Þg:

ð37Þ

The creation of the incident plane wave from all the shots is ∞ ∞ D D D D D ˜ P˜ 0 ðkD P0 ðkx ; ky ; zw ; q; zs Þ ¼ x ; ky ; zw ; q; xs ; ys ; zs Þ expfiðkx xs þ ky ys Þÿgdxs dys : ¹∞ ¹∞

ð38Þ The chosen upgoing scattered plane wave is the sum of the selected reflections of all the downgoing waves: the incident plane wave and all the scattered downgoing waves. Thus, using the boundary condition of (7) the upgoing wave can be expressed as ˜ˆ U U ˜ˆ D D U D D D D U˜ˆ ðkU x ; ky ; zw ; q; kx ; ky ; zs Þ ¼ Rðkx ; ky ; zw ; q; kx ; ky ; zs ÞP 0 ðkx ; ky ; zw ; q; zs Þ ∞ ∞ ˜ˆ ; k ; z ; q; kD ; kD ; z Þdk dk : ð39Þ ˜ˆ U ; kU ; z ; q; k ; k ; z Þ:Dðk þ ð2pÞ¹2 Rðk w x y s x y w s x y x y x y ¹∞ ¹∞

Recovery of the reflection coefficients In practice, the double integral in (39) is a summation over all wavenumbers. In fact, all the integrals in this paper must be handled as summations. At a given frequency q let

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 861

the discretization of shots and receivers give N independent plane-wave components, each defined by a (kx, ky) pair, and let us number these 1, 2, . . . , N. The incident wave is, say, the ith one of these, and the chosen upgoing wave is, say, the jth one. Then, following (6), (39) could be written as Uji ¼ Rji Pi þ

N X

Rjk Dki ;

ð40Þ

k¼1

which is one equation with N unknowns, Rjk, k ¼ 1, 2, . . . , N. The whole procedure must therefore be carried out for every possible incident plane wave, giving N simultaneous linear equations with N unknowns: Uji ¼ Rji Pi þ

N X

Rjk Dki ;

i ¼ 1; 2; . . . ; N:

ð41Þ

k¼1

These may then be solved for the Rjk using normal methods. Since this is valid for a single upgoing plane wave and for a single frequency, it must be repeated for all upgoing plane-wave components and for all frequencies.

Re-introduction of an incident field and transformation back to the space– time domain In order to look at the response in the space–time domain we now re-introduce an incident field, multiply by the reflection coefficient, and transform back to the space– time domain. There is some freedom in the choice of incident field. The directivity pattern of the source array as well as the source ghost is taken into account in the computation of the true incident field and both these factors have been removed, essentially deconvolved, from the data to recover the reflection coefficients. The reintroduction of the incident field permits a less directional source with a shorter source time function to be used. The new incident field is thus i ¹ D D D Z1 ðqÞ:Z2 ðkD P˜ˆ 0 ðzw ; kD x ; ky ; q; z0 Þ ¼ x Þ:Z3 ðky Þ: expfikz ðzw ¹ z0 Þg; 2kD z

ð42Þ

in which the superscript minus sign is introduced to show that the free surface has been removed. Ideally, the functions Z1, Z2 and Z3 would all be constants, thus placing the same delta-function monopole source at points at depth z0 below sea-level. (The seasurface is now absent.) In practice, because of the need to control the influence of additive noise, jZ1(q)j must be similar to the sum of the original notional source spectra plus their ghosts (Ziolkowski, Underhill and Johnston 1998), and it must be expected that the directivity of the source array cannot be removed entirely, thus forcing Z2(kD x ) ) to be non-constant. In any case, it is at this point that source signature and Z3(kD y deconvolution and removal of source directivity take place. The depth z0 of the chosen source is arbitrary, but it must be no greater than zw. A single plane-wave component of the reflection response of the earth at level z in

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

862

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

response to the chosen single incident plane wave is ¹ ˜ˆ U U U D D U D D P˜ˆ ðkU x ; ky ; z; kx ; ky ; q; z0 Þ ¼ expfikz ðzw ¹ zÞg:Rðkx ; ky ; zw ; kx ; ky ; qÞ ¹ D × :P˜ˆ 0 ðzw ; kD x ; ky ; q; z0 Þ;

ð43Þ

in which the superscript minus sign indicates the absence of the free surface. The multiple-free seismograms in the space–time domain are then obtained by the fivedimensional inverse Fourier transform: ∞ ∞ ∞ ∞ ∞ 1 U D D P˜ˆ ¹m ðkU p¹m ðx; y; z; t; xs ; ys ; z0 Þ ¼ x ; ky ; z; kx ; ky ; q; z0 Þ: ð2pÞ5 ¹∞ ¹∞ ¹∞ ¹∞ ¹∞ D U U × expf¹iðqt þ kD x xs þ ky ys ¹ kx x ¹ ky yÞg D U U × dkD x dky dkx dky dq

ð44Þ

This is the desired result. A summary of the processing steps The steps in the processing of the data may be summarized as follows: 1 Remove the incident field from the measurements in the space–time domain. This can be done if the incident field is measured at the source. This leaves the scattered field. 2 Sort the data into common-receiver gathers. D 3 Transform every gather into the frequency–wavenumber (q, kD x , ky ) domain with a triple Fourier transform. This gives the response at each receiver to every possible incident plane wave. D 4 Re-sort the data into incident plane-wave gathers; that is, for each (q, kD x , ky ) component, collect the response at every receiver. This is a two-dimensional spatial array for each frequency and for each of the N possible incident waves. 5 Transform each of these arrays to the wavenumber domain with a double Fourier transform, giving N wavenumbers for each array. 6 Separate each component into an upgoing part and a downgoing part using the zeropressure condition at the water-surface. 7 Each upgoing plane wave is then the sum of the reflection of the incident plane wave and the reflections of all the downgoing scattered plane-wave components, as given by (40). Construct equations (41) for the same upgoing plane-wave component for each incident plane wave for one frequency. Solve for the reflection coefficients. 8 Repeat step 6 for all upgoing wavenumbers and for all frequencies. ¹ D 9 Introduce a desired incident wavefield P˜ˆ 0 (zw, kD x , ky , q, z0). 10 Construct the plane-wave reflection responses for this incident field using (43) and transform back to the space–time domain with a five-dimensional Fourier transform over upgoing wavenumbers, all downgoing wavenumbers, and over all frequencies.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 863

An example for a one-dimensional earth model Equations (41) may be written in the following matrix form: 3 2 32 3 2 Rj1 Uj1 P1 þ D11 D12 : D1N 76 R 7 6U 7 6 D P2 þ D22 : D2N 21 76 j2 7 6 j2 7 6 7 6 76 7 6 76 : 7 6 : 7 6 : : 7¼6 76 7 6 76 : 7: 6 : 7 6 : : : : 7 6 76 7 6 7 6 76 7 6 54 : 5 4 : 5 4 : : UjN

DN1

DN2

:

PN þ DNN

ð45Þ

RjN

For a one-dimensional earth, the incident wave, the upgoing scattered wave and the downgoing scattered wave all have the same wavenumber. There are no off-diagonal terms and we simply have Ui ¼ ðPi þ Di ÞRi ;

ð46Þ

which is the same as (5), as expected. Equation (46) may be rearranged to give Ri ¼

Ui ; Pi þ Di

ð47Þ

or, using our frequency–wavenumber notation, as ˜ x ; ky ; zw ; qÞ ¼ Rðk

U˜ ðkx ; ky ; zw ; qÞ ; ˜ x ; ky ; zw ; qÞ P˜ 0 ðkx ; ky ; zw ; qÞ þ Dðk

ð48Þ

in which P˜0 is given by (25), which incorporates all the source signature information, and only a single-shot gather is required because they are all the same. For a monopole source it can be shown that this formulation is consistent with the forward problem ˆ is broad formulated by Fokkema and Ziolkowski (1987). The reflection response R bandwidth, signature deconvolution having been handled by the division in (49). Figure 10 shows the results of a test of this approach. The model is the same as in Fig. 2a. Figure 10a is the same synthetic seismogram as Fig. 4b; Fig. 10b shows the result of processing the synthetic data of Fig. 10a to remove the multiples. Figure 10c is the same as Fig. 4d: the response without a water-surface; Fig. 10d is the difference between Figs 10b and c and shows the numerical noise. Discussion of 2D split-spread measurements and source problems The ideal configuration of shots and receivers is as shown in Fig. 11. In practice 3D marine seismic data are obtained with a vessel towing one or more source arrays and multiple receiver cables. Normally the source is in front of the receivers. To permit the proposed scheme to be applied in practice, it would be necessary to tow the source at a considerable distance behind the vessel to permit the data to be obtained in a twodimensional split-spread configuration. This is feasible, but has not yet been done, as

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

864

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Figure 10. Test of method on synthetic data from a one-dimensional earth model. (a) Same seismogram as Fig. 4b; (b) result of processing data in Fig. 10a to remove the free surface; (c) same seismogram as Fig. 4d; (d) difference between (b) and (c) to show the numerical noise.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 865

Figure 11. Ideal layout (schematic) of sources (circles) and receivers (triangles). The shot interval is the same in both horizontal directions and equal to the receiver interval, also the same in both directions, and equal to half the shortest wavelength of interest to prevent spatial aliasing. Clearly this layout is easier to achieve today with sea-floor receivers than with towed streamers.

far as we are aware. A more expensive alternative is to use two vessels, one in front of the other, the one in front towing hydrophone cables and the one behind towing the source and hydrophone cables. With conventional cables the in-line sampling interval is much smaller than the cross-line interval. There are therefore interpolation problems cross-line. These are simply problems that have to be faced. The methods of Berkhout–Verschuur, Carvalho–Weglein–Stolt and Fokkema–van den Berg–van Borselen face exactly the same problems. It is important to have the correct sampling both in-line and cross-line, although, clearly, this is currently difficult to achieve. The vessel is assumed to sail in parallel lines, first in one direction and then in the opposite direction. The sail lines must be close together to obtain adequate spatial sampling for the common-receiver gathers (step 2 above). Also, to ensure that the source directivity remains the same when the vessel changes direction, the source array must have both in-line and cross-line symmetry. This method also assumes excellent shot-to-shot signature repeatability, as do all the other wave-theoretical methods.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

866

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

It is clearly not practical to implement this method at the time of writing. However, the method may easily be adapted to ocean-bottom cable data (White 1965; Barr and Sanders 1989), in which hydrophones and vertical geophones are used to obtain the separation of upgoing and downgoing scattered waves without any loss of bandwidth introduced by the receiver ghost. The layout of a 2D patch of 2-component or 4-component (3-component geophone plus a hydrophone) receivers is feasible at present, and the shooting boat can easily be designed to perform the simple tasks described here. The only remaining practical issue is shot-to-shot repeatability, which is a requirement for the superposition of sources in common-receiver gathers to simulate incident plane waves. Airguns are now very reliable and can be synchronized with sufficient precision. The major cause of shot-to-shot signature variations is variations in the depth at which the guns are towed. Most seismic contractors tow airgun arrays that are suspended from buoys or floats, as shown in Fig. 12a. The drag on the float is equal but not opposite to the force in the towing cable, and there is thus a couple acting on the system. Variations in this couple, caused by variations in the drag on the float, result in variations in the depth of the guns. Many years ago, T.-A. Haugland proposed towing the float rather than the guns themselves, as shown in Fig. 12b, ensuring that the force in the towing cable is equal and opposite to the drag on the float, and therefore minimizing the main cause of depth variations in the source. This system was implemented in 1981 on Seismic Profilers’ two vessels the Nina Profiler and the Liv Profiler (Seismic Profilers A/S 1982; Ziolkowski et al. 1982). It should be noted that proper spatial sampling of the data would give enormously high fold of coverage. A smaller signal-to-noise ratio of a single-shot record could therefore be tolerated, allowing the source array to be smaller and more compact, becoming a better approximation to a point source. The small source would be easier to control and, with fewer source elements, would be more reliable, giving better shotto-shot repeatability. Conclusions A new wave-theoretical method is proposed for removing multiples from marine 3D seismic reflection data. It is fundamentally a prestack method. Using measurements of the source wavefield made during data acquisition, three problems associated with earlier work are solved: (i) the new method can accommodate source arrays, rather than point sources; (ii) it can remove the incident field without simultaneously removing part of the scattered field; (iii) the minimum-energy criterion (to find a wavelet) is eliminated. The essence of the method is to obtain the reflection response at a plane above the sea-floor as the ratio of the upgoing wave to the downgoing wave in the frequency– wavenumber domain. For a given incident downgoing plane wave there are many downgoing scattered waves. In practice therefore this principle leads to a set of

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 867

(a)

VARIABLE DRAG

TOW

DRAG

(b)

TOW

VARIABLE DRAG

UNIFORM DRAG

Figure 12. (a) Airguns are towed while suspended from a float or paravane; (b) airguns are suspended from a float which is towed.

simultaneous equations for the reflection coefficients for each upgoing plane-wave component and for each frequency. The unknowns in our formulation are these reflection coefficients. The known constituents are the plane-wave components of the incident field, the upgoing scattered field and the downgoing scattered field. These components can be separated from the measurements, provided the source wavefield is measured and the data are properly sampled in time and space. To create the shot records that would have been measured if the sea-surface had been absent, each reflection coefficient is multiplied by complex amplitude and phase factors, for source and receiver terms, before the five-dimensional Fourier transformation back to the space–time domain. The source wavefield must be known in detail, as the signature, including the source

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

868

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

ghost, is required for every downgoing wavenumber. The directional response of the source array can be deconvolved with this approach. We have shown that the method works for the simple case of a line source in a water layer overlying an elastic half-space. However, this is only a partial validation of our method. Taylor has validated our method analytically for the 2D case for an arbitrary distribution of line diffractors below the source (Fast 2-D synthetic seismograms to test sea-surface multiple removal algorithms, submitted to Geophysical Prospecting). Taylor’s proof can readily be generalized to 3D. Taylor’s new algorithm can generate fast 2D synthetic seismograms on which we are currently testing the method. Twodimensional validation is required before the 3D problem can be tackled. As it stands, the method is unlikely to be implemented on towed streamer data for several years, for practical reasons. However, it may be applied now to arrays of sea-floor 2-component and 4-component data. Acknowledgements This paper has resulted after many false starts over the last six years. We thank our colleagues, particularly Sergei Zatsepin, for critically reading the manuscript and finding flaws in our previous attempts to find a solution. Anton Ziolkowski thanks Lasse Amundsen for an important conversation on 6 September 1996 (and for several subsequent discussions) and Roger Bilham and Clint Frasier who allowed him to think aloud the morning after the night before during the Denver SEG meeting. This paper has been reviewed by four reviewers, two editors, and two Associate Editors; we thank them all for their help and advice which has resulted in two major revisions and a minor revision. We hope the paper is now much clearer. This research has been partly funded by the Petroleum Science and Technology Institute through their funding of the Chair of Petroleum Geoscience, partly by the Natural Environment Research Council via project no. GR3/9064, partly by the Engineering and Physical Science Research Council via ROPA award no. GR/L73524, and partly through the financial support of Broken Hill Proprietary, Geco-Prakla (UK) Limited, Saga Petroleum A/S, and Texaco Britain Limited.

References Aki K. and Richards P.G. 1980. Quantitative Seismology. W.H. Freeman. Barr F.J. and Sanders J.I. 1989. Attenuation of water-column reverberations using pressure and velocity detectors in a water-bottom cable. 59th SEG meeting, Dallas, USA, Expanded Abstracts, 653–656. Berkhout A.J. 1982. Seismic Migration: Imaging of Acoustic Energy by Wavefield Extrapolation. A. Theoretical Aspects, 2nd edn, pp. 211–218. Elsevier Science Publishing Co. van Borselen R.G., Fokkema J.T. and van den Berg P.M. 1996. Removal of surface-related wave phenomena – the marine case. Geophysics 61, 202–210. Carvalho P.M. and Weglein A.B. 1994. Wavelet estimation for surface multiple attenuation using

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

Marine seismic wavefield measurement 869 a simulated annealing algorithm. 64th SEG meeting, Los Angeles, USA, Expanded Abstracts, 1481–1484. Carvalho P.M., Weglein A.B. and Stolt R.H. 1991. Examples of a nonlinear inversion method based on the T matrix of scattering theory: application to multiple suppression. 61st SEG meeting, Houston, USA, Expanded Abstracts, 1319–1322. Carvalho P.M., Weglein A.B. and Stolt R.H. 1992. Nonlinear inverse scattering for multiple suppression: application to real data. Part 1. 62nd SEG meeting, New Orleans, USA, Expanded Abstracts, 1093–1095. Deregowski S.M. 1971. Optimal digital filtering and inverse filtering in the frequency domain. Geophysical Prospecting 19, 729–768. Deregowski S.M. 1978. Self-matching deconvolution in the frequency domain. Geophysical Prospecting 26, 252–290. Fokkema J.T. and van den Berg P.M. 1990. Removal of surface-related wave phenomena: the marine case. 60th SEG meeting, San Francisco, USA, Expanded Abstracts, 1689–1692. Fokkema J.T. and van den Berg P.M. 1993. Seismic Applications of Acoustic Reciprocity. Elsevier Science Publishing Co. Fokkema J.T. and Ziolkowski A.M. 1987. The critical reflection theorem. Geophysics 52, 965–972. Fuchs K. and Mu¨ller G. 1971. Computation of synthetic seismograms with the reflectivity method and comparison with observations. Geophysical Journal of the Royal Astronomical Society 23, 417–433. Ikelle L.T., Roberts G. and Weglein A.B. 1997. Source signature estimation based on the removal of first-order multiples. Geophysics 62, 1904–1920. Kennett B.L.N. 1979. The suppression of surface multiples on seismic records. Geophysical Prospecting 27, 584–600. Mayne W.H. 1950. Seismic surveying. US Patent 2 732 906 (application 1950). Mayne W.H. 1962. Common reflection point horizontal data stacking techniques. Geophysics 27, 927–938. Parkes G.E., Ziolkowski A.M., Hatton L. and Haugland T. 1984. The signature of an airgun array: computation from near-field measurements including interactions – practical considerations. Geophysics 49, 105–111. Peacock K.L. and Treitel S. 1969. Predictive deconvolution: theory and practice. Geophysics 26, 754–760. Pieprzak A.W., Henman R.E., Saltaug R., Moldoveanu N., Parker J.A. and Vlasin J. 1998. Real time source signature estimator – a case study of performance in marginal weather from the Orca Basin experiment. 68th SEG meeting, New Orleans, USA, Expanded Abstracts, 111– 114. Riley D.C. and Claerbout J.F. 1976. 2-D multiple reflections. Geophysics 41, 592–620. Seismic Profilers A/S 1982. Areal Airgun Arrays: Brochure of Seismic Profilers A/S, Lokketangen 20, 1300 Sandvika, Norway. Taner M.T., O’Doherty R.F. and Koehler F. 1995. Long period multiple suppression by predictive deconvolution in the x–t domain. Geophysical Prospecting 43, 433–468. Tatham R.H. 1989. Tau-p filtering. In: Tau-p, a Plane Wave Approach to the Analysis of Seismic Data (ed. P.L. Stoffa). Kluwer Academic Publishers. Tatham R.H., Keeney J.W. and Noponen I. 1983. Application of the t-p transform (slant-stack) in processing seismic reflection data. Australian Society of Exploration Geophysicists 14, 163– 172.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

870

A.M. Ziolkowski, D.B. Taylor and R.G.K. Johnston

Vermeer G.O. 1990. Seismic Wavefield Sampling. Society of Exploration Geophysics. Verschuur D.J. 1991. Surface-related multiple elimination, an inversion approach. Doctoral dissertation, Delft University of Technology. Verschuur D.J. and Berkhout A.J. 1992. Surface-related multiple elimination: practical aspects. 62nd SEG meeting, New Orleans, USA, Expanded Abstracts, 1100–1103. Verschuur D.J., Berkhout A.J. and Wapenaar C.P.A. 1989. Wavelet estimation by prestack multiple elimination. 59th SEG meeting, Dallas, USA, Expanded Abstracts, 1129–1132. Verschuur D.J., Berkhout A.J. and Wapenaar C.P.A. 1992. Adaptive surface-related multiple elimination. Geophysics 57, 1166–1177. Weglein A.B., Gasparotto F.A., Carvalho P.M. and Stolt R.M. 1997. An inverse scattering series method for attenuating multiples in seismic reflection data. Geophysics 62, 1975–1989. White J.E. 1965. Seismic Waves: Radiation, Transmission and Attenuation. McGraw-Hill Book Co. Ziolkowski A. 1998. The measurement of airgun bubble oscillations. Geophysics 63, 2009–2024. Ziolkowski A.M. and Johnston R.G.K. 1997. Marine seismic sources: QC of wavefield computation from near-field pressure measurements. Geophysical Prospecting 45, 611–639. Ziolkowski A.M., Parkes G.E., Hatton L. and Haugland T.-A. 1982. The signature of an airgun array: computation from near field measurements including interactions. Geophysics 47, 1413–1421. Ziolkowski A.M., Underhill J.R. and Johnston R.G.K. 1998. Wavelets, well-ties, and the search for subtle stratigraphic traps. Geophysics 63, 297–313.

q 1999 European Association of Geoscientists & Engineers, Geophysical Prospecting, 47, 841–870

ZAPDOC.TIPS | To ensure the functioning of the site, we use **cookies**. We share information about your activities on the site with our partners and Google partners: social networks and companies engaged in advertising and web analytics. For more information, see the Privacy Policy and Google Privacy & Terms. Your consent to our cookies if you continue to use this website. Accept