WO2016155771A1

WO2016155771A1 - Deghosting method

Info

Publication number: WO2016155771A1
Application number: PCT/EP2015/056884
Authority: WO
Inventors: Lasse Amundsen; Johan Olof Anders Robertsson
Original assignee: Statoil Petroleum ASA
Current assignee: Equinor Energy AS
Priority date: 2015-03-30
Filing date: 2015-03-30
Publication date: 2016-10-06
Anticipated expiration: 2017-09-30

Abstract

A method of deghosting seismic data using a model, the seismic data having been captured at a location below the sea surface and the model comprising a wave-propagation space having the wave propagation properties of sea-water, the method comprising: a. injecting into the model a seismic wavefield based on the seismic data at an injection location in the model corresponding to the location at which the seismic data was captured, and allowing the injected seismic wavefield to propagate through the wave-propagation space; b. recording the propagating seismic wavefield at a recording location in the model, the recording location being located at a distance from the injection location of twice the depth of the location below the sea surface at which the seismic data was captured; c. scaling or convolving the recorded seismic wavefield with a reflection coefficient, R; and d. adding the scaled or convolved seismic wavefield to a subsequent seismic wavefield, which is based on a corresponding subsequent portion of the seismic data, the subsequent seismic wavefield being the seismic wavefield that is to be injected at the injection location at the same time as the propagating seismic wavefield passes the recording location.

Description

DEGHOSTING METHOD

The present invention relates to methods of deghosting seismic data and to methods of reducing the noise in deghosted seismic data.

Seismic ghosts are a long-standing issue in the marine seismic exploration industry.

A source ghost is an event starting its propagation upward from the seismic source, and a receiver ghost ends its propagation moving downward at the receiver. They both have a reflection at the sea surface, which leads to a reduction of the useful frequency bandwidth and therefore damages seismic resolution. Ghosts have traditionally proved difficult to eliminate through data processing, leading to attempts to eliminate them in data acquisition through methods, such as hydrophone-geophone streamers, dual-streamer or slanted- streamer towing (e.g., Ozdemir et al., 2008; Day et al., 2013; Soubaras and Lafet, 2013). Although significant progress has been made over the last several years, the problem is still not fully solved.

The ghost removal process is known as deghosting. Deghosting can be applied both at the receiver side and on the source side. The inventors have described a new method for deghosting and other applications of wavefield separation and prediction that relies on multicomponent recordings (Amundsen and Robertsson, 2014). However, multicomponent recordings require special consideration when they are recorded. It is therefore desirable to have a method of deghosting that does not rely upon multicomponent recordings, such as a method that uses pressure data.

Receiver-side wave equation based deghosting techniques for conventional pressure recorded data are described in, e.g., Amundsen (1993), Fokkema and van den Berg (1993), Robertsson and Kragh (2002), Amundsen et al. (2005), and Weglein et al. (2002).

Amundsen et al. (2013) show that receiver-side wave equation based deghosting can be derived and implemented in the frequency domain as spatial deterministic deconvolution of marine pressure recordings. The deghosting/deconvolution operator has however poles at fundamental frequencies equal to an integer multiple of a function of the receiver depth. Beasley et al. (2013a, 2013b) and Beasley and Coates (2014) present a time-space domain approach where the wave equation is used to simulate propagation of the up- and down- going wavefields between the receivers and the water surface to accomplish wavefield separation. They suggest that the up-going, deghosted wavefield can be computed iteratively by first evaluating a Kirchhoff integral over earlier (and previously computed) values of the up-going wavefield which is then subtracted from the recorded pressure wavefield which yields the next value of the up-going wavefield. They also mention the use of two-way finite-difference propagators such as the ones used in this application and by Amundsen and Robertsson (2014) and Robertsson and Amundsen (2014). Receiver-side deghosting techniques have specifically been developed for variable depth streamers. Soubaras and Lafet (2013) propose to use a dual-input, migration and mirror migration, and a joint deconvolution algorithm to produce a 3D deghosted output. Also integral equation solutions to deghosting for seismic pressure wavefields that have been recorded on non-horizontal streamers have been published (Riyanti et al., 2008; van Borselen et al., 2013; Amundsen and Reitan, 2014). Van Borselen et al. (2013) state that subject to the depth profile of the steamer (and the resulting notch diversity), the integral equation can be ill posed, because the integral kernel has zeroes at the notch frequencies. After receiver-side deghosting, the source ghost is still present in the seismic data. Most published approaches to source-side deghosting are based on vertical source array acquisition, in which sources are towed in an over/under fashion, such as the procedures suggested by Moldoveanu (2000, 2001 ), Moldoveanu et al. (2007), Robertsson et al. (201 1 ), and those published in Vaage (2005). Recently, a so-called ghost-free solution has been introduced which requires a time and depth distributed source, using subsources deployed at different depths and fired with specific time delays (Parkes and Hegna, 201 1 ). Further, it has been proposed to apply Green's theorem also for source-side deghosting (see Mayhan and Weglein, 2013).

Although new technologies are being introduced, there is still a pressing need to develop and improve source- and receiver-side deghosting solutions for conventional towed pressure streamer data, as most published deghosting solutions for such data suffer from stability-related issues. The concepts of claims 1 , 2, 6, 7, 14, 17, 18, 19, 22 and 27 have been published by the inventors in Robertsson and Amundsen (2014). The method is related to the method presented by Beasley et al. (2013a, 2013b) and Beasley and Coates (2014) in the respect that it assumes causality to predict the down-going wavefield from the up-going wavefield using a wave equation propagator. However, as is discussed below the present invention provides simplified methods that are more stable. Further, methods are disclosed that obtain a stable solution with minimal noise, avoiding some fundamental assumptions of the method disclosed in Robertsson and Amundsen (2014) (e.g., that the first arrival must be up-going). In a first aspect, the invention provides a method of deghosting seismic data using a model, the seismic data having been captured at a location below the sea surface, and the model comprising a wave-propagation space having the wave propagation properties of sea- water, the method comprising: a. injecting into the model a seismic wavefield based on the seismic data at an injection location in the model corresponding to the location at which the seismic data was captured, and allowing the injected seismic wavefield to propagate through the wave-propagation space; b. recording the propagating seismic wavefield at a recording location in the model, the recording location being located at a distance from the injection location of twice the depth of the location below the sea surface at which the seismic data was captured; c. scaling or convolving the recorded seismic wavefield with a reflection coefficient, R; d. adding the scaled or convolved seismic wavefield to a subsequent seismic wavefield, which is based on a corresponding subsequent portion of the seismic data, the subsequent seismic wavefield being the seismic wavefield that is to be injected at the injection location at the same time as the propagating seismic wavefield passes the recording location.

The method is advantageous as there is no need for multi-component seismic data to be used. Further, the method is stable for all frequencies of the seismic data, i.e. there are no frequencies that cannot be handled due to, for example, notching. Further still, the method can be used both for receiver-side deghosting and source-side deghosting (receiver- side deghosting takes place in the common shot domain, whereas source-side deghosting operates in the common receiver domain). Further still, the method can handle seismic data which has been captured at varying depths (i.e. not constant depth). In the method it is merely required that the depth of the captured seismic data is known, however this depth can vary over time. Further still, the model of the present invention is simple and robust. For instance, there is no need for a free surface to be present in the model.

The method assumes causality to predict the down-going wavefield from the up- going wavefield, and may use a wave-equation propagator in the wave-propagation space, preferably a two-way wave equation finite-difference propagator.

The seismic data may relate to primary wavefields or multiples.

The seismic data may be marine seismic hydrophone data. The method works for conventional hydrophone streamer data, even when the streamer is towed at arbitrary and variable depths. Apart from the captured seismic data, only the depth profile of the streamer recording the data beneath a sea surface with a known reflection coefficient as well as the propagation velocity in water above the streamer need be known. The method may comprise the step of capturing the seismic data.

The seismic data may consist of pressure data. The seismic data may consist of velocity data. The velocity data may be single component velocity data. The velocity data may be particle velocity data. The single component may the vertical velocity component (v_z).

As mentioned above, the method does not require multi-component seismic data, thus mere single-component seismic data can be used. This allows the method to be used for more seismic data sets and eases the capturing of the seismic data.

The wave-propagation space may be rectangular. The rectangle may be orientated such that it has a horizontal lower and upper edges and vertical side edges. The horizontal and vertical edges correspond to the horizontal and vertical directions of the sea in which the seismic data was captured.

The model may be a 2D model (e.g. rectangular). The model may be a 3D model (e.g. cuboid). In the preceding and foregoing discussion, the "rectangle" may also refer to cuboid-shaped models, and an "edge" of the model may refer to a surface of the 3D model.

The recording location may be vertically above or vertically below the injection location in the model. This is advantageous because this location of the recording location best correlates the geometry of the model to the real subsea geometry.

The wave-propagation space may comprise edges, at least some of which have absorbing boundary conditions. All of the edges may have absorbing boundary conditions.

Due to the use of a wave-propagation space in the model, after injecting single- component seismic data into the wave-propagation space the injected seismic wavefield will be a multi-component wavefield (e.g. it will consist of both pressure and velocity seismic data).

Since single-component, preferably directionless, seismic data may be used (e.g. pressure data, or vertical component velocity data), it should be appreciated that when such seismic data is injected into a wave-propagating model, the injected wavefield normally propagates in all directions. However, in the present case, the seismic data actually captured is mostly up-going seismic data (with down-going ghosts). Thus, if this is injected into the wave-propagating model, an inaccurate wavefield solution will be found. It is therefore important to compensate for the directionless nature of the captured seismic data.

The inventors have devised two solutions for doing so.

The first solution is that the injection location may be located at an edge portion of the wave-propagation space. In this case, the edge portion of the model may comprise rigid or free boundary conditions. The boundary conditions of the edge portion may be Neumann or Dirichlet boundary conditions. When single-component velocity seismic data (e.g. v_z data) is injected the boundary condition is preferably Dirichlet, when pressure data is used the boundary condition is preferably Neumann. The rigid or free boundary provides the correct propagation of the injected wavefield. The rigid or free boundary condition forces the seismic wavefield injected at the location to propagate away from the boundary and hence toward the recording location. At least some of the remaining edges of the wave- propagation space may have absorbing boundary conditions, for example convolutional perfectly matched layers (PML). This reduces reflections within the wave-propagation space.

The edge portion may be the lower edge of the rectangle. The remaining edges of the rectangle may form the remaining edges of the wave-propagation space. The edge portion may extend in a substantially rectilinear direction.

