Method for deghosting seismic data

Abstract

A method of at least partially deghosting recorded seismic s-waves, wherein recorded seismic data is provided, wherein the recorded seismic data has been recorded at a receiver located beneath the Earth's surface, and wherein the recorded seismic data includes s-wave data. The method may include the steps of finding a model of the Earth's crust for use in deghosting the recorded seismic data using the s-wave data, wherein the model includes at least one region and wherein the model includes the Earth's surface and the location of the receiver, using the model to find a deghosting operator that, when applied to the s-wave data, at least partially deghosts the s-wave data, and applying the deghosting operator to the s-wave data to at least partially deghost the s-wave data.

Claims

1. A method of at least partially deghosting recorded seismic s-waves, wherein recorded seismic data is provided, wherein said recorded seismic data is recorded at a receiver located beneath Earth's surface, and wherein said recorded seismic data comprises s-wave data, the method comprising: finding a model of the Earth's crust for use in deghosting the recorded seismic data using the s-wave data, wherein the model comprises at least one region, and wherein the model comprises the Earth's surface and the location of the receiver; using said model to find a deghosting operator that, when applied to the s-wave data, at least partially deghosts the s-wave data; and applying the deghosting operator to the s-wave data to at least partially deghost the s-wave data.

2. The method as claimed in claim 1, wherein the model comprises at least two regions, and wherein the model comprises a boundary between the two regions.

3. The method as claimed in claim 1, further comprising at least one of: recording the recorded seismic data; and selecting substantially only the s-wave data from the recorded seismic data.

4. The method as claimed in claim 1, wherein finding the model comprises at least one of: determining at least one of number of regions in the model, at least size of one of the number regions in the model and the location of the receiver in the model; finding the model using the s-wave data directly; assuming that all of the s-wave data is at least substantially purely vertically-propagating; and finding the model using geotechnical data, wherein geotechnical data is provided.

5. The method as claimed in claim 1, wherein the model comprises at least three regions.

6. The method as claimed in claim 1, wherein the model is at least partially defined by parameters comprising an s-wave reflection coefficient of the Earth's surface, the s-wave reflection coefficient at any boundaries or boundary between adjacent regions, an s-wave velocity, an attenuation factor and/or density.

7. The method as claimed in claim 6, wherein finding the model comprises determining the parameters.

8. The method as claimed in claim 1, wherein second recorded seismic data is provided, wherein the second recorded seismic data comprises second s-wave data, and wherein finding the model comprises finding the model using a comparison of the s-wave data and the second s-wave data.

9. The method as claimed in claim 1, wherein using said model to find said deghosting operator comprises finding an impulse response of the model, and finding the deghosting operator from the impulse response.

10. The method as claimed in claim 1, further comprising using the model to find an operator to bring the s-wave data from the receiver location to the Earth's surface.

11. The method as claimed in claim 1, wherein the model is an elastic model or a visco-elastic model.

12. The method as claimed in claim 1, wherein at least one of the s-wave velocity of the model and quality factor of the model generally decreases toward the Earth's surface.

13. The method as claimed in claim 1, wherein the recorded seismic data and the s-wave data comprises data from one or more sources having large offsets from the receiver.

14. The method as claimed in claim 13, wherein the large offsets from the receiver are greater than 500 m, 1000 m, 3000 m, or 7000 m.

15. The method as claimed in claim 1, wherein the recorded seismic data has been recorded at a receiver beneath in sea bed.

16. The method as claimed in claim 15, wherein the model is a model of the sea bed.

17. The method of at least partially deghosting recorded seismic s-waves, wherein recorded seismic data recorded at a plurality of receivers is provided, wherein said plurality of receivers are located beneath the Earth's surface, and wherein said recorded seismic data comprises s-wave data, the method comprising performing the method of claim 1 for each of the receivers.

18. The method of imaging a geological structure comprising performing the method of claim 1 to produce at least partially deghosted s-wave data, and using said at least partially deghosted s-wave data to image the geological structure.

19. The method as claimed in claim 18, comprising repeating the method of claim 18 to produce images of the geological structure at different times; and using said images to view or determine how the geological structure is changing over time.

20. A computer program product comprising computer readable instructions that, when run on a computer, is configured to perform the method of claim 1.

21. A method of prospecting for hydrocarbons, comprising: performing the method of claim 1; and using the at least partially deghosted s-wave data to prospect for hydrocarbons.

22. A method of producing hydrocarbons, comprising: performing the method of claim 1; and producing hydrocarbons.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) Certain preferred embodiments will now be described by way of example only and with reference to the accompanying drawings, in which

(2) FIG. 1 shows a comparison between an image of a reservoir produced from data measured at the seabed (FIG. 1b) and data measured beneath the seabed (FIG. 1a);

(3) FIG. 2 shows a schematic workflow of the present method;

