Method of processing seismic data

[0001] 1. Field of the Invention

[0002] This present invention relates to a method of processing seismic data and, in particular, to a method of processing marine seismic data to, reduce the effects of “ghost reflections”. The present invention can be applied to processing both marine seismic data acquired in a flat sea and marine seismic data acquired in a rough sea. The present invention also relates to processing seismic data so as to estimate the depth of a seismic receiver, or so as to monitor the accuracy of a seismic receiver.

[0003] 2. Description of the Related Art

[0004] FIG. 1 is a schematic diagram of a marine seismic survey in which seismic energy is emitted from a source 1 and detected by a seismic receiver 2 at a depth h below the surface 6 of the sea. Energy emitted from the source is reflected by the seabed 3 or by a reflector 4 below the seabed 3 and is then detected by the receiver. This path of seismic energy is labelled 5 in FIG. 1. Information about the geological structure of the earth's interior can be derived from the reflected seismic energy incident on the receiver.

[0005] The seismic receiver 2 shown in FIG. 1 is a streamer, which is a type of seismic receiver that is frequently used in marine seismic surveying. A streamer contains a plurality of sensors S₁, S₂ . . . S_N such as pressure sensors and/or particle velocity sensors distributed along its length, which may be some hundreds of meters, and is thus able to measure the reflected seismic energy at a number of points simultaneously. A streamer is generally suspended from one or more floats 8, so that all the receivers of the streamer are at the same depth in a flat sea.

[0006] In addition to the desired path 5 of seismic energy shown in FIG. 1, other seismic energy paths will occur as a result of seismic energy being reflected or scattered from the sea surface 6. These paths are known as “ghost reflections”. For example, reference 7 in FIG. 1 shows a ghost reflection in which the seismic energy reflected by the reflector 4 is not directly incident on the receiver 2, but undergoes a further reflection at the sea surface 6 before reaching the receiver. Down-going sea-surface ghost reflections are an undesirable source of contamination of seismic data, since they obscure the interpretation of the desired up-going reflections from the earth's interior.

[0007] Ghost reflections also arise from sea surface reflections or scattering near the source 1. The overall seismic energy signal received at a point below the source will contain a component that has come direct from the source and a component that has come via a reflection or scattering at the sea surface.

[0008] Removing the ghost reflections from seismic data is for many experimental configurations equivalent to up/down wavefield separation of the recorded data. In such configurations the down-going part of the wavefield represents the ghost signal and the up-going wavefield represents the desired signal. The up/down wavefield separation can be achieved by applying a suitable filter to the seismic data.

[0009] The present invention provides a method of de-ghosting seismic data, the method comprising the steps of:

[0010] a) obtaining seismic data using a seismic source and a seismic receiver disposed within a fluid medium;

[0011] b) determining a derivative with respect to time of the fluid pressure at a point in the vicinity of a selected one of the seismic source and the receiver;

[0012] c) determining a derivative with respect to a first direction of the fluid pressure at a point in the vicinity of the selected one of the seismic source and receiver;

[0013] d) determining a derivative with respect to a second direction of the fluid pressure at a point in the vicinity of the selected one of the seismic source and receiver from the results of steps (b) and (c); and

[0014] e) processing the seismic data using the derivative with respect to a second direction of the fluid pressure, thereby to reduce the effects in the processed seismic data of seismic energy reflected and/or scattered at the sea surface.

[0015] In a typical marine seismic survey, the time-derivatives of the pressure and the derivatives in one direction, for example the horizontal direction, are readily available, whereas the derivatives in other directions, for example the vertical direction, are not directly known and have to be estimated. The present invention uses the readily available derivatives to provide a better estimate of a derivative that is not directly known, and this allows more accurate de-ghosting of the seismic data.

[0016] In a preferred embodiment step (a) comprises obtaining seismic data at two different depths below the surface of the fluid medium, and in a particularly preferred embodiment step (a) comprises measuring the fluid pressure at two different depths below the surface of the fluid medium. This allows the derivative of the de-ghosting step to be further improved.