In this case, because the wavefield is injected at an edge portion as described above, the wavefield propagates only away from the edge portion, and hence in one direction away from the injection location. For example, when the edge portion is the edge of the rectangle, the wavefield propagates generally perpendicularly away from the edge portion. For example, when the edge portion is the lower edge of the rectangle, the wavefield generally propagates upwards. Thus, all of the injected seismic wavefield may propagate in one direction, for example the direction in which it is known for it generally to have been propagating when it was captured, i.e. upwards.

Since all the injected seismic wavefield propagates in the "correct" direction, in this embodiment it is possible to use a reflection coefficient at the recording location of 0 < R < 1 , preferably 0.95. See below for more discussion on the reflection coefficient.

The second solution is that the injection location may be at any location within the wave-propagation space (i.e. within the space itself, or at or near the edges). In this case, at least two edges of the wave-propagation space (e.g. the upper and lower edges, or the entire perimeter) may have absorbing boundaries. This prevents unwanted reflections.

However, in this embodiment, it will be appreciated that only half of the injected seismic wavefield propagates towards the recording location from the injection location. For example, when the recording location is vertically above the injection location, half of the injected wavefield propagates upwards and the other half propagates downwards. Thus, the recorded wavefield has an amplitude of approximately one half of the "correct" value. In order to compensate for this, a reflection coefficient of 0 < R < 2, preferably 2, preferably 1 .95, preferably 1 .9, preferably 1.8 may be used. By using a larger reflection coefficient, the reduction in amplitude can be compensated for.

The absorbing boundary conditions may be Perfectly Matched Layers (PML), for example convolutional PMLs.

The injected seismic data may drive the model.

The steps of a to d may be occurring continuously, i.e. seismic data may be injected as the wavefield propagates through the model. Steps a to d may therefore be considered to be repeated for every time sample of seismic data injected into the model. Here, "continuous" may include sequential discrete data separated in time (which is typically how the seismic data is captured and processed). A subsequent piece of seismic data may be injected into the model prior to the wavefield from the preceding piece of seismic data being recorded.

The complete deghosted seismic data produced by performing the method for all of a seismic data set captured over a period of time may be referred to as a solution.

The wave-propagation space may comprise two side edge portions opposite one another. The method may comprise performing steps a to d with at least one side edge portion of the model comprising a first boundary condition to provide first deghosted seismic data; performing steps a to d with the at least one side edge portion of the model comprising a second complimentary boundary condition to provide second deghosted seismic data; and summing the first and second deghosted seismic data.

Thus, the steps a to d may be effectively run through twice, with one of the side edge portions of the model having a first boundary condition on the first run and a second boundary condition on the second run, the second boundary condition being complimentary to the first boundary condition.

Likewise, the same boundary conditions may be applied to the other side edge portion during the first and second runs.

Complimentary boundary conditions may be Neumann and Dirichlet boundary conditions. The same seismic data may be used to calculate the first and second deghosted seismic data. The first and second boundary conditions may be complimentary.

The first and second deghosted seismic data calculated in such a way can be used to remove all of the reflections from the side edge portions of the wave-propagation space. This is at least true for wavefield energy that interacts once with the side edge portions. The first order reflection vanishes once the first and second deghosted seismic data is summed because of the complimentary boundary conditions used. This method of removing reflections is particularly useful when the seismic data has been captured close to the source of seismic energy.

The side edge portions may be vertically orientated.

The side edge portions may join the edge portion. The edge portion may be horizontally orientated. The wave-propagation space may comprise a fourth edge portion, opposite the edge portion and joining the two side edge portions. The fourth edge portion may have absorbing boundary conditions, e.g. PML.

The side edge portions may be joined by a third and a fourth edge portion, the third and fourth edge portions being opposite one another. The third and fourth edge portions may be parallel to one another. The side edge portions and the third and fourth edge portions may form the edges of the rectangle space (or the faces of the cuboid space). The third and fourth edges may have absorbing boundary conditions, e.g. PML.

In the case of a cuboid, the wave-propagation space may comprise fifth and sixth surfaces closing the cuboid (i.e. joining to the first, second, third and fourth surfaces, at opposite ends thereof). The fifth and sixth surfaces may have absorbing boundary conditions (e.g. PML), or complimentary boundary conditions (e.g. Neumann and Dirichlet).

The model may be a finite-difference model, a finite-element model, a spectral element model, or a pseudo-spectral model. The wave-propagation space may be a grid, e.g. a finite-difference wave-propagation grid. The method may be a finite-difference method, a finite-element method, a spectral element method, or a pseudo-spectral method.

The model may be in the space-time domain. The axes of the wave-propagation space may be spatial, i.e. depicting distances. The seismic data may be injected at a location corresponding to the location at which they were captured. The seismic data may be injected over a time and at a rate corresponding to the time period over and the rate at which the seismic data were captured.

Working in the space-time domain means that the method is unaffected by any number of notches in the frequency spectrum of the data and therefore will deghost "through notches". Notches are inherent in any seismic data frequency spectrum where ghosts are present.

The location at which the seismic wavefield is injected into the model corresponds to the location beneath the sea at which the seismic wavefield data were captured.

The seismic wavefield may be injected at a point source or a plurality of point sources. The point source(s) may be monopole(s). The plurality of point sources may correspond to the location at which the seismic data was captured using a streamer.

The respective point sources may be at the same depth as one another. The respective point sources may be spaced horizontally from one another.

The depth may be constant over time. The depth of all of the point sources may vary together over time. The depth of the point sources may vary with respect to each other and the sea surface over time.

The respective point sources may be at varying depths with respect to one another.

The respective depths may be constant over time. However, the respective depths may vary over time.

The configuration of the point sources merely depends on the configuration of the equipment used to capture the data, e.g. the streamer.

The propagating seismic wave may be recorded at a plurality of locations. The recording locations may comprise a plurality of point recording data corresponding to the plurality of point sources. The configuration of the recording datum merely depends on the configuration of the source(s) and/or the varying depth of the data-capturing locations, which may be due to a varying depth of the streamer and/or having a rough sea surface.

At step d, the recorded wavefield at a recording datum may be added to the subsequent seismic wavefield being injected at the corresponding source location.

The recording datum/data may be spaced at a distance of twice the depth of the captured seismic data substantially vertically away from the corresponding source. The recording datum/data may be spaced at a distance of twice the depth of the captured seismic data substantially perpendicularly to the edge portion.

The model may not include any free surface. A free surface is not required. This simplifies the method in comparison to the methods of the prior art.

The direct wave may have been removed from the seismic data.

The seismic data may comprise captured seismic data, the captured seismic data missing near offset seismic data, and reconstructed near offset seismic data, the method may comprise reconstructing the reconstructed near offset seismic data. It is common to have missing near offset seismic data, however it is best to have complete data sets. Thus, it is useful to reconstruct the missing near offset data. Such reconstruction can be performed by normal moveout (NMO).

The seismic data may comprise captured seismic data and reconstructed negative offset seismic data, the method comprising reconstructing the reconstructed negative offset seismic data. The use of negative offset seismic data is particularly important in regions of dipping geology. Negative seismic data can be reconstructed by mirroring positive offset seismic data.

A value of 0 < R < 1 may be used in the method. Using R of less than 1 is advantageous since it helps to stabilise the deghosting method, and is closer to the actual reflection coefficient of the sea-air interface. Preferably an R of 0.95 or 0.9 is used. As mentioned above, this value may be used when the wavefield is injected at the edge portion of the model.

A value of 0 < R < 2 may be used in the method. Using R of less than 2 is advantageous since it helps to stabilise the deghosting method, and is closer to the actual reflection coefficient of the sea-air interface. Preferably an R of 1.95 or 1.9 or 1 .8 may be used. As mentioned above, this value may be used when the injection location is not at the edge of the model and/or all of the edges of wave model are absorbing boundaries.

The reflection coefficient of the sea surface is approximately -1 . Thus, in order to remove the ghost reflections from the subsequent seismic wavefield, the reflection coefficient used in the present method should be of opposite polarity and should best match the amplitude of the scaled/convolved seismic wavefield to the amplitude of the ghost wavefield. Thus, an R of approximately 1 (or approximately 2, in the case where the amplitude of the recorded wavefield is one half of the "correct value", as discussed above) should be used. Using an R of 0.95 or 0.9 (or 1 .95, 1.9, or 1.8) can help to stabilise the solution.

The reflection coefficient, R, may be a function that may vary with time. Such a reflection coefficient may allow for changing surface reflection coefficients that vary in time.

Conformal mapping may be used to model the effects of a rough sea surface and/or a depth varying recording streamer. Specific examples of using conformal mapping for such a purpose are given in Fornberg (1988).

The method may comprise performing steps a to d using forward time extrapolation (FTE) to produce forward time extrapolated seismic data. In this case, it is assumed that the initial captured/injected seismic data corresponds to purely an up-going seismic wave. An example of this method is given in the specific description below.

The method may comprise performing steps a to d using reverse time extrapolation (RTE) to produce reverse time extrapolated seismic data. This is advantageous since the assumption that the initial captured/injected seismic data corresponds to purely an up-going seismic wave (i.e. and does not include any down-going component) is avoided.

This assumption is avoided by instead assuming that the last arrival from the total captured wavefield is a down-going ghost reflection. The fact that the last received wavefield is down-going can be used to predict the down-going wavefield at all times, via reverse-time extrapolation. This extrapolation is similar to the example of the FTE method discussed in the specific description, except that the first data used is assumed to be from a purely down-going reflection and the extrapolation is done in reverse in time.

Thus, using the RTE method, the down-going wavefield may be calculated which can then be, at the end, subtracted from the total recorded wavefield to estimate the up-going (deghosted) wavefield.

The assumption that the latest captured seismic wavefield arrival is down-going is much more robust than the assumption that the first captured seismic wavefield arrival is up- going as there is very little energy present at the end of the captured seismic record. The last energy that is recorded it will always be down-going as it will be the ghost of the latest arrival from the sub-surface.

Further, since it is performed in reverse time, it is the weakest arrivals (and hence weakest energy) of the captured data that are injected into the model first. This helps to increase the signal-to-noise ratio in RTE implementation, in comparison to the FTE implementation, since any noise generated from weaker (deeper) arrivals will interfere with stronger shallower arrivals. The inventors have found that the FTE and the RTE solutions have identical deghosted signal but different noise. Thus, a combination of the FTE and RTE can be used to calculate the noise, as is discussed further below.

In a second aspect, the invention provides a method of reconstructing a multicomponent wavefield from seismic data using a model, the seismic data having been captured at a location below the sea surface, and the model comprising a wave-propagation space having the wave propagation properties of sea-water and a reflecting surface at a location corresponding to the location of the sea surface, the method comprising: a. injecting into the model a seismic wavefield based on the seismic data at an injection location in the model corresponding to the location at which the seismic data was captured; b. introducing a boundary condition at the injection location; and c. allowing the wavefield to propagate through the wave-propagation space and reflect off the reflecting surface, thus