(4) FIG. 3 shows an example of data recorded at the seabed (FIG. 3a) and a dataset recorded by a buried receiver (FIG. 3b);

(5) FIG. 4 shows an example of the step of comparing the spectrum of recorded seismic data from a buried receiver to the spectrum of recorded seismic data from the seabed;

(6) FIG. 5 shows an example of a model used in the present method;

(7) FIG. 6 shows another example of a different model used in the present method;

(8) FIG. 7 shows a comparison between the magnitude part of the impulse response of a model and the recorded seismic data;

(9) FIG. 8 shows how recorded seismic data in the time domain can be improved using the present method; s

(10) FIG. 9 shows how an image of a reservoir can be improved using the present method;

(11) FIG. 10 illustrates a horizontally layered model with three layers;

(12) FIG. 11 shows an example of synthetic response at the receiver level and at the seabed (surface); and

(13) FIG. 12 shows sections of final migrated PS full stacks from the same location with different processing performed.

DETAILED DESCRIPTION

(14) With regard to FIG. 1, FIG. 1a shows an image of a reservoir produced using ps-wave seismic data measured using receivers that are buried beneath the seabed. FIG. 1b shows an image of the same reservoir produced using ps-wave seismic data measured using receivers that are placed on the seabed.

(15) The inventors observed that the image of the reservoir produced using buried receivers is more blurred and distorted than the image produced using receivers on the seabed. Further, the inventors observed that the subsurfaces are not as continuous for the image produced using buried receivers as it is for the image produced using receivers on the seabed.

(16) The inventors have investigated the causes of this worse image quality, and have developed the present method to improve the quality of the imaging using ps-waves recorded at buried receivers.

(17) FIG. 2 shows a schematic workflow of the present invention.

(18) The first step 1 is to gather seismic data from a seismic receiver that is buried beneath the Earth's surface. This data comprises horizontal particle velocity and may also comprise vertical particle velocity and/or pressure data.

(19) The second step 2 is to filter the recorded seismic data such that the resulting filtered data is substantially only ps-wave data. One way this is achieved is shown in FIG. 3. First of all, the horizontal particle data only is considered. The inventors have found that, because ps-waves are very slow (in comparison to p-waves), the ps-waves incident on the buried receiver are substantially only vertically-propagating. This means that the ps-waves are substantially only recorded on the horizontal particle velocity receiver. Secondly, the p-wave contribution to the horizontal particle velocity data that is recorded is removed. This is done by removing all data that arrives before a certain time for a given source-receiver offset. Because p-waves are faster than ps-waves, the initial part of the data will contain p-wave contributions and no or little ps-wave contributions. The amount of data to be ignored (or removed or filtered) can be determined by comparing the horizontal particle velocity data with the vertical particle velocity data and/pressure data, which will include corresponding data from the p-wave arrivals. In FIG. 3, the data from times before the dotted line is removed from the data set.

(20) Returning to FIG. 2, step 3 involves producing a model 10 for deghosting the recorded ps-wave data from the recorded ps-wave data. Two different examples of such a model are shown in FIGS. 5 and 6.

(21) Regarding FIG. 5, the model 10 is a 1D model that consists of two layers 11, 12. The layers 11, 12 comprise an upper layer 11 and a lower layer 12. There is a boundary 14 between the upper layer 11 and the lower layer 12. The model 10 also comprises the Earth's surface 13 on top of the upper layer 11. The model comprises a receiver location 15. The boundary 14 is at a depth d0 and the receiver 15 is at a depth d1. In FIG. 5, d1 is shown as being greater than d0, but d1 may be less than d0. The depth of the model extends to a depth below the receiver 15 depth d1 and below the boundary 14 depth d0. The ghost signal at the receiver 15 may include ringing from multiple reflections between boundaries 13 and 14.

(22) The upper layer 11 may have a depth d0 of greater than 0.05 m, greater than 0.2 m, or greater than 0.5 m. The upper layer 11 may have a depth d0 of less than 2 m, less than 1 m, or less than 0.5 m. The uppermost layer 11 may have a depth d0 of 0.05-3 m, 0.2-2 m or 0.5-1 m.

(23) The upper layer 11 may have an s-wave velocity of greater than 5 m/s, greater than 20 m/s or greater than 50 m/s. The upper layer 11 may have an s-wave velocity of less than 300 m/s, less than 100 m/s or less than 50 m/s. The uppermost layer 11 may have an s-wave velocity of 5-300 m/s or 5-100 m/s.

(24) The upper layer 11 may have a quality factor of greater than 1, greater than 2 or greater than 5. The uppermost layer 11 may have a quality factor of less than 2000, less than 1000, or less than 500. The upper layer 11 may have a quality factor of 1-2000, 2-1500 or 5-1000.