[0017] In an alternative preferred embodiment, step (a) comprises measuring the fluid pressure and the particle velocity. If the depth is known, this embodiment allows the pressure and particle velocity sensors to be monitored. For example, the particle velocity can be computed from the depth and the measured pressure values, and the computed value may then be compared with the measured values of the particle velocity. Alternatively, the depth may be estimated from the measured values of the particle velocity and pressure.

[0018] A second aspect of the present invention provides a method of calibrating a seismic receiver as claimed in claim 20. This method enables the accuracy of measurement of particle velocity to be monitored.

[0019] A third aspect of the present invention provides a method of calibrating a seismic receiver as claimed in claim 22, and a fourth aspect of the invention relates to a method of determining the depth of a seismic receiver as claimed in claim 24. These aspects of the invention relate to a receiver that can measure both the particle velocity and the fluid pressure.

[0020] Other features of the invention are set out in the dependent claims.

[0034] The equation relating pressure p to particle velocity ν n an acoustic medium with a bulk modulus K is:
{dot over (p)}=−K{overscore (∇)}.{overscore (ν)}  (1)

[0035] In equation (1), {dot over (p)} denotes the time-derivative of p. The pressure p and particle velocity v are also related by the equation of motion, which has the following form

[0036] $\begin{array}{cc}\stackrel{\_}{v}=-\frac{1}{\rho }\nabla p& \left(2\right)\end{array}$

[0037] Here, p is the density of the acoustic medium. The equation of motion (2) and equation (1) together complete the set of first order partial differential equations describing acoustic wave propagation. These equations can be combined to give the well-known, second-order, acoustic wave equation:
{umlaut over (p)}=α²∇²p  (3)

[0038] In equation (3), α=√{square root over (K/p)} is the wave propagation speed in the acoustic medium.

[0039] Since the curl of the gradient of an acoustic wavefield is zero (this being a fundamental vector relation), equation (2) yields:
∇×{overscore (ν)}={overscore (0)}  (4)

[0040] A preferred embodiment of the invention relating to de-ghosting a streamer located in the vicinity of the sea-surface will now be described. In this embodiment the streamer is provided with a plurality of pressure sensors distributed along its length. Although the propagation of seismic energy is three-dimensional in any real seismic surveying arrangement, this example will consider only 2D wave propagation in the x-z plane (z increases downwards, as indicated in FIG. 1) for clarity of explanation. The invention can, however, be applied to full three-dimensional propagation of seismic energy.

[0041] It has been reported, for example in co-pending UK patent application No 9906456.0, that de-ghosting of seismic data can be carried out if the vertical derivative of pressure ∂_zp is known. Thus, it is necessary to compute the vertical derivative of pressure along the streamer at a point (x₀z₀). This can be achieved by writing an estimate of the vertical derivative of pressure as a Taylor expansion in the z-direction of p around (x₀z₀). After rearranging the terms of the Taylor expansion we obtain:

[0042] $\begin{array}{cc}{\partial }_{z}p\left(x_0,z_0\right)=\frac{p\left(x_0,z_0\right)-p\left(x_0,z_0-h\right)}{h}+\frac{h}{2}{\partial }_{\mathrm{zz}}p\left(x_0,z_0\right)+O\left(h^2\right)& \left(5\right)\end{array}$

[0043] In equation (5), h is the vertical distance of the streamer at point (x₀,z₀) below the sea-surface, ∂_z denotes a first-order derivative in the z-direction and ∂_zz denotes a second order derivative. The 0(h²) term shows that the largest error terms have factor h² in them and can therefore be ignored to a reasonably accurate approximation.

[0044] In the case of a marine seismic survey that has a single streamer it is not possible to measure ∂_zzp directly. Pressure is measured only at the streamer, and this may be at a constant depth in a flat sea. However, the second-order derivative of p in the vertical direction, ∂_zzp(x₀z₀), can be related to the time derivative of pressure and to the derivative of pressure in a direction along the streamer, which is assumed to be horizontal, by using the wave equation. This allows us to increase the accuracy of a derivative estimate. A correction of this type is known as a Lax-Wendroff correction in the synthetic finite-difference literature, and is described by Lax and Wendroff in “Difference Schemes for Hyperbolic Equations with high order accuracy”, Comm. Pure Appl. Math Vol. 27 (1964), and by Blanch and Robertsson in “Geophysical Journal International”, Vol. 131, pp381–386 (1997).