reconstructing the total wave-field.

Such a reconstructed multicomponent wavefield can be used, for instance, to calibrate multicomponent seismic data or can be used in performing deghosting.

The method is advantageous as there is no need for multi-component seismic data to be used in order to produce a reconstruction of a multicomponent seismic wavefield.

Further, the method is stable for all frequencies of the seismic data, i.e. there are no frequencies that cannot be handled due to, for example, notching. Further still, the method can be used both for receiver-side deghosting and source-side deghosting (receiver-side deghosting takes place in the common shot domain, whereas source-side deghosting operates in the common receiver domain). Further still, the method can handle seismic data which has been captured at varying depths (i.e. not constant depth). In the method it is merely required that the depth of the captured seismic data is known, however this depth can vary over time. Further still, the model of the present invention is simple and robust. For instance, there is no need for a free surface to be present in the model.

The seismic data may relate to primary wavefields or multiples.

The seismic data may consist of velocity data. The velocity data may be single component velocity data. The velocity data may be particle velocity data. The single component may the vertical velocity component (v_z). In this case, the boundary condition at the injection location may preferably be Dirichlet boundary conditions. The seismic data may consist of pressure data. In this case, the boundary condition at the injection location may preferably be Neumann boundary conditions.

The reconstructed seismic wavefield may be used to calibrate seismic data as follows. It is usually the case that the captured pressure data are of much higher quality than the captured single component velocity (e.g. v_z) data and there is therefore a need to calibrate the single component velocity data, particularly at low frequencies. The method of this aspect can be used to generate total seismic wavefield data from single component velocity data. The total seismic wavefield data includes seismic pressure data. This reconstructed seismic pressure data can then be compared to captured pressure data corresponding to the captured velocity data. Any differences may indicate a need to (and can be used to) calibrate the captured velocity data.

The reflecting surface may be orientated horizontally and may be vertically above (or below) the injection location. In the case where there are more than one injection locations, the reflecting surface may be orientated generally parallel to the direction between adjacent injection locations. The purpose of the reflecting surface is to model the sea surface in the model.

The wave-propagation space may comprise edges that have absorbing boundary conditions.

Since single-component, directionless seismic data may be used (e.g. pressure data), it should be appreciated that when such seismic data is injected into a wave- propagating model, it propagates in all directions. However, in the present case, the seismic data actually captured is mostly up-going seismic data (with down-going ghosts). Thus, if this is injected into the wave-propagating model, an inaccurate wavefield solution will be found. It is therefore important to compensate for the directionless nature of the captured seismic data. This is done by introducing the boundary condition at the injection location. (Alternatively, it could be said that the model comprises a surface at the injection location, the surface having the boundary condition.)

Due to the use of a suitable wave-propagation space in the model, after injecting single-component seismic data into the wave-propagation space the injected seismic wavefield will be a multi-component wavefield (e.g. it will consist of both pressure and velocity seismic data).

In order to compensate for the directionless nature of the captured seismic data, the present aspect includes the step of introducing a boundary condition at the injection location. The boundary condition may be a boundary condition such that the injected seismic wavefield generally only propagates toward the reflecting surface (e.g. upwards in the model, when the reflecting surface is above the injection location). Thus, the injected seismic wavefield does not propagate away from the reflecting surface.

The boundary condition may be an artificial free-surface boundary condition.

The boundary condition may have a reflection coefficient of -1 .

The boundary condition may be a Dirichlet boundary condition.

The boundary condition may be a Neumann boundary condition.

The model may comprise one or more sources where the seismic wavefield is injected. The sources may be superimposed on top of the boundary condition. The boundary condition may be superimposed on top of the sources.

The boundary condition may form a second reflecting surface in the model at the injection location(s). The second reflecting surface may be parallel to the reflecting surface.

The method may comprise: d. allowing the propagating seismic wave that has reflected off the reflecting surface in step c to propagate back to the injection location where it is again reflected back towards the reflecting surface.

Thus, when single component seismic data (e.g. pressure data or single-component velocity data) is injected into the model, the resulting injected wavefield is effectively given a direction of propagation by the boundary condition at the injection location. The direction may be generally towards the reflecting surface. Thus, all of the injected wavefield may propagate towards the reflecting surface.

Thus, this method also relies on the assumption that, when using FTE, the first captured seismic data corresponds to a purely up-going wavefield and, when using RTE, the last captured seismic data corresponds to a purely down-going wavefield.

Due to the presence of the boundary condition at the injection location, a full wavefield can be reconstructed starting from only single-component seismic data (e.g.

pressure data or single-component velocity data). The boundary condition may force the wavefield to propagate toward the reflecting surface (e.g. upwards). This is acceptable for the initial seismic data, since the initial seismic data is upwards only (in FTE) or downwards only (in RTE). However, subsequent seismic data can have both upward and downward components. Forcing all of the subsequent injected wavefields towards the reflecting surface is therefore not itself accurate. However, due to the reflecting surface, the reflected wavefield in the propagating space of the model interferes with subsequently injected wavefields, thus cancelling out the portion of the injected wavefield that in the real situation (i.e. in the real sub-sea situation) was actually propagating downwards. Thus, between the injection location and the reflecting surface in the model, the only wavefields are those corresponding to the upward propagating wavefields in the real situation, and their downward propagating reflections. Thus the full wavefield present in reality between the capturing location (e.g. the streamer location, which corresponds to the injection location in the model) and the sea surface (which corresponds to the reflecting surface in the model) is reconstructed by repeating steps a to c for all the seismic data. From the full wavefield, deghosting can be performed using any known technique of deghosting using a full reconstructed wavefield.

Put another way, the reflected wavefield will arrive at the injection location at the time that the corresponding down-going seismic data is injected upwards. However, the reflected wavefield interferes with this injected wavefield so that no wavefields that are down-going in the real sub-sea situation propagate upwards in the model. (Whilst "up" and "down" are used here, it is clear that in relation to the model these can be any two opposite directions.) Thus, this method can construct the full multi-component wavefield between the injection location and the reflecting surface in the model by simply injecting single- component seismic data.

The injection location may be located at an edge portion of the wave-propagation space. Thus, the boundary condition may be introduced at the edge portion. At least some of the remaining edges of the wave-propagation space may have absorbing boundary conditions, for example convolutional perfectly matched layers (PML). This reduces reflections within the wave-propagation space.

Because the wavefield is injected at an edge portion as described above, and because an appropriate boundary condition is introduced at the injection location, the wavefield propagates only away from the edge portion, and hence in one direction away from the injection location. For example, when the edge portion is the edge of the rectangle, the wavefield propagates generally perpendicularly away from the edge portion. For example, when the edge portion is the lower edge of the rectangle, the wavefield generally propagates upwards. Thus, all of the injected seismic wavefield may propagate in one direction, for example the direction in which it is known for it generally to have been propagating when it was captured, i.e. upwards.

The reflecting surface may be a free surface in the model.

The reflection coefficient of the reflecting surface may be -1 < R < 0, preferably -1 , preferably -0.95.

The model may be a 2D model (e.g. rectangular). The model may be a 3D model (e.g. cuboid). In the preceding and foregoing discussion, the "rectangle" may also refer to cuboid-shaped models, and the "edge" may refer to a surface.

The injected seismic data may drive the model.

The steps of a to c may be occurring continuously, i.e. seismic data may be injected as the wavefield propagates through the model. Steps a to c may therefore be considered to be repeated for every time sample of seismic data injected into the model. Here,

"continuous" may include sequential discrete data separated in time (which is typically how the seismic data is captured and processed). A subsequent piece of seismic data may be injected into the model prior to the wavefield from the preceding piece of seismic data being reflected.

The wave-propagation space may comprise two side edge portions opposite one another. The method may comprise performing steps a to c with at least one side edge portion of the model comprising a first boundary condition to provide a first reconstructed wavefield; performing steps a to c with the at least one side edge portion of the model comprising a second complimentary boundary condition to provide a second reconstructed wavefield; and summing the first and second reconstructed wavefields.

Complimentary boundary conditions may be Neumann and Dirichlet boundary conditions. The same seismic data may be used to calculate the first and second reconstructed wavefields. The first and second boundary conditions may be complimentary.

Likewise, the same boundary conditions applied to the side edge portion above may also be applied to the other side edge portion during the first and second runs through steps a to c.

The first and second reconstructed wavefields calculated in such a way can be used to remove all of the reflections from the side edge portions of the wave-propagation space. This is at least true for wavefield energy that interacts once with the side edge portions. The first order reflection vanishes once the first and second reconstructed wavefields are summed because of the complimentary boundary conditions used. This method of removing reflections is particularly useful when the seismic data has been captured close to the source of seismic energy. As will be appreciated, the reconstructed wavefields with removed reflections can be then be used to produce more accurate deghosted seismic data.

The side edge portions may be vertically orientated.

The model may be in the space-time domain. The axes of the wave-propagation space may be spatial, i.e. depicting distances. The seismic data may be injected at a location corresponding to the location at which they were captured in the real subsea situation. The seismic data may be injected over a time and at a rate corresponding to the time period over and the rate at which the seismic data were captured.

Working in the space-time domain means that the method is unaffected by any number of notches in the frequency spectrum of the data and therefore will reconstruct "through notches". Notches are inherent in any seismic data frequency spectrum where ghosts are present.

The configuration of the point sources merely depends on the configuration of the equipment used to capture the data, e.g. the streamer. The distance between the injection location(s) and the reflecting surface merely depends on the depth of the equipment used to capture the data, e.g. the streamer.

The direct wave may have been removed from the seismic data.

As mentioned above, a value of -1 < R < 0 for the reflection coefficient of the reflecting surface may be used in the method. Using R of greater than -1 is advantageous since it helps to stabilise the reconstructing method, and is closer to the actual reflection coefficient of the sea-air interface. Preferably an R of -0.95 or -0.9 is used.

The reflection coefficient of the real sea surface is approximately -1. Thus, when a the reflecting surface in the model has -1 < R < 0 and the boundary condition imposed at the injection location has a reflection coefficient of -1 , the wavefield that is reflected off the boundary condition at the injection location and hence is upwardly propagating in the model (i.e. the wave that is injected upwards, then reflected downwards off the reflecting surface and reflected upwards again due to the boundary condition at the injection location) has opposite polarity compared to the injected data and will therefore destructively interfere to remove the erroneously upwards injected downgoing wavefield.

The reflection coefficient, R, may be a function that may vary with time. Such a reflection coefficient may allow for changing sea surface reflection coefficients that vary in time.

The method may comprise performing steps a to c using forward time extrapolation (FTE) to produce forward time extrapolated seismic data. In this case, it is assumed that the initial captured/injected seismic data corresponds to purely an up-going seismic wave.