(25) The upper layer 11 may have a density of greater than 500 kg/m.sup.3, greater than 1000 kg/m.sup.3, or greater than 2000 kg/m.sup.3. The uppermost layer 11 may have a density of less than 10000 kg/m.sup.3, less than 5000 kg/m.sup.3, or less than 2000 kg/m.sup.3. The uppermost layer 11 may have a density of 500-10000 kg/m.sup.3, 500-5000 kg/m.sup.3, or 500-2000 kg/m.sup.3.

(26) The lower layer 12 may have an s-wave velocity of greater than 20 m/s, greater than 50 m/s or greater than 100 m/s. The lower layer 12 may have an s-wave velocity of less than 600 m/s, less than 200 m/s or less than 100 m/s. The lowermost layer 12 may have an s-wave velocity of 20-600 m/s, 50-200 m/s or 70-150 m/s. The lower layer 12 may have a greater s-wave velocity than the upper layer 11.

(27) The lower layer 12 may have a quality factor of greater than 2, greater than 5, greater than 10 or greater than 100. The lower layer 12 may have a quality factor of less than 10000, less than 5000, or less than 1000. The lower layer 12 may have a quality factor of 2-10000, 10-5000 or 50-2000. The quality factor of the lower layer 12 may be greater than the quality factor of the upper layer 11.

(28) The lower layer 12 may have a density of greater than 500 kg/m.sup.3, greater than 1000 kg/m.sup.3, or greater than 2000 kg/m.sup.3. The lower layer 12 may have a density of less than 10000 kg/m.sup.3, less than 5000 kg/m.sup.3, or less than 2000 kg/m.sup.3. The lower layer 12 may have a density of 500-10000 kg/m.sup.3, 500-5000 kg/m.sup.3, or 500-2000 kg/m.sup.3. The lower layer 12 may have the same density as the upper layer 11.

(29) The reflection coefficient of the boundary 14 between the upper 11 and lower 12 layers may be at least partially defined by the s-wave velocity in those layers.

(30) The reflection coefficient of the Earth's surface 13 may be 0.1-1, preferably 0.4-1, preferably 0.6-1, preferably 0.6-0.9.

(31) The receiver 15 may have a depth d1 of greater than 0.1 m, greater than 0.2 m, greater than 0.5 m, greater than 1 m or greater than 1.5 m. The receiver 15 may have a depth d1 of less than 4 m, less than 2 m or less than 1 m. The receiver 15 may have a depth d1 of 0.1-4 m, 0.2-2 m or 0.5-1.5 m. The depth d1 of the receiver 15 may be greater than or less than the depth d0 of the upper layer 11.

(32) Regarding FIG. 6, the model 10 is a 1D model that consists of three layers 11, 12, 16. The layers 11, 12, 16 comprise an upper layer 11, an intermediate layer 12 and a lower layer 16. There is an upper boundary 14 between the upper layer 11 and the intermediate layer 12, and there is a lower boundary 17 that is between the intermediate layer 12 and the lower layer 16. The model 10 also comprises the Earth's surface 13 on top of the upper layer 11. The model comprises a receiver location 15. The upper boundary 14 is at a depth d0, the receiver 15 is at a depth d1 and the lower boundary is at depth d2. In FIG. 6, d1 is shown as being greater than d0 but less than d2, but d1 may be less than d0 or greater than d2. The depth of the model extends to a depth below the receiver 15 depth d1 and below the boundary 17 depth d2. Because the receiver 15 may be located between two reflective boundaries 14, 17, the ghost signal at the receiver 15 may include ringing from multiple reflections between boundaries 14 and 17 in addition to multiple reflections between boundaries 13 and 14 and combinations thereof.

(33) The upper layer 11 may have a depth d0 of greater than 0.05 m, greater than 0.2 m, or greater than 0.5 m. The upper layer 11 may have a depth d0 of less than 2 m, less than 1 m, or less than 0.5 m. The uppermost layer 11 may have a depth d0 of 0.05-3 m, 0.2-2 m or 0.5-1 m.

(34) The upper layer 11 may have an s-wave velocity of greater than 5 m/s, greater than 20 m/s or greater than 50 m/s. The upper layer 11 may have an s-wave velocity of less than 200 m/s, less than 100 m/s or less than 50 m/s. The uppermost layer 11 may have an s-wave velocity of 5-200 m/s or 20-100 m/s.

(35) The upper layer 11 may have a quality factor of greater than 1, greater than 2 or greater than 5. The uppermost layer 11 may have a quality factor of less than 2000, less than 1000, or less than 500. The upper layer 11 may have a quality factor of 1-2000, 2-1500 or 5-1000.

(36) The upper layer 11 may have a density of greater than 500 kg/m.sup.3, greater than 1000 kg/m.sup.3, or greater than 2000 kg/m.sup.3. The uppermost layer 11 may have a density of less than 10000 kg/m.sup.3, less than 5000 kg/m.sup.3, or less than 2000 kg/m.sup.3. The uppermost layer 11 may have a density of 500-10000 kg/m.sup.3, 500-5000 kg/m.sup.3, or 500-2000 kg/m.sup.3.