[0045] Thus, according to the invention Equation (3) is used to convert derivatives that are not known into derivatives that can be computed from the measured pressure data. From equation (3) the following relationship can be obtained:

[0046] $\begin{array}{cc}{\partial }_{\mathrm{zz}}p\left(x_0,z_0\right)=\frac{1}{{\alpha }^2}\ddot{p}\left(x_0,z_0\right)-{\partial }_{\mathrm{xx}}p\left(x_0,z_0\right)& \left(6\right)\end{array}$

[0047] By substituting the above formula for ∂_zzp, Equation (5) then becomes:

[0048] $\begin{array}{cc}{\partial }_{z}p\left(x_0,z_0\right)=\frac{p\left(x_0,z_0\right)}{h}+\frac{h}{2{\alpha }^2}\ddot{p}\left(x_0,z_0\right)-\frac{h}{2}{\partial }_{\mathrm{xx}}\left(x_0,z_0\right)+O\left(h^2\right)& \left(7\right)\end{array}$

[0049] Here we have made use of the fact that the pressure vanishes at the sea-surface:
p(x₀;z₀−h)=0.

[0050] The significance of equation (7) is it does not contain ∂_zzp. In equation (7), ∂_zp is expressed in terms of the derivative of the pressure with respect to time and the derivative of the pressure with regard to the x-direction (the horizontal direction). Both of these quantities can be readily estimated from the available data. The time-derivative of the pressure can be readily computed from pressure measurements made by a sensor over time. Moreover, as noted above a streamer contains a plurality of pressure sensors disposed along its length, and the derivative of the pressure in the x-direction can be readily computed from pressure measurements made at a given time by adjacent pressure sensors along the streamer.

[0051] Equation (7) also provides a means of computing both the x-component and the z-component of particle velocity along the streamer, using equation (2), provided that the streamer depth at each point along the streamer is known at all times.

[0052] The derivatives with respect to x,∂_xp, in the expression for ν_x derived from equation (2) and ∂_xxp in equation (7), can be computed by differentiating pressure measurements along the streamer. As described in co-pending UK patent application No 9906456.0, computing ∂_zp enables the seismic data to be de-ghosted. Applying the following filter to the seismic data provides de-ghosting, valid for all incidence angles:

[0053] $\begin{array}{cc}\tilde{p}=\frac{1}{2}\left(p+\frac{1}{{\mathrm{ik}}_z}{\partial }_{z}p\right)& \left(8\right)\end{array}$

[0054] The variable {overscore (p)} denotes the recorded up-going pressure (deghosted), i is the square root of −1, k₂=√{square root over ((w/α)²−k_x²)} is the vertical wavenumber, k_x is the horizontal wavenumber, and w is the angular frequency.

[0055] The method for estimating ∂_zp according to the invention in principle is independent of the roughness of the sea, as is the deghosting step described in equation (8).

[0056] Thus, according to the present invention, the derivative in the z (vertical) direction of the pressure is determined from the pressure measurements made along a horizontal line. Once the derivative of the pressure in the z direction has been calculated, the seismic data can be de-ghosted using equation (8).

[0057] FIGS. 3(a) to 6(b) compare results of the present invention with results obtained by a prior art method. These figures are simulations obtained using a flat sea surface or using a rough sea-surface having the profile shown in FIG. 2. It will be seen that the rough sea surface of FIG. 2 has a significant wave height of around 2 m. In the simulation it was assumed that the pressure is being recorded by a horizontal streamer at the points marked with an asterix, namely at x=−12.5 m, x=−137.5 m, x=−262.5 m and x=−387.5 m. The source is a 50 Hz zero-phase explosive source located at x=0 and at a depth of 50 m.

[0058] FIGS. 3(a) and 3(b) show a comparison between ∂_zp calculated explicitly using a finite-difference method (solid line) and an estimate calculated without a Lax-Wendroff correction (broken lines). The traces in each figure correspond to recording locations x=−12.5 m (uppermost trace), x=−137.5 m, x=−262.5 m, and x=−387.5 m (lowest trace).