The method may comprise performing steps a to c using reverse time extrapolation (RTE) to produce reverse time extrapolated seismic data. This is advantageous since the assumption that the initial captured/injected seismic data corresponds to purely an up-going seismic wave (i.e. and does not include any down-going component) is avoided.

This assumption is avoided by instead assuming that the last arrival from the total captured wavefield is a down-going ghost reflection. The fact that the last received wavefield is down-going can be used to predict the down-going wavefield at all times, via reverse-time extrapolation.

Further, since it is performed in reverse time, it is the weakest arrivals (and hence weakest energy) of the captured data that are injected into the model first. This helps to increase the signal-to-noise ratio in RTE implementation, in comparison to the FTE implementation, since any noise generated from weaker (deeper) arrivals will interfere with stronger shallower arrivals.

The inventors have found that the FTE and the RTE solutions have identical reconstructed wavefield signals but different noise. Thus, a combination of the FTE and RTE can be used to calculate the noise, as is discussed further below.

The method may also be a method of deghosting seismic data further comprising: d. deghosting the seismic data using the total reconstructed wave-field. Deghosting may be carried out when the seismic data is velocity seismic data, such as v_z, or pressure seismic data.

In a third aspect, the invention provides a method of deghosting single-component seismic data comprising using a reverse time extrapolator.

The single-component seismic data may be pressure data or single component velocity data. The velocity data may be particle velocity data. The single component may the vertical velocity component (v_z).

This is advantageous since the assumption that the initial captured/injected seismic data corresponds to purely an up-going seismic wave (i.e. and does not include any down- going component) is avoided.

Thus, using the RTE method, the down-going wavefield may be calculated which can then be at the end subtracted from the total recorded wavefield to estimate the up-going (deghosted) wavefield.

Further, since it is performed in reverse time, it is the weakest arrivals (and hence weakest energy) of the captured data that are injected into the model first. This helps to increase the signal-to-noise ratio in the RTE implementation, in comparison to the FTE implementation, since any noise generated from weaker (deeper) arrivals will interfere with stronger shallower arrivals.

The method may comprise calculating a down-going wavefield from the single- component seismic data using the reverse time extrapolator and subtracting the down-going wavefield from the total wavefield. This produces a deghosted wavefield.

The method may comprise using a model having at least one edge, the at least one edge having absorbing boundary conditions.

The method may comprise using a model having at least two edges. The method may comprise producing first deghosted seismic data using the single-component seismic data and a first set of boundary conditions for at least one of the two edges; producing second deghosted seismic data using the same single-component seismic data and a second set of boundary conditions for the at least one of the edges, the first and second sets of boundary conditions being complimentary to each other; and summing the first and second deghosted seismic data.

Likewise, the same boundary conditions that are applied to the side edge portion above may also be applied to the other side edge portion during the production of the first and second deghosted seismic data.

The direct wave may have been removed from the single-component seismic data.

The single-component seismic data may comprise captured seismic data, the captured seismic data missing near offset seismic data, and reconstructed near offset seismic data, the method comprising reconstructing the reconstructed near offset seismic data.

The single-component seismic data may comprise captured seismic data and reconstructed negative offset seismic data, the method may comprise reconstructing the reconstructed negative offset seismic data.

The reconstructed negative offset seismic data may be reconstructed by mirroring positive offset seismic data.

The method of this aspect can be combined with any of the appropriate features of the previous and subsequent aspects.

The inventors have found that the FTE and the RTE solutions have identical deghosted signal but different noise. Thus, a combination of the FTE and RTE can be used to calculate the noise, as is discussed further below.

In a fourth aspect, the invention provides a method comprising subtracting forward time extrapolated deghosted seismic data from reverse time extrapolated deghosted seismic data, or subtracting the reverse time extrapolated deghosted seismic data from the forward time extrapolated deghosted seismic data.

The method may be a method of reducing the noise in deghosted seismic data.

The method may be a method of identifying the noise in deghosted seismic data.

The forward time extrapolated deghosted seismic data may have been deghosted using a forward time extrapolator and the reverse time extrapolated seismic data has been deghosted using a reverse time extrapolator. The seismic data used for both extrapolations should be the same. The method may comprise producing the forward and reverse extrapolated deghosted seismic data. Any of the appropriate methods disclosed in this application may be used to produce the deghosted seismic data.

This method allows for the noise in the deghosted seismic data to be identified, and hence suppressed. The inventors have found that the deghosted seismic data (i.e. the purely up-going seismic wavefield) produced by RTE and FTE for the same input seismic data is identical, but that each produces different noise signals. Thus, by calculating an RTE and an FTE solution for the same seismic data, and subtracting one from the other, data is produced which only contains noise.

In a fifth aspect, the invention provides a method of reducing noise in deghosted seismic data, comprising: constructing a gather of deghosted seismic data by producing forward time extrapolated deghosted seismic data for a first set of locations and producing reverse time extrapolated seismic data for a second set of locations, wherein the first and second set of locations are interleaving locations; performing a transform on the gather; muting the noise signal on the transform; and performing the inverse transform on the transform.

The "locations" may be temporal or spatial. Thus, the interleaved gather may comprise interleaving traces from adjacent spatial locations, such that every second trace comes from the FTE solution and every other trace comes from the RTE solution.

Alternatively, the interleaved gather may comprise interleaving seismic data from different subsequent times, such that every second sample in a trace may come from an RTE solution and every other sample may come from an FTE solution.

The noise signal on the transform can be identified because its location(s) will be different from the location(s) of the remainder of the signal. In the claimed gather, since the FTE and RTE produce the same deghosted seismic data, except for noise, (i.e. the purely up-going wavefield) this signal will be continuous across the gather. However, since the FTE and RTE produce the noise signals that are at least partially uncorrelated with each other, adjacent noise data in the gather will be at least partially uncorrelated. Thus, when for instance the Fourier transform of the gather is taken, noise peaks will mostly appear at around the Nyquist frequency. These peaks can be deleted. The desired noise-free deghosted signal will appear around zero.

In the example above a Fourier transform was suggested. This can be applied both in time and space (FK transform). As an alternative to an FK transform a tau-p (e.g., Radon) transform can also be used.

In this aspect, the forward/reverse time extrapolated deghosted seismic data may be produced by performing any of the appropriate methods of the other aspects of the invention. In a sixth aspect, the invention provides a method of reducing noise in deghosted seismic data, comprising: producing first deghosted seismic data using a model requiring one or more input parameters and seismic data; producing second deghosted seismic data using the same model and seismic data used to produce the first deghosted seismic data but using a different value of at least one of the one or more different input parameter(s);

calculating first impulse response data using the same model and same value(s) of the one or more input parameters as for the first deghosted seismic data; calculating second impulse response data using the same model and same value(s) of the one or more input parameters as for the second deghosted seismic data; using the first deghosted seismic data, the second deghosted seismic data, the first impulse response and the second impulse response to calculate the noise in at least one of the first and second deghosted seismic data.

Deghosted seismic data can contain noise. The present aspect reduces this noise by considering the recorded noise as originating from secondary sources in the model in accordance with Huygens' Principle. It is important to note that, because the first and second deghosted solutions and the first and second impulse responses are produced in the same model, they all have data points at the same locations in the model.

The essence of this aspect is that in each sample of the deghosted seismic data (i.e. the deghosted seismic data at a particular data point) there is an element of noise. In accordance with Huygens' Principle, this element of noise will have an impact on the noise of all data samples at subsequent locations. The aim is therefore to estimate and remove this impact. The above-described aspect achieves this.

Once the impact of the noise from one data point has been calculated and removed from all subsequent data points, there may still be some noise remaining in the subsequent data points. In order to remove the impact of the remaining noise, the same method can be repeated for calculating and removing the impact of the noise remaining in the subsequent data points. The one or more input parameter may be any parameter of the model that the user of the model can input and hence vary. It may preferably be the reflection coefficient, as discussed below. The first/second deghosted seismic data may be produced by any of the appropriate methods disclosed in the various aspects of this application.

An impulse response of the deghosted seismic data is a synthetic simulation of the response of the model due to an input of a spike (such as a delta function). The impulse response data may be calculated by inputting a spike into the model at a location. As is discussed further below, an impulse response may be calculated at each location in the model, by inputting a spike of appropriate amplitude into the model at each location.

Alternatively, one impulse response can be calculated by inputting a spike at a given location in the model (can be any location in the model). The amplitude of this spike may be any amplitude, but preferably 1. The values/amplitude of this impulse response can then be scaled by an appropriate factor for each location in the model, as described further below. The impulse response may be a function of the depth at which the seismic data was captured, the water velocity, the sea surface reflection coefficient and the distance from the spike input location.

Note that for the particularly relevant case where the seismic data capturing depth is constant, the water velocity is constant and the sea surface reflection coefficient does not vary with offset, the impulse response is only a function of distance from the spike input location. The model may advantageously comprise a constant separation between the recording location and the seismic data injection location, a constant water velocity and/or a constant reflection coefficient.

The first and second deghosted seismic data may be calculated separately. The up- going wavefield is identical in each. Thus, only the noise differs between the first and second deghosted seismic data. Thus, taking the difference between the first and second deghosted seismic data results in a gather that contains no desired signal (i.e. the up-going wavefield); but rather only noise. The objective is now to identify which part of this noise came from the first deghosted solution and which part came from the second deghosted solution.

This may be done by calculating a difference seismic data set by calculating the difference between the first and second deghosted solutions. At each location, the contribution to the difference seismic data from each of the first and second deghosted seismic data is known.

Thus, for a given location a first impulse response may be calculated by using the same model and same value(s) of the one or more input parameters as were used for calculating the first deghosted seismic data, by inputting at that location a spike with an amplitude corresponding to the first deghosted seismic data's contribution to the difference seismic data at that location. Also, for that same given location, a second impulse response may be calculated by using the same model and same value(s) of the one or more input parameters as were used for calculating the second deghosted seismic data inputting, by inputting at that location a spike with an amplitude corresponding to the second deghosted seismic data's contribution to the difference seismic data at that location.

Alternatively, the first impulse response may be calculated by using the same model and same value(s) of the one or more input parameters as were used for calculating the first deghosted seismic data, by inputting a spike at any location in the model. The amplitude of this spike may preferably be 1 . For a given location, the first impulse response can be scaled as if it were caused by a spike with an amplitude corresponding to the first deghosted seismic data's contribution to the difference seismic data at that location. The scaled first impulse response can then be treated as if it originated from the given location (i.e. since the value/amplitude of the impulse response is known as a function of distance and time from the spike input location, when calculating the impact of the noise at a given location on all subsequent locations, the scaled first impulse response is treated as if its spike input location is at the given location, and so the value/amplitude of the input response becomes a function of distance and time away the given location). Also, the second impulse response may be calculated by using the same model and same value(s) of the one or more input parameters as were used for calculating the second deghosted seismic data, by inputting a spike at any location in the model. For a given location, the second impulse response can be scaled as if it were caused by a spike with an amplitude corresponding to the second deghosted seismic data's contribution to the difference seismic data at that location. The scaled second impulse response can then be treated as if it originated from the given location (i.e. since the value/amplitude of the impulse response is known as a function of distance and time from the spike input location, when calculating the impact of the noise at a given location on all subsequent locations, the scaled impulse response is treated as if its spike input location is at the given location, and so the value/amplitude of the input response becomes a function of distance and time away from the given location).