(37) The intermediate layer 12 may have a depth d2 of greater than 0.2 m, greater than 0.4 m, or greater than 1 m. The intermediate layer 12 may have a depth d2 of less than 10 m, less than 5 m, or less than 2 m. The intermediate layer 12 may have a depth d2 of 0.2-10 m or 1-5 m.

(38) The intermediate layer 12 may have an s-wave velocity of greater than 20 m/s, greater than 50 m/s, or greater than 100 m/s. The intermediate layer 12 may have an s-wave velocity of less than 400 m/s, less than 200 m/s or less than 100 m/s. The intermediate layer 12 may have an s-wave velocity of 20-400 m/s, 50-200 m/s or 70-150 m/s. The intermediate layer 12 may have a greater s-wave velocity than the upper layer 11.

(39) The intermediate layer 12 may have a quality factor of greater than 2, greater than 5, greater than 10 or greater than 100. The intermediate layer 12 may have a quality factor of less than 3000, less than 2000, or less than 1000. The intermediate layer 12 may have a quality factor of 2-3000, 10-2000 or 50-1000. The quality factor of the intermediate layer 12 may be greater than the quality factor of the upper layer 11.

(40) The intermediate layer 12 may have a density of greater than 500 kg/m, greater than 1000 kg/m.sup.3, or greater than 2000 kg/m.sup.3. The intermediate layer 12 may have a density of less than 10000 kg/m.sup.3, less than 5000 kg/m.sup.3, or less than 2000 kg/m.sup.3. The intermediate layer 12 may have a density of 500-10000 kg/m.sup.3, 500-5000 kg/m, or 500-2000 kg/m.sup.3. The intermediate layer 12 may have the same density as the upper layer 11.

(41) The reflection coefficient of the boundary 14 between the upper 11 and intermediate 12 layers may be at least partially defined by the s-wave velocity in those layers.

(42) The lower layer 16 may have an s-wave velocity of greater than 40 m/s, greater than 100 m/s or greater than 200 m/s. The lower layer 16 may have an s-wave velocity of less than 1000 m/s, less than 600 m/s or less than 500 m/s. The lower layer 16 may have an s-wave velocity of 40-1000 m/s, 100-600 m/s or 200-500 m/s. The lower layer 16 may have a greater s-wave velocity than the upper layer 11 and the intermediate layer 12.

(43) The lower layer 16 may have a quality factor of greater than 4, greater than 10, greater than 20 or greater than 100. The lower layer 16 may have a quality factor of less than 5000, less than 1000, or less than 500. The lower layer 16 may have a quality factor of 4-5000, 20-1000 or 100-500. The quality factor of the lower layer 16 may be greater than the quality factor of the upper layer 11 and the intermediate layer 12.

(44) The lower layer 16 may have a density of greater than 500 kg/m.sup.3, greater than 1000 kg/m.sup.3, or greater than 2000 kg/m.sup.3. The lower layer 16 may have a density of less than 10000 kg/m.sup.3, less than 5000 kg/m.sup.3, or less than 2000 kg/m.sup.3. The lower layer 16 may have a density of 500-10000 kg/m.sup.3, 500-5000 kg/m.sup.3, or 500-2000 kg/m.sup.3. The lower layer 16 may have the same density as the upper layer 11 and the intermediate layer 12.

(45) The reflection coefficient of the boundary 17 between the lower 16 and intermediate 12 layers may be at least partially defined by the s-wave velocity in those layers.

(46) The reflection coefficient of the Earth's surface 13 may be 0.1-1, preferably 0.4-1, preferably 0.6-1, preferably 0.6-0.9.

(47) The receiver 15 may have a depth d1 of greater than 0.1 m, greater than 0.2 m, greater than 0.5 m, greater than 1 m or greater than 1.5 m. The receiver 15 may have a depth d1 of less than 4 m, less than 2 m or less than 1 m. The receiver 15 may have a depth d1 of 0.1-4 m, 0.2-2 m or 0.5-1.5 m. The depth d1 of the receiver 15 may be greater than or less than the depth d0 of the upper layer 11, and greater than or less than the depth d2 of the intermediate layer 12.

(48) FIGS. 5 and 6 are just exemplary models, and others with differing layer numbers could also be used.

(49) The first step in producing such a model 10 is determining the number of layers 11, 12, 16 that should be present in the model 10. This may be achieved by the user selecting the number of layers 11, 12, 16 in the model 10 and by inverting the s-wave data. Once the number of layers 11, 12, 16 has been selected, the parameters of the model 10 can be found (see more on this below) by inverting the s-wave data, and then the impulse response of the model 10 can be found. If this impulse response does not accurately model the observed ghosting in the recorded s-wave data, then a different number of layers 11, 12, 16 can be selected, and the process can start again.