[0059] In FIGS. 3(a) and 3(b) ∂_zp was calculated simply by dividing the pressure recorded by the sensor on the streamer by the vertical distance between the recording location and the sea-surface (as noted above, the pressure becomes zero at the sea-surface). Thus, this approximation assumes that ∂_zp is constant between the sea-surface and the streamer. FIG. 3(a) shows the results of a streamer depth of 3.5 m, and FIG. 3(b) shows the results for a streamer depth of 5.5 m. It will be seen that computing ∂_zp by the approximation of dividing the recorded pressure at the streamer by the vertical distance between the recording location and the sea-surface is not particularly accurate.

[0060] FIGS. 4(a) and 4(b) again relate to a flat sea and show (in dashed lines) results of a method of the present invention in which ∂_zp was calculated using equation (7) above compared with the results (solid lines) of computing ∂_zp using a finite-difference calculation. It will be seen that this method of the present invention provides a much more accurate estimate of ∂_zp than does the simple geometric averaging used in FIGS. 3(a) and 3(b).

[0061] FIGS. 5(a) to 6(b) show results of computing ∂_zp for the rough-sea shown in FIG. 2. FIGS. 5(a) and 5(b) show results for the prior art method. In these figures, the results in broken lines are for a calculation of ∂_zp by the prior art method of dividing the pressure recorded by the sensor on the streamer by the vertical distance between the recording location and the sea-surface, while the results in solid lines are a finite-difference calculation for ∂_zp.

[0062] FIGS. 6(a) and 6(b) show results obtained by using the present invention. In these Figures, the solid lines again show the result of a finite-difference calculation for ∂_zp, whereas the broken lines show the results of computing ∂_zp using equation (7). It will again be seen that the method of the present invention provides much more accurate estimates of ∂_zp than does the simple geometric averaging method.

[0063] In each of FIGS. 4(a), 4(b), 5(a), 5(b), 6(a) and 6(b), the traces correspond to recording locations x=−12.5 m (uppermost trace), x=−137.5 m, x=−262.5 m, and x=−387.5 m (lowest trace).

[0064] Since the present invention provides more accurate values of ∂_zp, the accuracy of the de-ghosting using the filter of equation (8) will also be improved.

[0065] In a further embodiment of the invention, the accuracy of the estimates of ∂_zp is further increased by eliminating the next higher order term in the Taylor expansion, and increasing the accuracy to 0(h³). This requires solving an implicit system of equations for ∂_zp at each frequency:

[0066] $\begin{array}{cc}{\partial }_z\left(1-\frac{h^2}{6}\left(\frac{w^2}{{\alpha }^2}+{\partial }_{\mathrm{zz}}\right)\right)p\left(x_0,z_0\right)=\frac{p\left(x_0,z_0\right)}{h}-\frac{h{w}^2}{2{\alpha }^2}p\left(x_0,z_0\right)-\frac{h}{2}{\partial }_{\mathrm{xx}}p\left(x_0,z_0\right)+O\left(h^3\right)& \left(9\right)\end{array}$

[0067] One way of solving equation (9) is in the f−x domain using a second order accurate derivative estimate to approximate the term ∂_zz on the left-hand side. It may also be preferable to find filter approximations to equation (9) since equation (9) has effectively introduced a pole that could cause instabilities or amplify noise.

[0068] FIGS. 7(a) and 7(b) shows the results of this higher-order de-ghosting method of the present invention. These figures show the result of calculating ∂_zp using a finite-difference method (solid lines) and the results according to the invention of calculating ∂_zp using equation (9) (broken lines). The four traces in each figure again correspond to recording locations of x=−12.5 m (uppermost trace), x=−137.5 m, x=−262.5 m, and x=−387.5 m (lowest trace). FIGS. 7(a) and 7(b) both relate to a streamer depth of 3.5 m, and FIG. 7(a) shows results for a flat sea and FIG. 7(b) shows the results for the rough sea surface of FIG. 2. It will be seen that estimating ∂_zp using equation (9) produces generally good results.

[0069] The present invention is applicable to a seismic receiver in the vicinity of the sea-surface, such as a streamer that is towed at a shallow depth below the sea-surface. The maximum allowable depth of the streamer below the sea-surface is generally similar to the maximum separation used to separate the two streamers in a twin streamer configuration, since the largest inaccurate term in both cases is of O(h²), or of O(h³) if equation (9) is used. It is preferable that the depth of the streamer below the sea surface does not exceed half the wavelength of seismic energy emitted by the source.