Thus, in this case, the first and second impulse responses are calculated only once, but the scaled first and second impulse responses are calculated for each location.

In this case, once the first/second impulse response is calculated, values of the impulse response are known as a function of offset (distance and time) from the spike injection location. These values can be stored in a look-up table. Thus, if these values are scaled appropriately, and subtracted from the corresponding data samples in the first/second deghosted seismic data set (the "corresponding" data samples being those of the

corresponding offset from the "given location" discussed above), noise can be removed from the first/second deghosted seismic data set by having calculated only one impulse response for the model. The first/second impulse response can therefore be calculated by inputting a spike, preferably of amplitude 1 , into the model at any particular location since it is only the value of the impulse response as a function of offset from the input location that is important.

The impulse responses give estimates for the impact of the noise present in the sample at the given location on all subsequent locations.

The (scaled) first impulse response may be subtracted from first deghosted seismic data at subsequent locations. The (scaled) second impulse response may be subtracted from the second seismic data at subsequent locations. This (scaling and) subtraction may occur at each subsequent location in the seismic data. Because the (scaled) first/second impulse response gives an estimate of the impact of the noise present at the current location on all subsequent locations in the first/second deghosted seismic data, when the (scaled) first/second impulse response is subtracted from the deghosted seismic data the impact of the noise from the current location is removed from the first/second deghosted seismic data. If this method is then repeated for all locations subsequent to the given location, then the impact of the noise present in all locations in the first/second deghosted seismic data set is removed. Thus, the noise can be removed from the first/second deghosted seismic data set.

It should be noted that once the impact of the noise at the current sample is removed from all subsequent samples, then a new difference seismic data set may need to be calculated from the new first and second deghosted seismic data sets (i.e. where the impact of noise from the current sample has been removed) before calculating the impact of noise due to noise at a subsequent sample. Doing so ensures that the noise is iteratively removed from the first and second deghosted seismic data sets.

In other words, as can be appreciated, by subtracting an impulse response

(generated by inputting a spike with an appropriate amplitude at one location or scaling a previously calculated impulse response by an appropriate factor calculated from the difference seismic data set) from the first deghosted seismic data, the noise caused by noise in a sample of the first deghosted seismic data at that location can be removed from the first deghosted seismic data. Further, by subtracting an impulse response (generated by inputting a spike with an appropriate amplitude at one location or scaled by an appropriate factor calculated from the difference seismic data set) from the second deghosted seismic data, the noise caused by noise in a sample of the second deghosted seismic data at that location can be removed from the second deghosted seismic data. The calculation of the difference seismic data set, the calculation of (scaled) impulse responses, and the subtraction of the impulse response from the deghosted seismic data, can then be repeated for all other locations in the seismic data. Thus, the method may comprise iteratively subtracting the (scaled) first and/or second impulse response data from the respective first and/or second deghosted seismic data.

Thus, noise-reduced first and/or second deghosted seismic data can be produced.

This noise-reduced first and/or second deghosted seismic data can then be scaled using any appropriate factor. This scaling may be advantageous since it may return the noise-reduced deghosted seismic data to its "correct" value (e.g. correct amplitude). The deghosted seismic data may not be at its "correct" value at this stage because of the prior input parameter manipulation required by the method.

One possible way of reducing the noise in deghosted seismic data is now discussed in more detail. First, the first and second deghosted seismic data sets may be calculated.

Second, the first and second impulse responses may be calculated. The first and second impulse responses are generated by inputting a spike, preferably with an amplitude 1 , into the model and any particular location.

Third, the difference seismic data set may be calculated, as mentioned above. This is done, for instance, by subtracting the first/second deghosted seismic data set from the second/first deghosted seismic data set. As mentioned above, this difference seismic data set is pure noise, or may possibly also include residual ghosts that have not been successfully removed from the deghosted seismic data. As will be appreciated, the difference seismic data set thus comprises numerous data at particular time and space locations. For each sample in the difference seismic data set, a known fraction of the amplitude has come from the first deghosted solution and a known fraction of the amplitude has come from the second deghosted solution. The fractions are known, or can be calculated, from the one or more input parameters that were varied between producing the first and second deghosted seismic data.

Fourth, a data sample of the difference seismic data set may be taken at a first time and space location, preferably the earliest location. A known fraction of the amplitude of this data sample came from the first deghosted solution. A known fraction of the amplitude of this data sample came from the second deghosted solution.

Fifth, the first impulse response may be scaled such that the input spike has an amplitude equal to the known fraction of the data sample that came from the first deghosted solution. Further, the second impulse response may be scaled such that the input spike has an amplitude equal to the known fraction of the data sample that came from the second deghosted solution.

Sixth, the scaled first impulse response can be subtracted from the first deghosted seismic data at all locations including and subsequent to the first location (the spike injection location being aligned with the first location). The scaled second impulse response can be subtracted from the second deghosted seismic data at all locations including and

subsequent to the first location (the spike injection location being aligned with the first location),

These steps can be repeated for all time and space locations. For example, first the steps may be performed over all subsequent spatial locations at the same time location, and then over all subsequent time locations. Alternatively, first the steps may be performed over all subsequent time location at the same spatial location, and then over all subsequent spatial locations. It should be appreciated that once the impact of the noise in the first sample has been removed from the first and second deghosted seismic data sets, then the difference seismic data set should preferably be recalculated using the new first and second deghosted seismic data sets. In this way, when the method is repeated over all subsequent locations, the method iteratively removes the impact of noise from the first and second deghosted seismic data sets.

Once repeated over all data samples (i.e. over all time and space locations), the noise will have been removed from (or at least reduced in) the first and second deghosted data sets.

Thus, the above method may repeat for each spatial location and then for each time interval for which deghosted seismic data has been calculated.

The input parameter may be a reflection coefficient, R, and the first deghosted seismic data may be calculated using R-i and the second deghosted seismic data may be calculated using R₂. The reflection coefficient may be the reflection coefficient applied in step c of the first aspect or the reflection coefficient of the reflecting surface in the second aspect.

In this case, the amplitude of the spike inputted at a given location to calculate the impulse response for the noise present in the first deghosted seismic data, or the amplitude of the spike for the scaled first impulse response, is {Ri (Ri+ R₂) times the amplitude of the difference seismic data sample at the given location}

As discussed above, corresponding resulting data may be calculated for all other time and space locations and then subtracted from the first deghosted seismic data.

As is also discussed above, a corresponding method may be performed for the second deghosted seismic data. In this case, the amplitude of the spike inputted at a given location to calculate the impulse response for the noise present in the second deghosted data, or the amplitude of the spike of for the scaled second impulse response, is {R2 (Ri+ R2) times the amplitude of the difference seismic data sample at the given location}.

The method may further comprise producing noise-reduced seismic data by subtracting the calculated noise (i.e. the impulse response(s)) from the first/second deghosted seismic data.

The method may further comprise rescaling the noise-reduced seismic data using a different reflection coefficient, R₃.

The first/second deghosted seismic data may be produced by performing any of the appropriate methods from the other aspects of the invention. Embodiments of the invention will now be discussed, by way of example only, with reference to the accompanying drawings in which:

Figure 1 shows an example of the model used in the method of the present invention;

Figure 2 shows an exemplary captured data set used to test the present methods; Figure 3 shows reconstructed data based on the data of Figure 2;

Figure 4 shows a portion of the reconstructed data of Figure 3;

Figure 5 shows deghosted reference data for comparison with the deghosted data to be produced by the present method;

Figure 6 shows deghosted data produced by an embodiment of the present invention;

Figure 7 shows deghosted data produced by another embodiment of the present invention;

Figure 8 shows deghosted data produced by an embodiment of the present invention where the noise has been attenuated.

Example of deghosting method

In one embodiment, marine seismic hydrophone data is deghosted by letting the data drive a finite-difference simulation. Similar to Amundsen and Robertsson (2014), by injecting (e.g., Robertsson and Chapman, 2000) seismic data on a finite-difference model it is possible to separate the wavefield into various constituents such as up- and down-going wavefields, the direct wave, multiples, etc. However, whereas Amundsen and Robertsson (2014) require multicomponent seismic data, the present method only requires hydrophone data as well as the knowledge of where the sea surface is relative to the recording locations.

There are several ways that one can configure the model used in this finite-difference simulation to deghost the data. Figure 1 illustrates the basic principles behind the

deghosting method for one such model. For simplicity, assume that data are recorded beneath a flat sea surface at z=0 at a streamer towed at a constant depth z=z_r. However, it should be noted that these assumptions are not necessary and that the method naturally conforms with deghosting hydrophone data acquired under a rough sea surface or with streamers at variable depth. Also assume that the first arrival in the seismic record is an up- going wave. For most marine seismic data this means that the direct wave has been removed from the recorded data which is a common requirement in deghosting of single pressure component data. In a later embodiment, this limiting assumption can be avoided.

The model in Figure 1 has the wave propagation properties of seawater. The top, left and right edges of the model have absorbing boundaries attenuating all waves incident on them. The bottom edge has a Dirichlet boundary condition with monopole point sources coinciding with it. Note that the Dirichlet boundary coincides with the location of the physical streamer in the acquired data at z=z_r. The wavefield emitted from the monopole point sources has two superimposed parts:

1 . the acquired wavefield recorded using hydrophones on the physical streamer at z=z_r; and

2. the wavefield recorded during the time-space finite-difference computation on a mirror streamer at depth z= -z_r, which is injected at z=z_r instantaneously during the finite-difference computation.

Note that the "surface" at z=-z_r is a recording datum and that no boundary conditions are enforced. Also note that there is no free surface at z=0 in the finite-difference computation.

Consider as an example what will happen during the finite-difference computation for a single recorded up-going primary arrival (i.e. the first arrival in the acquired data, which is thus assumed to have only up-going waves, and not down-going waves) and its associated down-going ghost reflection:

1 . First, the primary is injected from the point sources on the physical streamer at z=z_r. This event will radiate upwards with correct amplitude due to the enforced Neumann boundary condition.