(50) The next step in finding the model 10 comprises determining the depth d0, d2 to the base of the layer(s) 11, 12 and determining the values of the parameters of the model 10. The parameters are the values that define the s-wave propagation in the model 10. The parameters comprise: s-wave reflection coefficient of the Earth's surface 13, the s-wave reflection coefficient at any (or every) boundary 14, 17 between adjacent layers 11, 12, 16 (these reflection coefficients may be different to or the same as one another, where there are multiple boundaries 14, 17), the s-wave velocity of the different layers 11, 12, 16 (which may be different to or the same as each other), quality factor of the different layers 11, 12, 16 (which may be different to or the same as each other) and/or the density of the different layers 11, 12, 16 (which may be different to or the same as each other).

(51) Finding the depths of the base of the layers 11, 12, and the values of the parameters of the model 10 for the different layers 11, 12, 16 is achieved by inverting the s-wave data.

(52) Step 3 can be achieved by using the s-wave data directly. However, with a view of FIGS. 3 and 4, the seismic data recorded at a buried receiver may be referred to as first seismic data and a second seismic data set may be recorded at the Earth's surface (e.g. the sea bed) near to the location of the receiver that records the first seismic data. When this is the case, step 3 can be achieved by comparing the first and second data sets as follows.

(53) With regard to FIG. 3, FIG. 3a shows exemplary second seismic data recorded at the seabed. FIG. 3b shows exemplary first seismic data recorded beneath the sea bed at a location (in latitude and longitude) close to where the second seismic data is recorded. The depth of the receiver that recorded the first seismic data in FIG. 3b is around 0.8 m. The depth of the receiver that recorded the second seismic data in FIG. 3a is around 0.0 m.

(54) As mentioned above, the p-wave arrivals are excluded from the first and second seismic data by removing or muting the data before a certain arrival times at given offsets (before the dotted line on each data set).

(55) FIG. 3 shows the first and second seismic data in the time domain. Both data sets are converted into the frequency domain, and then the ratio of the amplitude of the first seismic data to the amplitude of the second seismic data is taken. This ratio (i.e. the amplitude of the first seismic data/the amplitude of the second seismic data) in the frequency domain is displayed in FIG. 4.

(56) As can be seen in FIG. 4, there are two large notches in the frequency spectrum of the first seismic data in comparison to the second seismic data. The inventors have found that this is because of s-wave ghosts in the first seismic data.

(57) To find the correct number of layers 11, 12, 16, the depths of the layers 11, 12, 16 and the parameters in the model 10, the model 10 is found that gives a similar (or as near as possible) frequency spectrum to the frequency spectrum shown in FIG. 4. This is achieved by inverting the frequency spectrum of FIG. 4 to find the model 10.

(58) Thus, the model 10 can be found.

(59) It should be noted that there is no need to use the ratio of the first and second seismic data sets to identify notches and for the inversion calculation. Notches may be clear simply from the seismic data recorded at the buried receiver itself, in which case there is no need to use seismic data recorded at the surface to find the model 10. The location of these identified notches in the frequency spectrum, can be used directly in the inversion to determine the parameters in the model.

(60) In order to ensure the model 10 found by inverting the s-wave data is reasonable, the model may be compared to measured geotechnical data. This geotechnical data could also be used as an additional constrain on the inversion calculation when finding the model 10.

(61) Returning to FIG. 2, once the model has been found in step 3, the model can be used in step 4 to compute a deghosting operator. This can be achieved by producing an impulse response using the model.

(62) For one receiver, an example of an impulse response (the light grey line marked “Modeled with one layer ghost”) compared to measured seismic s-wave data (the dark grey line marked “buried receiver/seabed receiver”) can be seen in FIG. 7. This impulse response has been calculated by taking a ratio of the impulse response at the receiver location 15 to the impulse response at the Earth's surface 13 in a one-layer model 10, such as that shown in FIG. 5. This ratio can then be compared to the ratio of the measured s-wave data at the buried receiver to the measured s-wave data at a receiver on the Earth's surface. The comparison can be seen in FIG. 7. As can be seen, the modeled impulse response closely models the main two notches shown in the measured s-wave data. This is an indication that the found model 10 is a reasonable model.

(63) From the impulse response, the inverse of the impulse response can be found. The inverse of the impulse response may be the deghosting operator.