[0070] The present invention requires that the vertical distance below the sea surface of each of the sensors distributed along the streamer must be known. A conventional streamer is generally provided with a depth sensor so that an operator can monitor the depth of the sensor. Provided that this depth sensor is sufficiently accurate, it will be possible to use the average depth of the streamer as measured by this sensor as the depth of each of the recording elements. This should provide good results in flat sea conditions. For improved accuracy in rough sea conditions, however, it is preferable that a plurality of depth sensors, such as hydrostatic sensors, are provided on the streamer, so that the distance of each pressure sensor below the instantaneous sea-surface can be recorded. Alternatively, it is possible that the depth of each sensor can be estimated from the acquired seismic data.

[0071] The embodiments of the present invention described above relate to a seismic receiver that is provided with sensors for measuring the pressure p of incident seismic energy. As previously noted, however, it is also possible to provide a seismic receiver with sensors for measuring particle velocity {overscore (ν)}. An alternative embodiment of the invention is suitable for use with a seismic receiver, for example a streamer, that has both pressure sensors and particle velocity sensors. Such a streamer can directly measure both the fluid pressure p and the particle velocity {overscore (ν)}. The quantities p and {overscore (ν)} are not independent from one another, however, but are related by equation (1) or (2), so that there is an extra degree of freedom in the subsequent processing of the seismic data acquired by the streamer.

[0072] In one embodiment of the invention, it is assumed that the depth of the streamer is known. This embodiment allows the output from the particle velocity sensors to be compared with the output from the fluid pressure sensors. This comparison can be made, for example, by calculating the particle velocity from the measured values for pressure and from the depth of the streamer, and comparing these calculated values of the particle velocity with the directly measured values for the particle velocity. These values should, of course, be the same if the particle velocity and the fluid pressure have been measured accurately. Alternatively, the pressure can be calculated from the measured values for the particle velocity and from the depth of the streamer, and these calculated values of the pressure can be compared with the directly measured values for the pressure. This process can be considered as calibrating the pressure sensors relative to the particle velocity sensors, and allows the accuracy of the sensors to be monitored.

[0073] This process for monitoring the accuracy of the sensors can be carried out while acquired seismic data is being processed by a de-ghosting method of the present invention. Alternatively, it can be carried out on its own, for example to check the sensors at the start of a seismic survey. This allows the extra degree of freedom provided by measuring both p and {overscore (ν)} to be used to calibrate the sensors.

[0074] In this embodiment, the depth of the streamer can be estimated, for example from knowledge of the arrangement for suspending the streamer from floats. For improved accuracy, however, it is preferable that the depth of the receiver is measured directly.

[0075] In another embodiment of the invention, the extra degree of freedom resulting from the measurement of both p and {overscore (ν)} is used to calculate the distance between each pressure sensor and the sea-surface. In the case of a flat sea, this gives the depth and shape of the streamer.

[0076] A further embodiment of the invention again relates to monitoring the accuracy of a streamer that is provided with both fluid pressure and particle velocity sensors. In the embodiment the accuracy of the particle velocity sensors is monitored by using equation (4), which states that the curl of the particle velocity should be identically zero. This in principle allows the particle velocity sensors to be monitored since, if {overscore (ν)} has been measured accurately, then calculating {overscore (∇)}x{overscore (ν)} should give zero so that {overscore (∇)}x{overscore (ν)} acts as a measure of the quality of the measurement of {overscore (ν)}. In practice, however, it has hitherto been difficult to compute {overscore (∇)}x{overscore (ν)} with sufficient accuracy for this to be a reliable method of monitoring the particle velocity sensors. Measurements of {overscore (ν)} are available at many points along the streamer but are only available at one depth, and this makes it difficult to compute the derivatives with respect to depth of components of {overscore (ν)}.