2. Next, in the finite-difference computation after some time the primary wavefield reaches z=-z_r where it is recorded as function of offset on the mirror streamer. It will continue to propagate upwards but this is of no interest since it will be absorbed when reaching the top boundary in the finite-difference model.

3. Simultaneously as the wave is recorded along the mirror streamer at z=-z_r it is added instantaneously (scaled with an assumed sign-reversed free-surface reflection coefficient R), time sample for time sample, to the wavefield that is continuously being injected at z=z_r where it will coincide with the ghost that erroneously could be injected upwards from the Dirichlet boundary (if we only have access to hydrophone recordings the injected wavefield will radiate recorded up- and down-going waves in both directions when injected in a finite-difference homogeneous model without vertical component of particle velocity data available).

4. Both the ghost being injected from the acquired data on the physical streamer and the upward propagating wavefield injected at z=z_r from the mirror streamer at z=-z_r are unphysical events if they could travel in the upward direction (i.e. they do not occur in reality, only in the model). The magnitude of their amplitudes are identical but they have opposite polarity since a sign- reversed reflection coefficient R is used in the extrapolation in step 3. As soon as they are injected on the physical streamer, they will destructively interfere.

5. Following the injection of the primary and the ghost on z=z_r there will not be any more energy propagating above z=z_r and therefore we will not be injecting further data that are related to the primary event from the mirror streamer at z=-z_r.

If data are recorded at locations above z=z_r during the finite-difference computation only an up-going primary with the correct arrival time and amplitude will be detected. In other words, the data have been deghosted. Note that the description above was for a primary reflection arriving from the sub-surface but the description applies equally well to multiples.

Finally, it is straightforward to extend this method to cases where the physical streamer has a variable depth. The only difference is that the injection and recording datums at z=z_r and z=-z_r no longer are flat. Instead, the mirror-streamer is positioned as far above the sea surface as the actual streamer is located below it.

Reverse time extrapolation

The deghosting method is based on injecting acquired conventional hydrophone streamer data into finite-difference (FD) grids. In the embodiment described above, the method is implemented using a forward time extrapolator (FTE) method. In this case, it is assumed that the very first energy is up-going as the corresponding down-going wave generated by this arrival is then predicted. In practice this means that the direct wave should be removed.

However, a more robust solution can be achieved by time reversing the data and using a reverse time extrapolator (RTE) instead. This has three significant advantages:

1 . The assumption that the first energy must be up-going is avoided.

2. The deghosting start with the weakest arrivals. Signal-to-noise is better in the RTE implementation compared to the FTE implementation as any noise generated from weaker (deeper) arrivals will interfere with stronger shallower arrivals (deghosting is started using the weakest energy first).

3. It turns out that the FTE and the RTE solutions have identical deghosted signal but different noise. As discussed below separately, noise estimates for the FTE and the RTE solutions from the difference between the FTE and RTE solutions can be predicted. Boundary condition considerations for the model

In a preferred embodiment, the data are injected on a lower boundary in the FD grid onto a Dirichlet boundary condition. The remaining boundaries on the FD grid should be absorbing boundaries so that all energy incident on them vanishes. The inventors found that excellent absorbing boundaries (e.g., convolutional PML boundaries) are critical as spurious reflections otherwise lead to ringing persistent noise in the deghosting solution.

For the side boundaries that are vertically oriented a different boundary condition approach that is almost exact can be used. By computing two solutions once with Dirichlet boundaries and once with Neumann boundaries we can perfectly remove all energy that has interacted with the boundaries once. Note that this is not true for energy that interacts twice with the boundaries but as opposed to general absorbing boundary FD applications, this is a very robust assumption. The fact that energy that has interacted with the boundary once disappears after summation of the two solutions is because Neumann and Dirichlet boundaries are so-called complimentary boundary conditions. This idea turns out to be particularly useful at the edge of the streamer closest to the ship where the source is usually located as energy propagating here has the largest amplitude, very steep grazing angles with respect to the vertical boundary. In addition, this edge is particularly critical due to the missing near offsets that degrades the deghosting solution (see below).

Data reconstruction

The deghosting method predicts down-going energy (i.e. ghosts) from energy that first is up-going. It is therefore essential that the energy producing the ghost has been captured. For the near offsets this can be problematic as data are never recorded down to zero offset in practice (apart from in OBC geometries). In addition, if geology is dipping it may be necessary to also acquired negative offset data. Methods to fill in the missing near offsets are important to mitigate significant noise problems at the near offsets.

Reconstructing missing near offset data prior to deghosting can for instance be done by simple NMO interpolation (although other techniques may be preferred as NMO interpolation may not be accurate enough due to NMO stretching issues).

Reconstructing missing negative offset data prior to deghosting can for instance be done through the method of mirroring positive offset data to the corresponding negative offset locations.

Deghosting stabilization and data conditioning

Apart from inputting the recorded hydrophone data the method relies on the knowledge of the streamer depth profile, the propagation velocity in the water between the streamer and the sea surface as well as the sea surface reflection coefficient. The method tends to be somewhat robust to small errors in these parameters. However, any noise tends to be amplified significantly if a sea surface reflection coefficient with a magnitude close to 1 is used. It is preferable in practice to use a sea surface reflection coefficient magnitudes of 0.95 or 0.9. This is not only to stabilize the deghosting method but also because practical experience shows that such values of the reflection coefficient are more realistic on real data.

Noise suppression methods

As noted above having access to both the FTE and the RTE deghosting solutions have some significant advantages. In particular, it opens up possibilities to estimate and remove noise present in the data after deghosting.

The inventors have found that the predicted up-going energy in the FTE and RTE embodiments of the method is identical. The noise, however, is different. In other words by subtracting the RTE solution from the FTE solution we obtain a gather that contains noise only. Several methods exist to separate the noise corresponding to the RTE and FTE solutions in such a gather. In the following we give one example.

A method for removing noise from the FTE and RTE solutions is based on constructing a gather where every second trace comes from the FTE solution and every second trace comes from the RTE solution. In such a gather the signal is continuous but the noise (if it is uncorrelated) will be unpredictable from trace to trace. By taking a wavenumber transform of such a gather, the noise will appear close to Nyquist as only every second trace is repeatable from the other. The signal on the other hand will map at positive and negative wavenumbers around zero wave number. After the Fourier transform the noise and signal have been separated so that noise (or signal) can be muted and then inverse Fourier transformed for a signal (or noise) estimate.

The same idea can be applied in time as well so that every second sample in a trace come from the RTE solution and every second sample from the FTE solution. The same procedure as for the wavenumber transform procedure is then repeated.

Alternatively, a more sophisticated approach to predict the noise in an FTE or RTE solution can be used. Suppose that there are two solutions from the RTE method but calculated using for instance two different reflection coefficients (the same method can be applied to two FTE solutions as well or an FTE and an RTE solution), e.g. using reflection coefficients of magnitudes 0.9 and 0.95 respectively. The difference between these two solutions results in a pure noise gather since the up-going energy will be identical. What is left is a combination of residual ghosts that have not been successfully removed in the two deghosting solutions plus the recursive impact that this noise has as it propagates and gets fed back through the deghosting method. Next the impulse response from the two deghosting solutions are computed (a synthetic simulation where the input is all zeros except for one spike at a location in the input data).

In summary there are therefore four input data sets:

1 . Deghosted data using reflection coefficient 1.

2. Deghosted data using reflection coefficient 2.

3. Synthetic impulse response for reflection coefficient 1.

4. Synthetic impulse response for reflection coefficient 2.

In order to isolate the two noise contributions the below workflow can be followed where we use the difference gather (the difference between deghosted data using reflection coefficient 1 and the deghosted data using reflection coefficient 2) as input data set:

1 . For each sample at all horizontal locations for time zero of the difference data set there is a value. A factor of 0.95/(0.9+0.95) of the amplitude must have come from deghosting solution 1 and this can therefore add this to a buffer matrix containing the estimate of deghosting solution 1 .

2. The impulse response for deghosting solution 1 due to a spike with the

amplitude 0.95/(0.9+0.95) of the value in the current sample is subtracted from the input difference data set.

3. Likewise the impulse response for deghosting solution 2 due to a spike with the amplitude 0.9/(0.9+0.95) of the value in the current sample is subtracted from the input difference data set.

4. We now proceed to the next spatial location at time 0 and go back to (i)

5. We now step on to t+dt and continue from (i)

In the end an estimate of the error (i.e. noise plus residual ghosts) in deghosting solution 1 (using reflection coefficient 1 ) will have been summed up in the buffer matrix. The output will be a data set that contains up-going waves only of deghosting solution 1 . Since it is up-going waves only, the data can be rescaled so that the result due to any desired sea surface reflection coefficient is obtained, even though for instance a smaller reflection coefficient was used in the deghosting process for stabilization processes.

Example: Application of deghosting method to North Sea data

In the following, applications of the method to one shot from the Mariner survey

(water depth 1 15m) are presented. These data contain recordings along 8 streamers. The two middle streamers (number 4 and 5) are closest in line with the source which is why we focused on these. These data have offsets as close as 60m from the source. The hydrophone data recorded along Streamer 5 are shown in Figure 2. All tests will be applied to these data only. The deghosting method predicts down-going energy (i.e. ghosts) from energy that first is up-going. It is therefore essential that the energy producing the ghost have been captured. For the near offsets this can be problematic as data are never recorded down to zero offset in practice (apart from in OBC geometries). In addition, if geology is dipping it may be necessary to also acquire negative offset data. Methods to fill in the missing near offsets are important to mitigate significant noise problems at the near offsets.

In the example, the missing near offset data were reconstructed through a simple NMO interpolation as follows:

1 . Carry out NMO correction for the closest offset available. An NMO velocity profile as function of two-way travel-time is needed.

2. Copy the trace after moveout to all desired offsets (from zero offset up to the first offset available).

3. Apply an inverse NMO for all the new offsets.

4. The resulting gather now includes data down to zero offset.

5. In addition we then mirrored all data to negative offsets so that we have a symmetric data set around zero offset. This should mitigate all problems at the front of the streamer.

The resulting data are plotted in Figures 3 and 4.

From the survey, there is also reference data showing the desired up-going

(deghosted) wavefield. This is available since the Mariner survey was acquired using a multicomponent streamer enabling deghosting. In contrast the present method does not need multicomponent data and utilizes pressure data only. The reference data are shown in Figure 5. Note that the direct wave and guided wave have been removed in the reference solution (very prominent arrival in a narrow cone in Figure 2). The deghosting estimates will therefore not look the same in this part of the data but this can be ignored for the purpose of demonstrating the current invention.