(64) Returning to FIG. 2, once the deghosting operator has been found in step 4, the deghosting operator can be applied to the recorded s-wave data to deghost the s-wave data in step 5. This comprises applying the deghosting operator to the s-wave data by convolving the measured s-wave data with the deghosting operator. With regard to FIG. 7, this effectively means that the inverse function of the light grey line (which is the ratio of the amplitude of the response function at the receiver location 15 to the amplitude of the response function at the Earth's surface 13) is multiplied with (or added in case where a dB scale is used) the dark grey line (which is the ratio of the amplitude of the measured s-wave data at the buried receiver to the amplitude of the measured s-wave data at the Earth's surface). Whilst the inverse of light grey line in FIG. 7 only gives the magnitude of the deghosting operator, the de-ghosting operator also includes a phase term.

(65) Applying such a deghosting operator should remove the notches and hence flatten the frequency spectrum. Further, it should correct the phase, shifting the events to their correct positions.

(66) The output of step 5 is deghosted s-wave data for the s-wave data that was measured at a buried receiver.

(67) As a demonstration of the improved data, FIG. 8a shows one particular buried receiver's raw s-wave data gather in the time domain. The horizontal axis is source offset and the vertical axis is time. As can be seen, the data is quite blurred and unclear. FIG. 8b shows the same buried receiver's s-wave data gather after the present method has been used to deghost the s-wave data. In comparison to FIG. 8a, the data is much sharper, and more details can be discerned, as is highlighted in the black circles.

(68) Steps 1 to 5 can be repeated for different buried receivers, when there is a plurality of buried receivers. A different model 10 may be constructed for each receiver since the environment in the location of each buried receiver may be different.

(69) Once the deghosted seismic s-wave data has been produced using the above steps, it can be used to more accurately image a reservoir.

(70) As a demonstration of the improved imaging, FIG. 9a shows an image of a reservoir produced using raw s-wave data gathered using buried receivers. The horizontal axis is receiver number (which equates to horizontal displacement) and the vertical axis is depth. As can be seen, the image is quite blurred and unclear so some of the details of reservoir cannot be seen, and there are several places where the continuity of subsurface reflectors is interrupted. Some of these issues are highlighted using the circles in FIG. 9a. FIG. 9b shows an image of the same reservoir produced using the same s-wave data as FIG. 9a, but said s-wave data having been deghosted using the present method. As can be seen, in comparison to FIG. 9a, the image is much sharper, more details can be discerned, and the continuity of the subsurface reflectors has been improved, as is highlighted in the black circles.

(71) An example of the application of an embodiment of the invention to permanent reservoir monitoring data will now be presented.

(72) A pure data-dependent method based on flattening of the spectra will not be sufficient to achieve deghosting, since the phase cannot be properly corrected. Instead, the method works under the assumption that the signal distortion, observed as notches, is created by the seabed itself and by strong contrasts in the shear wave velocity between a limited number of layers, typically less than 10 metres below the seabed.

(73) In this example, the following assumptions are made: The shear wave is travelling predominantly vertically in the very shallow so that the simple expressions for normal incidence reflection coefficients can be used. Mode conversions are neglected. Only S-S reflections and transmissions are considered. A horizontally layered model with three layers is considered. This is illustrated in FIG. 10.

(74) FIG. 10 illustrates the horizontally layered model with three layers that is used in this example. Here, the thickness of the first two layers are L.sub.1 and L.sub.2, while the third layer is a half-space. v.sub.1, v.sub.2 and v.sub.3 are the shear velocities. A shear wave u.sub.3 is propagating from below, and the response s.sub.r is measured at the receiver (red triangle). This response is a combination of the direct wave and all intra-layer multiples, but can in principle be described as the sum of up-going and down-going waves. u.sub.1 and d.sub.1 are the total up-going and down-going fields at the top of layer 1, and so forth. The reflection coefficients R.sub.1, R.sub.2 and R.sub.3 are defined by the velocity contrasts at the top of each layer.

(75) For a plane shear wave, the normal incidence reflection and transmission coefficients for the horizontal particle velocity is:

(76) $\begin{array}{cc}{R}_v=\frac{I_1-I_2}{I_1+I_2},T_v=1+R_v,& \left(1\right)\end{array}$

(77) where I.sub.1 is the shear impedance of the medium of the incoming wave, and I.sub.2 is the shear impedance of the medium on the other side of the interface. Neglecting density changes, we use
R.sub.1≈1, R.sub.2=(v.sub.2−v.sub.1)/(v.sub.1+v.sub.2), R.sub.3=(v.sub.3−v.sub.2)/(v.sub.2+v.sub.3)   (2)

(78) Note that R.sub.v=−R. Hence, for up-going and down-going shear waves impinging on the top of a layer i, the reflection coefficient is R.sub.i and −R.sub.i (where i=1, 2 or 3) and the transmission coefficient is 1+R.sub.i. and 1−R.sub.i, respectively. A propagation angle different from vertical may be allowed, as well as attenuation given by a Q-value in each layer. The propagation angle in a layer is given by Snell's law and is determined by the velocity and angle in the bottom layer as input to the calculation. For plane waves propagating with an angle θ to the vertical, the propagation of down-going and up-going waves is (since {right arrow over (k)}.Math.{right arrow over (Δz)}=kΔz cos θ and Δt=Δz/cos θ)