[0077] According to this embodiment of the invention, a Lax-Wendroff correction is used to improve the accuracy of the computation of the derivatives with respect to depth of the components of {overscore (ν)}, and thereby improve the accuracy with which {overscore (∇)}x{overscore (ν)} can be computed. The wave equation (3) can alternatively be written in terms of {overscore (ν)} rather than p, and using a Lax-Wendroff correction enables the derivative with respect to depth of the components of {overscore (ν)} to be expressed in terms of derivatives with respect to time and derivatives with respect to distance along the streamer. This enables the derivative with respect to depth of the components of {overscore (ν)} to be computed with greater accuracy; this is analogous to the increase in accuracy of computing ∂_zzp in previous embodiments.

[0078] A further embodiment of the invention is for use with a configuration in which two streamers are towed one above the other in the vicinity of the sea-surface. In this configuration, the pressure is measured at the level of each streamer, and it is also known that the pressure at the sea-surface is zero. This receiver configuration therefore provides sufficient information to enable the second-derivatives of pressure in both the horizontal direction and the vertical direction to be calculated from the measured pressure data. It is thus possible to determine ∂_zp directly from equation (5), since all the terms on the right-hand side of equation (5) can be computed directly from the measured pressure data. In this embodiment, the wave equation (3) which relates spatial derivatives of p to time derivatives of p, is not required in the computation of ∂_zp, and can be used to provide a constraint on the streamer depth or on the streamer separation in a de-ghosting step, since limits on the streamer depth or streamer separation can be determined.

[0079] Alternatively in this embodiment of the invention, a Lax-Wendroff relation can be used to increase the accuracy of the estimate of the vertical derivative between the streamers. This will allow the streamers to be towed further apart from each other while properly handling higher frequencies in the data. For this situation, the depth derivative is centred at a depth midway between the depths of the two streamers so that the Taylor expansion of the approximation takes the form:

[0080] $\begin{array}{cc}{\partial }_{z}p\left(x_0,z_0\right)=\frac{p\left(x_0,z_0+h/2\right)-p\left(x_0,z_0-h/2\right)}{h}-\frac{h^2}{24}{\partial }_{\mathrm{zzz}}p\left(x_0,z_0\right)+O\left(h^4\right)& \left(10\right)\end{array}$

[0081] In Equation (10), h is the vertical separation between the two streamers.

[0082] Again, by using equation (3) to write ∂_zzzp in terms of derivatives with respect to time or the x-direction, it is possible to obtain an implicit system of equations to solve for ∂_zp.

[0083] The embodiments of the present invention described above have related to the de-ghosting the signal recorded at the receiver. The invention can, however, also be used for source-side de-ghosting. In principle this could be done by directly measuring the fluid pressure at the source, but in practice it is possible to make use of the principle of reciprocity to avoid the need to do this.

[0084] The principle of reciprocity is a fundamental principle of wave propagation, and states that a signal is unaffected by interchanging the location and character of the sources and receivers. For example, if a surveying arrangement with a seismic source at point A and a receiver array at point B gives a certain signal at the receiver array, then using a single receiver at point A and a source array at point B would lead to the same signal, provided that the source array corresponds to the receiver array. (By “corresponds”, it is meant that the source array contains the same number of sources as the receiver array has receivers, and that the sources in the source array are arranged in the same positions relative to one another as the receivers in the receiver array).

[0085] In consequence of the principle of reciprocity, the signal emitted by an array of seismic sources can be de-ghosted by making measurements using a receiver array, whether the seismic sources are actuated simultaneously or consecutively. The recorded signal is analogous to that recorded in a reciprocal arrangement where the source and receiver locations are interchanged. The method outlined hereinabove can therefore also be applied to the source array that it is desired to de-ghost. The signal produced at the receiver array by an array of a plurality of seismic sources is measured, and a de-ghosting filter is derived from this measured signal. By the principle of reciprocity, this filter can be used to de-ghost the signal emitted by the source array.

[0086] The present invention can provide good results in source-side de-ghosting, since the depth of the source(s) in a typical marine seismic survey is generally well-known, for example because the source(s) is/are suspended from a boat by a known length of cable. The present invention can particularly be applied to a full three-dimensional VSP (vertical seismic profiling) acquisition pattern, where a large number of seismic sources are generally employed in a relatively small area. The full three-dimensional version of the present invention can then be applied efficiently.

[0087] The methodology of the present invention can also be used to design seismic source arrays that have two or more marine vibrators, air guns, or water guns.