The first embodiment outlined above assumes that the first arrival is up-going (referred to as the FTE implementation of the deghosting method). For typical towed marine seismic data this is not necessarily the case as the direct wave will be coming from the side or will be down-going. The method assumes that all data are consistent with the wave equation and that up-going waves generate down-going waves. Some care must therefore be taken if muting the data. The inventors have found that in deep water environments it will generally be much easier to identify and remove the direct wave as it separates well in time from the sea bottom reflection and refractions.

However, by inspecting the Mariner P and Z data, the inventors concluded that there were not much left of the direct wave in the data. This is perhaps not a great surprise given the tow depth (17m) and offset to the source (60m or 160m on the first data that we received). The direct wave would be coming directly from the side and suffering badly from interference with its ghost. In practice it can be considered to be gone. The method can therefore proceed by simply not doing anything to the data and by assuming that the first arrival is up-going.

Figure 6 shows a result of applying the FTE deghosting method to the data in Figure

2. First of all the deghosted data mostly resembles what would be expected. In particular, the reflectors between 1.6s and 2.0 s at 2km offset (the reservoir reflections) have been collapsed to narrower events consistent with the removal of the receiver ghost.

A different implementation of the method does not rely on a first arrival assumption: the RTE implementation described above. The result of applying this method to the data is shown in Figure 7. The deghosted signal looks highly repeatable between the FTE and the RTE method (in fact the signal is identical) whereas the noise is different in the two solutions (Figures 6 and 7). There is also close correspondence between the signal in the reference solution in Figure 5.

Although the deghosted solutions appear to have predicted the deghosted signal quite well, there is a large amount of noise in the output. A lot of this noise is different in character between the two deghosting methods. Because the noise between the two solutions is somewhat uncorrelated from trace to trace it is possible to apply the noise attenuation methodology described above that exploits this feature to mitigate significant portions of the noise. Figure 8 shows the result after combining the two deghosting solutions. It should be noted that the method has not resulted in any loss of signal.

Discussion

Beasley et al. (2013a, 2013b) have identified a number of pros and cons for time- space domain approaches to receiver-side deghosting. These pros and cons are valid also for the present deghosting method.

The method can naturally accommodate for streamers with a spatially varying shape. However, by using the method of conformal mapping (Fornberg, 1988) the streamer can be flattened at the expense of an irregular ghost recording location. Precise ghost recording locations are less problematic in the modelling methodology as output data at precise depth locations can be generated using accurate interpolators. Such an approach would allow for the use of sparse finite-difference grids where the grid-spacing is governed by the tolerance of numerical dispersion and not the ability to model the streamer shape.

The method can also be used for source side deghosting. By sorting the data into common receiver gathers and using reciprocity we can run the same methodology on the source side (source locations become receivers where data are acquired). In such a scenario the wavefield will be injected at source locations that now coincide with the Dirichlet boundary. Seismic data are often less well spatially sampled on the source side so care must be taken to ensure that the wavefield is appropriately densely sampled for the finite- difference computation (e.g., by interpolating the wavefield on the source side prior to deghosting).

The Mariner data set is a very challenging data set for the proposed deghosting methods. If such a data set had variable streamer depth that might also help in producing better results due to notch diversity in the data.

References

Amundsen, L, 1993, Wavenumber-based filtering of marine point-source data: Geophysics, 58, 1335-1348, doi: 10.1 190/1 .1443516.

Amundsen, L, and A. Reitan, 2014, Deghosting of marine pressure recordings along curved streamers, EAGE Conference and Exhibition, doi: 10.3997/2214-4609.20141588.

Amundsen, L, and J. O. A. Robertsson, 2014, Wave equation processing using finite- difference propagators, Part I: Wavefield dissection and imaging of marine multicomponent seismic data: Geophysics, in press. Amundsen, L, T. R0sten, J. O. A. Robertsson, and E. Kragh, 2005, On rough-sea deghosting of streamer seismic data using pressure gradient approximations: Geophysics, 70, no. 1 , V1-V9, doi: 10.1 190/1 .1852892.

Amundsen, L, H. Zhou, A. Reitan, and A. B. Weglein, 2013, On seismic deghosting by spatial deconvolution: Geophysics, 78, V267-V271.

Beasley, C. J., R. T. Coates, Y. Ji, and J. Perdomo, 2013a, Wave equation receiver deghosting: A provocative example: SEG Annual Meeting Technical Program Expanded Abstracts, 4226-4230.

Beasley, C. J., R. T. Coates, and C. Lapilli, 2013b, Wave equation receiver deghosting: IEEE 5th international workshop on Computational Advances in Multi-Sensor Adaptive Processing, 280 - 283. Beasley, C, J. and Coates, R. T., 2014, Deghosting measured survey data: US Patent No. US 20140140172 A1. Carlson, D. H., A. Long, W. Sollner, H. Tabti, R. Tenghamn and N. Lunde, 2007, Increased resolution and penetration from a towed dual-sensor streamer: First Break, 25, 71 -77.

Day, A., T. Kluver, W. Sollner, H. Tabti and D. Carlson, 2013, Wavefield separation methods for dual-sensor towed-streamer data: Geophysics, 78, WA55-WA70.

Fokkema, J. T., and P. M. van den Berg, 1993, Seismic applications of acoustic reciprocity: Elsevier Science. Fornberg, B., 1988, The pseudospectral method: Accurate representation of interfaces in elastic wave calculations: Geophysics, 53, 625-637.

Mayhan, J. D., and A. B. Weglein, 2013, First application of Green's theorem-derived source and receiver deghosting on deep-water Gulf of Mexico synthetic (SEAM) and field data: Geophysics, 78, no. 2, WA77-WA89, doi:0.1 190/geo2012-0295.1.

Moldoveanu, N., 2000, Vertical source array in marine seismic exploration: 70th Annual International Meeting, SEG, Expanded Abstracts, 53-56. Moldoveanu, N., 2001 , A seismic source, a marine seismic surveying arrangement, a method of operating a marine seismic source, and a method of deghosting seismic data: Australian Patent WO 01/75481 .

Moldoveanu, N., L. Combee, M. Egan, G. Hampson, L. Sydora, and W. Abriel, 2007, Over/under towed streamer acquisition: A method to extend seismic bandwidth to both higher and lower frequencies: The Leading Edge, 26, 41-58, doi: 10.1 190/1 .2431831 .

Ozbek, A., M. Vassallo, K. Eggenberger, D. J. van Manen, K. Ozdemir, S. Zeroug and J. O. A. Robertsson, 2013, New Marine Towed-streamer Acquisition Technology: OTC24001 , OTC conference, Houston.

Ozdemir, A. K., P. Caprioli, A. Ozbek, E. Kragh and J. O. A. Robertsson, 2008, Optimized deghosting of over/under towed streamer data in the presence of noise: The Leading Edge, 27, 190-199.

Parkes, G., and S. Hegna, 201 1 , A marine seismic acquisition system that provides a full "ghost-free" solution: 81 st Annual International Meeting, SEG, Expanded Abstracts, 37-41 . Riyanti, C. D., R. G. van Borselen, P. M. van den Berg, and J. T. Fokkema, 2008, Pressure wave-field deghosting for non-horizontal streamers: SEG

Robertsson, J. O. A. and L. Amundsen, 2014, Prediction of wavefield constituents by modelling, Part II: Deghosting of marine hydrophone streamer data in the time-space domain: Geophysics, in press.

Robertsson, J. O. A. and J. O. Blanch, 201 1 , Numerical methods, finite difference, in H. Gupta, ed., Encyclopaedia of Solid Earth Geophysics (ESEG), 2nd Edition: Springer Verlag, pp. 883-892.

Robertsson, J. O. A. and C. H. Chapman, 2000, An efficient method for calculating finite- difference seismograms after model alterations: Geophysics, 65, 907 - 918. Robertsson, J. O. A., and E. Kragh, 2002, Rough sea deghosting using a single streamer and a pressure gradient approximation: Geophysics, 67, 2005-201 1 .

Robertsson, J. O. A., I. Moore, M. Vassallo, A. K. Ozdemir, D. J. van Manen and A. Ozbek, 2008, On the use of multicomponent streamer recordings for reconstruction of pressure wavefields in the crossline direction: Geophysics, 73, A45-A49.

Robertsson, J. O. A., D.-J. van Manen, D. Halliday, and R. Laws, 201 1 , Seismic data acquisition and source-side derivatives generation and application: U. S. Patent 7,876,642. Soubaras, R. and Y. Lafet, 2013, Variable-depth streamer acquisition: Broadband data for imaging and inversion, 78, WA27-WA39. doi: 10.1 190/geo2012-0297.1

Vaage, S. T., 2005, Method and system for acquiring marine seismic data using multiple seismic sources: U. S. Patent 6,906,981. van Borselen, R. G., J. Fokkema, and P. van den Berg, 2013, Wavefield decomposition based on acoustic reciprocity: Theory and applications to marine acquisition: Geophysics, 78, WA41 -WA54. Virieux, J., 1984, SH-wave propagation in heterogeneous media: Velocity-stress finite- difference method: Geophysics, 49, 1933-1957. Weglein, A. B., S. A. Shaw, K. H. Matson, J. L. Sheiman, R. H. Stolt, T. H. Tan, A. Osen, G. P. Correa, K. A. Innanen, Z. Guo, and J. Zhang, 2002, New approaches to deghosting towed-streamer and ocean-bottom pressure measurements: 72nd Annual International Meeting, SEG, Expanded Abstracts, 1016-1019.

Claims

Claims:

A method of deghosting seismic data using a model, the seismic data having been captured at a location below the sea surface and the model comprising a wave- propagation space having the wave propagation properties of sea-water, the method comprising:

a. injecting into the model a seismic wavefield based on the seismic data at an injection location in the model corresponding to the location at which the seismic data was captured, and allowing the injected seismic wavefield to propagate through the wave-propagation space;

b. recording the propagating seismic wavefield at a recording location in the model, the recording location being located at a distance from the injection location of twice the depth of the location below the sea surface at which the seismic data was captured;

c. scaling or convolving the recorded seismic wavefield with a reflection coefficient, R;

d. adding the scaled or convolved seismic wavefield to a subsequent

seismic wavefield, which is based on a corresponding subsequent portion of the seismic data, the subsequent seismic wavefield being the seismic wavefield that is to be injected at the injection location at the same time as the propagating seismic wavefield passes the recording location.

A method as claimed in claim 1 , wherein the recording location is vertically above or vertically below the injection location in the model.

A method as claimed in claim 1 or 2, wherein the magnitude of reflection coefficient, R, is between 0 and 2, preferably between 0 and 1.

A method as claimed in any preceding claim, wherein the reflection coefficient, R, is a function that may vary with time.

A method as claimed in any preceding claim, wherein the method of conformal mapping is used to model the effects of a rough sea surface and/or a depth varying recording streamer.

A method as claimed in any preceding claim, wherein the wave-propagation space comprises edges, at least some of which have absorbing boundary conditions.