(79) $\begin{array}{cc}d\left(z+\Delta z\right)=d\left(z\right)e^{-i\stackrel{.fwdarw.}{k}.Math.\stackrel{.fwdarw.}{\Delta z}}{e}^{-\frac{\omega \Delta t}{2Q}}=e^{-\mathrm{ik}\Delta z.Math.P\left(\theta ,Q\right)},u\left(z+\Delta z\right)=u\left(z\right)e^{\mathrm{ik}\Delta z.Math.P\left(\theta ,Q\right)}& \left(3\right)\end{array}$

(80) where k=k.sub.1=2πf/v.sub.1 and so forth. The propagation correction factor in each layer is

(81) $\begin{array}{cc}P\left(\theta ,Q\right)=\cos \theta \left(1-\frac{i}{2Q}\left(1+{\tan }^2\theta \right)\right)& \left(4\right)\end{array}$

(82) Hence, the equations connecting the up-going and down-going fields at the top of the three layers are:
u.sub.1=u.sub.2(1+R.sub.2)e.sup.−ik.sup.1.sup.L.sup.1.sup..Math.P(θ.sup.1.sup.,Q.sup.1.sup.)−R.sub.2d.sub.1e.sup.−2ik.sup.1.sup.L.sup.1.sup..Math.P(θ.sup.1.sup.,Q.sup.1.sup.)  (5)
u.sub.2=u.sub.3(1+R.sub.3)e.sup.−ik.sup.2.sup.L.sup.2.sup..Math.P(θ.sup.2.sup.,Q.sup.2.sup.)−R.sub.3d.sub.2e.sup.−ik.sup.2.sup.L.sup.2.sup..Math.P(θ.sup.2.sup.,Q.sup.2.sup.)  (6)
d.sub.1=R.sub.1u.sub.1  (7)
d.sub.2=d.sub.1(1−R.sub.2)e.sup.−ik.sup.1.sup.L.sup.1.sup..Math.P(θ.sup.1.sup.,Q.sup.1.sup.)+R.sub.2u.sub.2  (8)
d.sub.3=d.sub.2(1−R.sub.3)e.sup.−ik.sup.2.sup.L.sup.2.sup..Math.P(θ.sup.2.sup.,Q.sup.2.sup.)+R.sub.3u.sub.3  (9)

(83) The last terms in equation (5) and (6) represent the internal multiples. A solution of equation (5)-(9) in terms of the upcoming wave u.sub.3 can then be written as an explicit scheme starting at the bottom and first moving upwards and then downwards in the layers:
u.sub.2=M.sub.MM.sub.2(1+R.sub.3)e.sup.−ik.sup.2.sup.L.sup.2.sup..Math.P(θ.sup.2.sup.,Q.sup.2.sup.)u.sub.3  (10)
u.sub.1=M.sub.1(1+R.sub.2)e.sup.−ik.sup.1.sup.L.sup.1.sup..Math.P(θ.sup.1.sup.,Q.sup.1.sup.)u.sub.2  (11)
where the internal multiple terms are

(84) $\begin{array}{cc}{M}_1=\frac{1}{1+R_1{R}_2{e}^{-2{\mathrm{ik}}_1{L}_1.Math.P\left({\theta }_1,Q_1\right)}}& \left(12\right)\\ M_2=\frac{1}{1+R_2{R}_3{e}^{-2{\mathrm{ik}}_1{L}_1.Math.P\left({\theta }_1,Q_1\right)}}& \left(13\right)\\ M_M=\frac{1}{\begin{array}{c}1+M_1{M}_2{R}_1{R}_3\left(1-R_1\right)\left(1+R_2\right)\\ e^{-2{\mathrm{ik}}_1{L}_1.Math.P\left({\theta }_1,Q_1\right)}{e}^{-2{\mathrm{ik}}_2{L}_2.Math.P\left({\theta }_2,Q_2\right)}\end{array}}& \left(14\right)\end{array}$

(85) The total field at the receiver depth z.sub.r is the sum of the (propagated) up-going and down-going fields in the respective layer. Finally, normalisation is performed against the up-going field at the receiver in a model with no internal multiples.

(86) The modelling scheme is used in a deghosting strategy for the recorded horizontal components, i.e. radial component after rotation and designature. The P-events (early arrivals) are muted from the data, and one average spectrum is calculated for each receiver. An initial three-layered model is established that gives a fair match (number of notches, notches in approximately the same place) to the receiver spectra. Notches are picked both on synthetic (s.sub.r) and real spectra, and the layered model is then obtained by minimizing an L.sup.2 objective function for the differences in notch position. Weighting factors are used allowing to e.g. put more weight on the first notch, since it is usually less uncertain as well as the most important to correct for. The inversion problem is solved by a local search.

(87) During the inversion, the model is re-parameterized into travel times and reflection coefficients. This allows for more robust inversion with fewer inversion parameters. First, the initial depth model is transformed to travel times and reflection coefficients: T.sub.r (two-way travel time from receiver to surface), T.sub.1 and T.sub.2 (two-way travel times in layers 1 and 2), and reflection coefficients R.sub.1, R.sub.2 and R.sub.3. This transform between depth/velocities and travel times/reflection coefficients does not have to be exact (a zero-incidence assumption is used for this), as long as the forward and inverse transform is exact. Typically, in the travel time/reflection coefficient domain two to three parameters can be inverted for and still get an acceptable (sometimes perfect) match to the picked notches. This parameterization is not used in the direct calculation of the synthetic s in the objective function. The model is instead transformed back to the depth/velocity domain and the synthetic response s.sub.r is calculated with initial angle and Q-factor.

(88) An example of synthetic response at the receiver level and at the seabed (surface) can be seen in FIG. 11.

(89) FIG. 11 shows typical synthetic PS spectra for a shallow layered model, showing the receiver response s.sub.r (solid line), as well as up-going u.sub.1 at the surface (dashed line), and the ratio s.sub.r/u.sub.1 (dash-dot), which is the chosen response for deghosting filter design. FIG. 11(a) shows the dB level and FIG. 11(b) shows the unwrapped phase.

(90) Although the measured receiver burial depth is available, it is inverted for as it has been found that it is not correct in all positions, in particular for very shallow receivers. The output of the inversion process is a shallow earth model consisting of thicknesses and shear wave velocities at each receiver.

(91) The subsequent deghosting can be performed in two different ways.

(92) The first is to inverse filter the recorded data with the synthetic receiver response (s.sub.r). This should yield the upcoming field (.sub.u3) from the subsurface.

(93) However, instead it was chosen to inverse filter the upcoming field with s.sub.r/u.sub.1. This has two advantages. First, it is slightly less model-dependent, as this ratio depends only on the properties above the receiver. Secondly, it brings the data to the surface, allowing the result to be compared with the previous non-buried OBS data recorded at the field. This choice results in the application of a 180-degree phase change at each notch. This is obviously a drastic modification of the data, but, if done correctly, it is what is desired to be achieved.

(94) This deghosting method was applied to permanent reservoir monitoring (PRM) data from the Grane field in the North Sea.

(95) Initially, the data was limited to two notches and a maximum frequency of 50 Hz. A weight factor of 3 for the first notch, and 1 for the second was used, and T.sub.r, T.sub.1 and R.sub.2 were inverted for.

(96) Using these inversion results as an initial model, a final inversion used three notches, a maximum frequency of 42 Hz and inverted for T.sub.r, T.sub.1, T.sub.2, R.sub.2 and R.sub.3 simultaneously.

(97) Q-values were not inverted for since they do not influence the notch positions, only the notch strength. Q=100 was used during the inversion, as a higher Q in this step increased accuracy, and Q=10 was used in the final filters, as this produced synthetic spectra closer to real spectra and stabilized the deghosting filter.

(98) The inversion predicted a shallow (0-2 m) earth model with S-velocities in the uppermost layer in the order of 20-40 m/s. Available geotechnical shallow soil sampling and cone penetration tests (CPT) in the Grane area reports very loose sand and very soft clay close to the seabed that could be consistent with such low S-velocities.

(99) Deghosting filters were designed and applied to radial gathers before demultiple. After a necessary update of statics, deghosted data was subsequently run through a full processing flow, including whitening. The result was compared to data with exactly the same processing but where deghosting had not been applied. It was also compared with images from previous non-buried ocean bottom data. It was found that deghosting clearly improved continuity and resolution in the final image, as well as better phase behaviour at major reflected events. The improved continuity and resolution in the final image is shown in FIG. 12.

(100) FIG. 12 shows sections of final migrated PS full stacks from the same location. The vertical axes are TWT (ms). FIG. 12 (a) shows previous non-buried survey (different processing), FIG. 12(b) shows PRM PS without deghosting, and FIG. 12(c) shows PRM PS with deghosting. Events in FIG. 12(c) are more focused, laterally continuous, and with improved phase in FIG. 12(c).

(101) Thus, using a method according to an embodiment of the invention, the deghosted data has more similar character to data from non-buried receivers, both in the overburden and at the reservoir level. The continuity of reflection events are improved, the holes in the spectrum are filled, the de-ghosted data responds better to whitening and the resolution is improved.

(102) It should be apparent that the foregoing relates only to the preferred embodiments of the present application and the resultant patent. Numerous changes and modification may be made herein by one of ordinary skill in the art without departing from the general spirit and scope of the invention as defined by the following claims and the equivalents thereof.