7. A method as claimed in any preceding claim, wherein the injection location is located at an edge portion of the wave-propagation space.

8. A method as claimed in claim 7, wherein the edge portion of the model comprises a Dirichlet boundary condition.

9. A method as claimed in claim 8, wherein at least some of the remaining edges of the wave-propagation space have absorbing boundary conditions. 10. A method as claimed in any of claims 1 to 6, wherein the injection location is within the wave-propagation space.

1 1 . A method as claimed in claim 10, wherein the at least two of the edges of the wave- propagation space have absorbing boundary conditions.

12. A method as claimed in claim 6, 9 or 1 1 , wherein the absorbing boundary conditions are Perfectly Matched Layers.

13. A method as claimed in any preceding claim, wherein the wave-propagation space comprises two side edge portions opposite one another, and the method comprises performing steps a to d with at least one of the two side edge portions of the model comprising a first boundary condition to provide first deghosted seismic data;

performing steps a to d with the at least one of the two side edge portions of the model comprising a second complimentary boundary condition to provide second deghosted seismic data; and summing the first and second deghosted seismic data.

14. A method as claimed in any preceding claim, wherein the seismic data consists of pressure data. 15. A method as claimed in any of claims 1 to 13, wherein the seismic data consists of particle velocity data.

16. A method as claimed in any preceding claim, wherein the method is a finite- difference method, a finite-element method, a spectral element method, or a pseudo- spectral method.

17. A method as claimed in any preceding claim, wherein the modelling is in the space- time domain.

18. A method as claimed in any preceding claim, wherein the seismic wavefield is injected at a point source or a plurality of point sources.

19. A method as claimed in claim 18, wherein respective point sources are at the same depth as one another.

20. A method as claimed in claim 18, wherein respective point sources are at varying depths with respect to one another.

21 . A method as claimed in any of claims 18 to 20, wherein the propagating seismic wave is recorded a plurality of locations.

22. A method as claimed in any preceding claim, wherein the model does not include any free surface.

23. A method as claimed in any preceding claim, wherein the direct wave has been

removed from the seismic data.

24. A method as claimed in any preceding claim, wherein the seismic data comprises captured seismic data, the captured seismic data missing near offset seismic data, and reconstructed near offset seismic data, the method comprising reconstructing the reconstructed near offset seismic data.

25. A method as claimed in any preceding claim, wherein the seismic data comprises captured seismic data and reconstructed negative offset seismic data, the method comprising reconstructing the reconstructed negative offset seismic data.

26. A method as claimed in claim 25, wherein the reconstructed negative offset seismic data is reconstructed by mirroring positive offset seismic data.

27. A method as claimed in any preceding claim, comprising performing steps a to d using forward time extrapolation to produce forward time extrapolated seismic data.

28. A method as claimed in any of claims 1 to 26, comprising performing steps a to d using reverse time extrapolation to produce reverse time extrapolated seismic data.

29. A method of reconstructing a multicomponent wavefield from seismic data using a model, the seismic data having been captured at a location below the sea surface, and the model comprising a wave-propagation space having the wave propagation properties of sea-water and a reflecting surface at a location corresponding to the location of the sea surface, the method comprising:

a. injecting into the model a seismic wavefield based on the seismic data at an injection location in the model corresponding to the location at which the seismic data was captured;

b. introducing a boundary condition at the injection location; and

c. allowing the wavefield to propagate through the wave-propagation space and reflect off the reflecting surface, thus reconstructing the total wave-field.

30. A method as claimed in claim 29, wherein the seismic data consists of pressure data. 31 . A method as claimed in claim 29, wherein the seismic data consists of particle

velocity data.

32. A method as claimed in claim 29, 30 or 31 , wherein the boundary condition is

configured such that the injected seismic wavefield only propagates generally toward the reflecting surface.

33. A method as claimed in claim 29, 30, 31 or 32, wherein the boundary condition has a reflection coefficient of -1 . 34. A method as claimed in any of claims 29 to 33, wherein the boundary condition is a

Dirichlet boundary condition.

35. A method as claimed in any of claims 29 to 33, wherein the boundary condition is a Neumann boundary condition.

36. A method as claimed in any of claims 29 to 35, wherein the reflection coefficient of the reflecting surface is -1 < R < 0, preferably -1 , further preferably -0.95.

37. A method as claimed in any of claims 29 to 36, wherein the reflection coefficient of the reflecting surface is a function that may vary with time.

38. A method as claimed in any of claims 29 to 37, wherein conformal mapping is used to model the effects of a rough sea surface and/or a depth varying recording streamer.

39. A method as claimed in any of claims 29 to 38, wherein the injection location is located at an edge portion of the wave-propagation space.

40. A method as claimed in any of claims 29 to 40, wherein the wave-propagation space comprises two side edge portions opposite one another, and the method may comprise performing steps a to c with at least one of the two side edge portions of the model comprising a first boundary condition to provide a first reconstructed wavefield; performing steps a to c with the at least one of the two side edge portions of the model comprising a second complimentary boundary condition to provide a second reconstructed wavefield; and summing the first and second reconstructed wavefields.

41 . A method as claimed in any of claims 29 to 40, comprising:

d. deghosting the seismic data using the total reconstructed wave-field.

42. A method as claimed in any of claims 29 to 41 , wherein steps a to c are performed using forward time extrapolation.

43. A method as claimed in any of claims 29 to 41 , wherein steps a to c are performed using reverse time extrapolation.

44. A method of deghosting single-component seismic data comprising using a reverse time extrapolator.

45. A method as claimed in claim 44, wherein the single-component seismic data

consists of pressure data.

46. A method as claimed in claim 44, wherein the seismic data consists of particle

velocity data. 47. A method as claimed in claim 28, 44, 45 or 46, comprising calculating a down-going wavefield from the single-component seismic data using the reverse time

extrapolator and subtracting the down-going wavefield from the total wavefield.

48. A method as claimed in any of claims 44 to 47, comprising using a model having at least one edge, the at least one edge having absorbing boundary conditions.

49. A method as claimed in any of claims claim 44 to 48, comprising using a model having at least two edges; producing first deghosted seismic data using the single- component seismic data and a first set of boundary conditions for at least one of the two edges; producing second deghosted seismic data using the same single- component seismic data and a second set of boundary conditions for the at least one of the edges, the first and second sets of boundary conditions being complimentary to each other; and summing the first and second deghosted seismic data. 50. A method as claimed in any of claims 44 to 49, wherein the direct wave has been removed from the single-component seismic data.

51 . A method as claimed in any of claims 44 to 50, wherein the single-component

seismic data comprises captured single-component seismic data, the captured single-component seismic data missing near offset single-component seismic data, and reconstructed near offset single-component seismic data, the method comprising reconstructing the reconstructed near offset single-component seismic data.

52. A method as claimed in any of claims 44 to 51 , wherein the single-component

seismic data comprises captured single-component seismic data and reconstructed negative offset single-component seismic data, the method comprising reconstructing the reconstructed negative offset single-component seismic data.

53. A method as claimed in claim 52, wherein the reconstructed negative offset single- component seismic data is reconstructed by mirroring positive offset single- component seismic data. 54. A method as claimed in any of claims 44 to 53, comprising any of the features of claims 1 to 26, claims 28 to 41 or claim 43.

55. A method comprising subtracting forward time extrapolated deghosted seismic data from reverse time extrapolated deghosted seismic data, or subtracting the reverse time extrapolated deghosted seismic data from the forward time extrapolated deghosted seismic data.

56. A method as claimed in claim 55, wherein the method is a method of identifying noise in deghosted seismic data.

57. A method as claimed in claim 55 or 56, wherein the method is a method of reducing the noise in deghosted seismic data.

58. A method of reducing noise in deghosted seismic data, comprising:

constructing a gather of deghosted seismic data by producing forward time extrapolated deghosted seismic data for a first set of locations and producing reverse time extrapolated seismic data for a second set of locations, wherein the first and second set of locations are interleaving locations,

performing a transform on the gather;

muting the noise signal on the transform; and

performing the inverse transform on the transform.

59. A method as claimed in claim 58, wherein the locations are temporal or spatial.

60. A method as claimed in claim 58 or 59, wherein the transform is a Fourier transform.

61 . A method as claimed in claim 58 or 59, wherein the transform is a Radon transform.

62. A method as claimed in any of claims 58 to 61 , wherein producing the forward time extrapolated deghosted seismic data comprises performing the method of claim 26 or 41 and wherein producing the reverse time extrapolated deghosted seismic data comprises performing the method of any of claims 27 or 42 to 51.

63. A method of reducing noise in deghosted seismic data, comprising: producing first deghosted seismic data using a model requiring one or more input parameters and seismic data;

producing second deghosted seismic data using the same model and seismic data used to produce the first deghosted seismic data but using a different value of at least one of the one or more different input parameter(s);

calculating first impulse response data using the same model and same value(s) of the one or more input parameters as for the first deghosted seismic data;

calculating second impulse response data using the same model and same value(s) of the one or more input parameters as for the second deghosted seismic data; using the first deghosted seismic data, the second deghosted seismic data, the first impulse response and the second impulse response to calculate the noise in at least one of the first and second deghosted seismic data.

64. A method as claimed in claim 63, comprising calculating difference seismic data by subtracting the first/second deghosted seismic data from the second/first seismic data.

65. A method as claimed in claim 64, comprising scaling the first and second impulse response data and subtracting the scaled first impulse response data and the scaled second impulse response data from the difference seismic data.

66. A method as claimed in any of claims 63 to 65, wherein the input parameter is a reflection coefficient, R, and wherein the first deghosted seismic data is calculated using R-i and the second deghosted seismic data is calculated using R₂.

67. A method as claimed in claim 66, wherein the impact of the noise in one sample in time and space is calculated from the first impulse response data scaled such that the input spike of the first impulse response has an amplitude of Ri ( Ri+ R₂) times the value of the difference seismic data at that time and the same space location.

68. A method as claimed in claim 66 or 67, wherein the impact of the noise in one

sample in time and space is calculated from the second impulse response data scaled such that the input spike of the second impulse response has an amplitude of R₂/( Ri+ R₂) times the value of the difference seismic data at the same time and the same spatial location.

69. A method as claimed in any of claims 63 to 68, comprising producing noise-reduced seismic data by subtracting the calculated noise from the first/second deghosted seismic data.

70. A method as claimed in claim 69, when dependent on any of claims 66 to 68,

comprising rescaling the noise-reduced seismic data using a different reflection coefficient, R₃.

71 . A method as claim in any of claims 63 to 70, wherein producing the first deghosted seismic data comprises performing the method of any of claims 1 to 54 and producing the second deghosted seismic data comprises performing the method of any of claims 1 to 54.