Methods and systems that determine a velocity wavefield from a measured pressure wavefield

0Associated
Cases 
0Associated
Defendants 
0Accused
Products 
0Forward
Citations 
0
Petitions 
1
Assignment
First Claim
1. In a process for generating an image of a subterranean formation using marine seismic techniques in which one or more sources are activated to generate acoustic energy that is reflected from the subterranean formation and recorded as pressure data by pressure sensors located in one or more streamers, the specific improvement comprising:
 computing an approximate frozen freesurface profile of a free surface above a streamer based on pressure data generated by pressure sensors located at receiver locations along the streamer, the pressure data represents the measured pressure wavefield;
computing reflectivity and normal derivative of the reflectivity at the receiver locations based on the approximate frozen freesurface profile;
computing a normal derivative of the measured pressure wavefield at the receiver locations based on the reflectivity and the normal derivative of the reflectivity at the receiver locations;
computing an approximate normalvelocity wavefield at the receiver locations along the streamer based on the normal derivative of the measured pressure wavefield at the receiver locations, the approximate normalvelocity wavefield does not include lowfrequency, streamer vibrational noise; and
generating an image of the subterranean formation based at least in part on the measured pressure wavefield and the approximate normalvelocity wavefield.
1 Assignment
0 Petitions
Accused Products
Abstract
Methods and systems that compute an approximate verticalvelocity wavefield based on a measured pressure wavefield and knowledge of the freesurface shape when the pressure wavefield was measured are described. The measured pressure wavefield is used to compute an approximate frozen freesurface profile of the free surface. The approximate frozen freesurface profile and the measured pressure wavefield are then used to compute an approximate verticalvelocity wavefield that does not include lowfrequency streamer vibrational noise. The approximate verticalvelocity wavefield and measured pressure wavefield may be used to separate the pressure wavefield into upgoing and downgoing pressure wavefields.
22 Citations
No References
Method for imaging a seasurface reflector from towed dualsensor streamer data  
Patent #
US 7,872,942 B2
Filed 10/14/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Method of wavefield extrapolation for singlestation, dualsensor towed streamer signals  
Patent #
US 7,929,373 B2
Filed 11/19/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

System for combining signals of pressure sensors and particle motion sensors in marine seismic streamers  
Patent #
US 7,684,281 B2
Filed 04/14/2008

Current Assignee
PGS Americas Inc.

Sponsoring Entity
PGS Americas Inc.

Method for attenuation of multiple reflections in seismic data  
Patent #
US 7,675,812 B2
Filed 06/30/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Method for imaging a seasurface reflector from towed dualsensor streamer data  
Patent #
US 20100091610A1
Filed 10/14/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Method of summing dualsensor towed streamer signals using seismic reflection velocities  
Patent #
US 20100027375A1
Filed 08/01/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Method for attenuating particle motion sensor noise in dual sensor towed marine seismic streamers  
Patent #
US 7,835,225 B2
Filed 10/11/2006

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

System for combining signals of pressure sensors and particle motion sensors in marine seismic streamers  
Patent #
US 7,359,283 B2
Filed 03/03/2004

Current Assignee
PGS Americas Inc.

Sponsoring Entity
PGS Americas Inc.

Apparatus and methods for multicomponent marine geophysical data gathering  
Patent #
US 7,239,577 B2
Filed 08/30/2002

Current Assignee
Petroleum GeoServices ASA

Sponsoring Entity
PGS Americas Inc.

Method for imaging of prestack seismic data  
Patent #
US 7,286,690 B2
Filed 10/23/2003

Current Assignee
PGS Americas Inc.

Sponsoring Entity
PGS Americas Inc.

Method for processing dual sensor seismic data to attenuate noise  
Patent #
US 6,894,948 B2
Filed 01/29/2003

Current Assignee
PGS EXPLORATION US LIMITED

Sponsoring Entity
PGS EXPLORATION US LIMITED

Angle dependent surface multiple attenuation for twocomponent marine bottom sensor data  
Patent #
US 6,654,693 B2
Filed 10/22/2001

Current Assignee
PGS Americas Inc.

Sponsoring Entity
PGS Americas Inc.

Method of reverberation removal from seismic data and removal of dual sensor coupling errors  
Patent #
US 5,754,492 A
Filed 02/12/1996

Current Assignee
PGS TENSOR INC.

Sponsoring Entity
PGS TENSOR INC.

Method and device for attenuating water column reverberations using colocated hydrophones and geophones in ocean bottom seismic processing  
Patent #
US 5,774,416 A
Filed 04/07/1995

Current Assignee
PGS TENSOR INC.

Sponsoring Entity
PGS TENSOR INC.

Method of reverberation removal from seismic data and removal of dual sensor coupling errors  
Patent #
US 5,825,716 A
Filed 05/07/1997

Current Assignee
PGS Tensor Inc. Houston TX, PGS Tensor Inc.

Sponsoring Entity
PGS TENSOR INC.

Method and apparatus for improved seismic prospecting  
Patent #
US 5,193,077 A
Filed 10/02/1990

Current Assignee
B P Exploration Company Limited, Atlantic Richfield Company Incorporated

Sponsoring Entity
BP Exploration, Atlantic Richfield Company Incorporated

Method of summing dualsensor towed streamer signals using crossghosting analysis  
Patent #
US 8,089,825 B2
Filed 08/29/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Combined impulsive and nonimpulsive seismic sources  
Patent #
US 8,427,901 B2
Filed 12/21/2009

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Array grouping of seismic sensors in a marine streamer for optimum noise attenuation  
Patent #
US 8,553,490 B2
Filed 11/09/2007

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

METHODS AND SYSTEMS FOR RECONSTRUCTION OF LOW FREQUENCY PARTICLE VELOCITY WAVEFIELDS AND DEGHOSTING OF SEISMIC STREAMER DATA  
Patent #
US 20140016436A1
Filed 07/10/2012

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

Attenuating seismic interference noise using a dual sensor recording system  
Patent #
US 8,811,115 B2
Filed 08/14/2008

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

METHODS AND SYSTEMS FOR ATTENUATING NOISE IN SEISMIC DATA  
Patent #
US 20150063064A1
Filed 12/16/2013

Current Assignee
PGS Geophysical As

Sponsoring Entity
PGS Geophysical As

29 Claims
 1. In a process for generating an image of a subterranean formation using marine seismic techniques in which one or more sources are activated to generate acoustic energy that is reflected from the subterranean formation and recorded as pressure data by pressure sensors located in one or more streamers, the specific improvement comprising:
computing an approximate frozen freesurface profile of a free surface above a streamer based on pressure data generated by pressure sensors located at receiver locations along the streamer, the pressure data represents the measured pressure wavefield; computing reflectivity and normal derivative of the reflectivity at the receiver locations based on the approximate frozen freesurface profile; computing a normal derivative of the measured pressure wavefield at the receiver locations based on the reflectivity and the normal derivative of the reflectivity at the receiver locations; computing an approximate normalvelocity wavefield at the receiver locations along the streamer based on the normal derivative of the measured pressure wavefield at the receiver locations, the approximate normalvelocity wavefield does not include lowfrequency, streamer vibrational noise; and generating an image of the subterranean formation based at least in part on the measured pressure wavefield and the approximate normalvelocity wavefield.  View Dependent Claims (2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13)
 14. A computer system for generating an image of a subterranean formation, the system comprising:
one or more processors; one or more datastorage devices; and machinereadable instructions stored in one or more of datastorage devices and executed by the one or more processors, the machinereadable instructions directed to computing an approximate frozen freesurface profile of a free surface above a streamer based on pressure data generated by pressure sensors located at receiver locations along the streamer, the pressure data represents the measured pressure wavefield; computing reflectivity and normal derivative of the reflectivity at the receiver locations based on the approximate frozen freesurface profile; computing a normal derivative of the measured pressure wavefield at the receiver locations based on the reflectivity and the normal derivative of the reflectivity at the receiver locations; computing an approximate normalvelocity wavefield at the receiver locations along the streamer based on the normal derivative of the measured pressure wavefield at the receiver locations, the approximate normalvelocity wavefield does not include lowfrequency streamer vibrational noise; and generating an image of the subterranean formation based at least in part on the measured pressure wavefield and the approximate normalvelocity wavefield.  View Dependent Claims (15, 16, 17, 18, 19, 20, 21, 28)
 22. A nontransitory computerreadable medium having machinereadable instructions encoded thereon for generating an image of a subterranean formation and enabling one or more processors of a computer system to perform the operations of
computing an approximate frozen freesurface profile of a free surface above a streamer based on pressure data generated by pressure sensors located at receiver locations along the streamer, the pressure data represents the measured pressure wavefield; computing reflectivity and normal derivative of the reflectivity at the receiver locations based on the approximate frozen freesurface profile; computing a normal derivative of the measured pressure wavefield at the receiver locations based on the reflectivity and the normal derivative of the reflectivity at the receiver locations; computing an approximate normalvelocity wavefield at the receiver locations along the streamer based on the normal derivative of the measured pressure wavefield at the receiver locations, the approximate normalvelocity wavefield does not include lowfrequency, streamer vibrational noise; and generating an image of the subterranean formation based at least in part on the measured pressure wavefield and the approximate normalvelocity wavefield.  View Dependent Claims (23, 24, 25, 26, 27, 29)
1 Specification
This application claims the benefit of Provisional Application 62/043915, filed Aug. 29, 2014.
In recent years, the petroleum industry has invested heavily in the development of improved marine survey techniques and seismic data processing methods in order to increase the resolution and accuracy of seismic images of subterranean formations. Marine surveys illuminate a subterranean formation located beneath a body of water with acoustic signals produced by one or more submerged seismic sources. The acoustic signals travel down through the water and into the subterranean formation. At interfaces between different types of rock or sediment of the subterranean formation a portion of the acoustic signal energy may be refracted, a portion may be transmitted, and a portion may be reflected back toward the formation surface and into the body of water. A typical marine survey is carried out with a survey vessel that passes over the illuminated subterranean formation while towing elongated cablelike structures called streamers. The streamers may be equipped with a number of collocated, dual pressure and particle motion sensors that detect pressure and particle motion wavefields, respectively, associated with the acoustic signals reflected back into the water from the subterranean formation. The pressure sensors generate pressure data that represents the pressure wavefield and the particle motion sensors generate particle motion data that represents the particle motion wavefield. The survey vessel receives and records the seismic data generated by the sensors.
A wavefield that propagates upward from the subterranean formation and is detected by the pressure or particle motion sensors is called an upgoing wavefield. Ideally, a pressure or particle motion upgoing wavefield alone may be used to compute a seismic image of the subterranean formation. However, the surface of the water acts as a nearly perfect acoustic reflector. As a result, the sensors also detect a downgoing wavefield created by reflection of the upgoing wavefield from the water surface. The downgoing wavefield is essentially the upgoing wavefield with a time delay that corresponds to the amount of time it takes for acoustic signals to travel up past the streamers to the water surface and back down to the streamers. The downgoing wavefield creates “ghost” effects in seismic images of the subterranean formation. Typical seismic data processing techniques use both the pressure wavefield and particle motion wavefield to capture the upgoing pressure and particle motion wavefields. However, the particle motion sensor may also record lowfrequency noise, such as lowfrequency noise created by streamer vibrations, that contaminates seismic images produced from the upgoing wavefields. Those working in the field of marine exploration seismology continue to seek methods and systems that improve seismic image quality.
Methods and systems that compute an approximate velocity wavefield based on a measured pressure wavefield and knowledge of the freesurface profile when the pressure wavefield was measured are described. The measured pressure wavefield is composed of pressure data generated by pressure sensors of one or more streamers. The measured pressure wavefield is used to compute an approximate frozen freesurface profile of the free surface. A frozen freesurface profile represents the frozenintime profile or shape of the free surface above the streamers at the time the pressure wavefield was measured. The approximate frozen freesurface profile and the measured pressure wavefield may be used to compute an approximate verticalvelocity wavefield which does not include lowfrequency noise typically measured by particle motion sensor. The approximate verticalvelocity wavefield and measured pressure wavefield may also be used to separate the pressure wavefield into upgoing and downgoing pressure wavefields. The upgoing pressure wavefield may, in turn, be used to compute seismic images with significantly higher resolution and deeper signal penetration than seismic images produced from seismic data contaminated with the downgoing wavefield and lowfrequency noise. High resolution seismic images may improve the accuracy of identifying hydrocarbon reservoirs and monitoring production of hydrocarbon reservoirs. In other words, typical wavefield separation is performed with measured pressure and verticalvelocity wavefields measured by collocated pressure and particle motion sensors. But methods and systems described herein eliminate the need for particle motion sensors and a measured verticalvelocity wavefield, because an approximate verticalvelocity wavefield may be obtained from the measured pressure wavefield alone, and wavefield separation may be performed with the measured pressure wavefield and the approximate verticalvelocity wavefield.
Streamer depth below the free surface 112 can be estimated at various locations along the streamers using depthmeasuring devices attached to the streamers. For example, the depthmeasuring devices can measure hydrostatic pressure or utilize acoustic distance measurements. The depthmeasuring devices can be integrated with depth controllers, such as paravanes or water kites that control and maintain the depth and position of the streamers (or portions thereof) as the streamers are towed through the body of water. The depthmeasuring devices are typically placed at intervals (e.g., about 300 meter intervals in some implementations) along each streamer. Note that in other implementations buoys may be attached to the streamers and used to maintain the orientation and depth of the streamers below the free surface 112.
The secondary waves may be generally emitted at different times within a range of times following the initial acoustic signal. A point on the formation surface 122, such as the point 138, may receive a pressure disturbance from the primary wavefield more quickly than a point within the subterranean formation 120, such as points 140 and 142. Similarly, a point on the formation surface 122 directly beneath the source 104 may receive the pressure disturbance sooner than a more distantlying point on the formation surface 122. Thus, the times at which secondary and higherorder waves are emitted from various points within the subterranean formation 120 may be related to the distance, in threedimensional space, of the points from the activated source.
Acoustic and elastic waves, however, may travel at different velocities within different materials as well as within the same material under different pressures. Therefore, the travel times of the primary wavefield and secondary wavefield emitted in response to the primary wavefield may be functions of distance from the source 104 as well as the materials and physical characteristics of the materials through which the wavefields travel. In addition, the secondary expanding wavefronts may be altered as the wavefronts cross interfaces and as the velocity of sound varies in the media are traversed by the wave. The superposition of waves emitted from within the subterranean formation 120 in response to the primary wavefield may be a generally complicated wavefield that includes information about the shapes, sizes, and material characteristics of the subterranean formation 120, including information about the shapes, sizes, and locations of the various reflecting features within the subterranean formation 120 of interest to exploration seismologists.
Each receiver 118 may be a multicomponent sensor composed of one or more particle motion sensors that detect particle motion, velocities, or accelerations over time and a pressure sensor that detects variations in water pressure over time.
The particle motion sensors are typically oriented so that the particle motion is measured in the vertical direction (i.e., {right arrow over (n)}=(0, 0, z)) in which case v_{z}({right arrow over (x)}_{r}, t) is called the verticalvelocity data and a_{z}({right arrow over (x)}_{r}, t) is called the vertical acceleration data. Alternatively, each receiver may include two additional particle motion sensors that measure particle motion in two other directions, {right arrow over (n)}_{1 }and {right arrow over (n)}_{2}, that are orthogonal to {right arrow over (n)}(i.e., {right arrow over (n)}·{right arrow over (n)}_{1}={right arrow over (n)}·{right arrow over (n)}_{2}=0, where “·” is the scalar product) and orthogonal to one another (i.e., {right arrow over (n)}_{1}·{right arrow over (n)}_{2}=0). In other words, each receiver may include three particle motion sensors that measure particle motion in three orthogonal directions. For example, in addition to having a particle motion sensor that measures particle velocity in the zdirection to give v_{z}({right arrow over (x)}_{r}, t), each receiver may include a particle motion sensor that measures the wavefield in the inline direction in order to obtain the inline velocity data, v_{x}({right arrow over (x)}_{r}, t), and a particle motion sensor that measures the wavefield in the crossline direction in order to obtain the crossline velocity data, v_{y}({right arrow over (x)}_{r}, t). In certain implementations, the particle motion sensors may be omitted and the receivers may be composed of only pressure sensors.
The streamers 106111 and the survey vessel 102 may include sensing electronics and dataprocessing facilities that allow seismic data generated by each receiver to be correlated with the time the source 104 is activated, absolute positions on the free surface 112, and absolute threedimensional positions with respect to an arbitrary threedimensional coordinate system. The pressure data and particle motion data may be stored at the receiver, and/or may be sent along the streamers and data transmission cables to the survey vessel 102, where the data may be stored electronically or magnetically on datastorage devices located onboard the survey vessel 102. The pressure data and particle motion, velocity, or acceleration data represent pressure and particle motion wavefields, and therefore, may also be referred to as the measured pressure wavefield and measured particle motion, velocity, or acceleration wavefield, respectively.
Returning to
As explained above, each pressure sensor 204 and particle motion sensor 206 generates seismic data that may be stored in datastorage devices located at the receiver or onboard the survey vessel. Each pressure sensor and particle motion sensor may include an analogtodigital converter that converts timedependent analog signals into discrete time series that consist of a number of consecutively measured values called “amplitudes” separated in time by a sample rate. The time series generated by a pressure or particle motion sensor is called a “trace,” which may consist of thousands of samples collected at a typical sample rate of about 1 to 5 ms. A trace is a recording of a subterranean formation response to acoustic energy that passes from an activated source into the subterranean formation where a portion of the acoustic energy is reflected and/or refracted and ultimately detected by a sensor as described above. A trace records variations in a timedependent amplitude that represents acoustic energy in the portion of the secondary wavefield measured by the sensor. The coordinate location of each time sample generated by a moving sensor may be calculated from global position information obtained from one or more global positioning devices located along the streamers, survey vessel, and buoys and the known geometry and arrangement of the streamers and sensors. A trace generated by a pressure sensor or particle motion sensor is wavefield data that may be represented as a set of timedependent amplitudes denoted by:
tr_{r}(t)={a_{r}(t_{j})}_{j=1}^{J} (1)
 where
 j is a time sample index;
 J is the number of time samples; and
 a_{r}(t_{l}) is the pressure or particle motion amplitude at time sample
For example, p({right arrow over (x)}_{r}, t) is the trace generated by a pressure sensor and v_{{right arrow over (n)}}({right arrow over (x)}_{r}, t) is the trace generated by a particle motion sensor. Each trace may also include a trace header not represented in Equation (1) that identifies the specific receiver that generated the trace, receiver GPS coordinates, and may include time sample rate and the number of samples.
 where
As explained above, the secondary wavefield typically arrives first at the receivers located closest to the sources. The distance from the sources to a receiver is called the “sourcereceiver offset,” or simply “offset,” which creates a delay in the arrival time of a secondary wavefield from an interface within the subterranean formation. A larger offset generally results in a longer arrival time delay. The traces are collected to form a “gather” that can be further processed using various seismic data processing techniques in order to obtain information about the structure of the subterranean formation.
The arrival times versus sourcereceiver offset is longer with increasing sourcereceiver offset. As a result, the wavelets generated by a formation surface or an interface are collectively called a “reflected wave” that tracks a hyperbolic curve. For example, hyperbolic curve 420 represents the hyperbolic distribution of the wavelets 410414 reflected from the formation surface 122, which are called a “formationsurface reflected wave,” and hyperbolic curve 422 represents the hyperbolic distribution of the wavelets 415419 from the interface 124, which are called an “interface reflected wave.”
A gather is a collection of traces that represents the pressure or verticalvelocity wavefield measured by corresponding pressure sensors or particle motion sensors. The gather shown in
In practice, however, measured pressure and verticalvelocity wavefields do not share the same broad frequency spectrum. For example, pressure sensors typically have a high signaltonoise ratio over a broad frequency range, but particle motion sensors often detect lowfrequency, streamer vibrational motion that contaminates the lowfrequency part of the verticalvelocity data. As a result, particle motion sensors may have a low signaltonoise ratio over the lowfrequency part of the frequency range.
Lowfrequency noise (e.g., vibrational noise) may be observed in the frequency spectra of the verticalvelocity wavefield. The pressure data represented by each trace in the measured pressure wavefield may be transformed from the spacetime (“st”) domain to the spacefrequency (“sf”) domain using a Fast Fourier Transform (“FFT”) or a discrete Fourier transform (“DFT”):

 where ω is the angular frequency.
Likewise, the verticalvelocity data represented by each trace in the verticalvelocity wavefield may be transformed from the st domain using an FFT or a DFT to the sf domain:
 where ω is the angular frequency.
Methods and systems described in greater detail below compute an approximate verticalvelocity wavefield from the measured pressure wavefield. In particular, the pressure data p({right arrow over (x)}_{r}, t) of each trace in the pressure wavefield is transformed from the st domain to wavenumberfrequency (“kf”) domain as follows:

 where
 k_{x }is the xdirection or inline wavenumber; and
 k_{y }is the ydirection or crossline wavenumber.
Methods and systems compute approximate verticalvelocity wavefield from the measured pressure wavefield in the kf domain:
P(k_{x}, k_{y}, ωz_{r})→{tilde over (V)}_{z}(k_{x}, k_{y}, ωz_{r}) (5)
 where “˜” represents approximate.
Approximate verticalvelocity data may replace the verticalvelocity data over the lowfrequency range 711 provided (1) the pressure data has a satisfactory signaltonoise ratio over the lowfrequency range, (2) the pressure spectrum of the pressure data has no notches over the lowfrequency range, and (3) the depth of the pressure sensors are known. As shown inFIG. 7 , the relative amplitude of the pressure spectrum 706 exhibits notches 712, 714 and 716 that depend on the depth of the streamer. The notches 712, 714, and 716 are shifted toward lower frequencies as streamer depth increases and shifted toward higher frequencies as streamer depth decreases. For the example spectra shown inFIG. 7 , the pressure spectrum 706 does not have notches in the lowfrequency range 711, indicating that the pressure data over the lowfrequency range may be used to calculate an estimated verticalvelocity wavefield over the lowfrequency range 711.
 where
In one implementation, the lowfrequency part of the verticalvelocity data may be replaced by the approximate verticalvelocity data in the kf domain using a combined verticalvelocity wavefield given by:

 where
 ω_{n }is an upper limit on the lowfrequency range;
 ω_{c }is the cutoff frequency; and
 W_{L }is a lowpass filter and W_{H }is a highpass filter that satisfy a condition
W_{L}+W_{H}=1 (7)
The cutoff frequency ω_{c }is less than the lowest frequency notch in the pressure spectrum. For example, with reference toFIG. 7 , the cutoff frequency should be less than about 60 Hz, which corresponds to the first notch 712. For ω_{n}<ω≤ω_{c}, the low and highpass filters may be frequency dependent:
 where
In another implementation in which recording of the verticalvelocity wavefield is omitted or the verticalvelocity wavefield is contaminated with noise over large portions of the frequency range, the approximate verticalvelocity wavefield may be computed over the full frequency range of the pressure data and used to substantially remove receiver ghost effects as follows. In the kf domain, the pressure data may be represented as a sum of upgoing pressure data and downgoing pressure data as follows:
P(k_{x}, k_{y}, ωz_{r})=P^{up}(k_{x}, k_{y}, ωz_{r})+P^{down}(k_{x}, k_{y}, ωz_{r}) (9)
 where
 P^{up}(k_{x}, k_{y}, ωz_{r}) represents the upgoing pressure data in the kf domain; and
 P^{down}(k_{x}, k_{y}, ωz_{r}) represents the downgoing pressure data in the kf domain.
The pressure data and approximate verticalvelocity wavefield may be used to separate the pressure data into upgoing and downgoing pressure data in the kf domain as follows:
 where

 where
 ρ is the density of water; and
 where
is the zdirection or vertical wavenumber with c_{0 }the speed of sound in water.
The separate upgoing and downgoing pressure data may be transformed from the kf domain back to the st domain using an inverse FFT (“IFFT”), or inverse (“IDFT”), to obtain corresponding separate upgoing and downgoing pressure data in the st domain. The upgoing wavefield may then be used to compute seismic images of the subterranean formation substantially free of the receiver ghost effects contained in the downgoing pressure data.
Methods and systems compute an approximate velocity wavefield based on the measured pressure data and acoustic reflectivity of the free surface above the streamers when the pressure data is measured. In particular, methods and systems first compute the normal derivative of the pressure data at receiver locations along a streamer based on the measured pressure data and acoustic reflectivity of the free surface. The normal derivative may then be used to compute the normal derivative of the particle velocity at the same receiver locations along the streamer which, in turn, may be used to compute the verticalvelocity of the wavefield at the receiver locations along the streamer. Normal derivatives of the pressure data at receiver locations along the streamer are computed from a Green'"'"'s second identity formulation of a sourcefree, closed surface that relates the pressure data measured at the receivers located along the streamer to the normal derivative of the measured pressure data and includes a Green'"'"'s function representation of acoustic reflectivity from a spatiotemporally varying free surface above the streamer at the time the pressure data is measured.
In order to relate the pressure wavefield to normal derivatives of the pressure wavefield at receiver locations along the streamer in the actual state shown in FIG. 9A, the two states represented in
∫_{V}d{right arrow over (x)}[P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω)∇^{2}G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω)−G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs},ω)∇^{2}P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω)]=∫_{S}dS{right arrow over (n)}·[P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω)∇G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω)−G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω)∇P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω)] (11)
 where
 P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω) represents the pressure data generated by a receiver at coordinate location {right arrow over (x)}_{r }resulting from a source at coordinate location {right arrow over (x)}_{s};
 G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω) represents a Green'"'"'s function at coordinate location {right arrow over (x)}_{r }and virtual source coordinate location {right arrow over (x)}_{vs}; and
 V represents the volume of the space enclosed by the surface S.
The pressure data, P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω), satisfies the Helmholtz wave equation given by
 where

 where
 A(ω) results from a Fourier transformation of a sourcetime function a(t); and
 c^{2}(χ)=c_{0}^{2}(1−α(χ)) with c_{0 }the speed of sound in water and α(χ) the refractive index.
The Green'"'"'s function, G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω), characterizes reflections from the free surface S_{fs }and is a solution of the acoustic wave equation for a Dirac delta pulse represented by a Dirac delta function as follows:
 where
Utilizing the Sommerfeld radiation condition (i.e., R→∞) and substituting the Helmholtz wave Equation (12a) and the acoustic wave Equation (12b) into Equation (11), Equation (11) reduces to a surface integral equation over the surface of the streamer S_{r }as follows:

 where
represents the normal derivative directed orthogonal to the streamer S_{r}; and
The unit normal vector {right arrow over (n)} in Equation (13) is given by

 where
and
 z(x, y) represents streamer depth.
The streamer depth z(x, y) at points only along the streamer S_{r }may be interpolated based on the depths generated by depthmeasuring devices located along the streamer.
 z(x, y) represents streamer depth.
In order to relate the measured pressure data P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω) to the normal derivative of the pressure data {right arrow over (n)}·∇P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω), Equation (13) is solved for the case α=½. In other words, the coordinates {right arrow over (x)}_{vs }of the virtual sources (i.e. the Green'"'"'s function sources) are moved to the streamer to obtain the following integral equation:
The surface integral given by Equation (14) numerically as follows. The surface S_{r }is broken into K small area elements dxdy. The pressure data and its normal derivative are considered constant over each of the small area elements and are equal to their respective values at the center of each element. Mapping the area elements onto a flat surface using
dS_{r}=γdxdy
and moving the virtual sources to coincide with actual receivers along the streamer S_{r}, Equation (14) may be rewritten as the following system of equations:

 where
 M_{vs,r}=G({right arrow over (x)}_{r}, {right arrow over (x)}_{vs}, ω)γdxdy;
 where


 “·” represents the scalar or dot product;
 {right arrow over (x)}_{vs }is a vsth virtual source coordinate along the streamer S_{r }(i.e. both the virtual sources ({right arrow over (x)}_{vs}) and the actual receivers ({right arrow over (x)}_{r}) are at the same location); and
 vs=1, . . . , K.
The normal derivative of the pressure wavefield at each of the K receiver locations may be determined from the set of Equations (15) by first rewriting the set of Equations (15) in matrix form as follows:


 where

M is called a “monopole matrix” with matrix elements M_{vs,r}; 
D is called a “dipole matrix” with matrix elements D_{vs,r}; and  {right arrow over (P)}=[P({right arrow over (x)}_{1}, {right arrow over (x)}_{s}, ω) P({right arrow over (x)}_{2}, {right arrow over (x)}_{s}, ω) . . . P({right arrow over (x)}_{K}, {right arrow over (x)}_{s}, ω)]^{T}.
Note that because the virtual sources coincide with the receivers located along the streamer S_{r}, P({right arrow over (x)}_{r}, {right arrow over (x)}_{s}, ω) equals P({right arrow over (x)}_{vs}, ω) for the index r equal to the index vs in Equation (16) and the source coordinates may be suppressed. The diagonal elements the monopole and dipole matricesM andD are singular and may be replaced by estimates over the discretization path as follows:

 where
The normal derivative of the measured pressure wavefield at each {right arrow over (x)}_{r}, {right arrow over (n)}·∇P({right arrow over (x)}_{r}, ω), may be computed as a function of the measured pressure data at each {right arrow over (x)}_{r}, P({right arrow over (x)}_{r}, ω), by solving for ∂{right arrow over (P)}/∂n in Equation (16) to obtain:

 where Ī is a K×K identity matrix.
Equation (18) may be used to compute the normal derivative of the measured pressure wavefield at each receiver location {right arrow over (x)}_{r}, {right arrow over (n)}·∇P({right arrow over (x)}_{r}, ω). Approximate normalvelocity data at each receiver location {right arrow over (x)}_{r }along the streamer may be calculated according to:
 where Ī is a K×K identity matrix.
The approximate normalvelocity data computed for a number of receiver locations forms an approximate normalvelocity wavefield.
The approximate normalvelocity data {tilde over (V)}_{{right arrow over (n)}}({right arrow over (x)}_{r}, ω) approximates the normal particle velocity that a particle motion sensor would have measured at the receiver location {right arrow over (x)}_{r}. However, unlike the normal particle velocity data that an actual particle motion sensor would generate, the approximate normalvelocity data {tilde over (V)}_{{right arrow over (n)}}({right arrow over (x)}_{r}, ω) does not include the lowfrequency noise (e.g., does not include streamer vibration noise). Approximate verticalvelocity data {tilde over (V)}_{z}({right arrow over (x)}_{r}, ω) may be computed at the receiver location {right arrow over (x)}_{r }from the approximate normalvelocity data {tilde over (V)}_{{right arrow over (n)}}({right arrow over (x)}_{r}, ω) based on the orientation of the receiver as shown in
The resulting approximate verticalvelocity data at the receiver location {right arrow over (x)}_{r }is given by:
{tilde over (V)}_{z}({right arrow over (x)}_{r}, ω)=cos ϕ·{tilde over (V)}_{{right arrow over (n)}}({right arrow over (x)}_{r}, ω) (21)
The approximate verticalvelocity data {tilde over (V)}_{z}({right arrow over (x)}_{r}, ω) also does not include the lowfrequency noise typically recorded by a particle motion sensor. The approximate verticalvelocity data {tilde over (V)}_{z}({right arrow over (x)}_{r}, ω) may be transformed from the spacefrequency domain to the kf domain using an FFT or DFT to obtain {tilde over (V)}_{z}(k_{x}, k_{y}, ωz_{r}), which may be used to compute combined verticalvelocity wavefield given by Equation (6) and separate approximate upgoing and downgoing pressure wavefields according to Equations (10a) and (10b). As a result, the combined verticalvelocity wavefield and the separate approximate upgoing and downgoing pressure wavefields do not include the lowfrequency noise.
Calculation of the Green'"'"'s function and normal derivative of the Green'"'"'s function used in Equation (18) to compute the normal derivative of the pressure wavefield at the receiver location {right arrow over (x)}_{r}, {right arrow over (n)}·∇P({right arrow over (x)}_{r}, ω), are now described. The Green'"'"'s function is the reflectivity of the free surface S_{fs }and the normal derivative of the Green'"'"'s function is the normal derivative of the reflectivity at the free surface. The following description presents a technique for computing the Green'"'"'s function and normal derivative of the Green'"'"'s function used in Equation (18).
A Green'"'"'s function for the model shown in
G({right arrow over (x)}′, {right arrow over (x)}_{vs}, ω)_{{right arrow over (x)}′={right arrow over (x)}}_{fs}=0
in order to obtain

 where
 G is the Green'"'"'s function or reflectivity of the free surface;
 where
represents the normal derivative directed orthogonal to the free surface S_{fs};
and
 G^{0 }is a free space Green'"'"'s function.
The freespace Green'"'"'s functions appearing in Equation (22a) and in subsequent Equations below is represented in general by:
 G^{0 }is a free space Green'"'"'s function.

 where
 i is the imaginary unit √{square root over (−1)};
 “∥·∥” is Euclidean distance;
 k_{0}=ω/c_{0}; and
 {right arrow over (y)} and {right arrow over (y)}′ represent general coordinate locations of two different points in the space represented in
FIGS. 11A11B .
The unit normal vector {right arrow over (n)}′ at any point along the free surface in Equation (22a) is given by
 where

 where
and
 f(x, y) approximates the profile or shape of the free surface S_{fs }above the streamer (i.e., f(x, y) is the approximate frozen freesurface profile).
Methods for computing the approximate frozen freesurface profile f(x, y) of the free surface S_{fs }are described below with reference toFIGS. 1221 . Note that the Green'"'"'s function G({right arrow over (x)}′, {right arrow over (x)}_{vs}, ω) represents the reflectivity resulting from two sources. The first source is the Dirac delta pulse and the second source is the free surface. A technique for computing the Green'"'"'s function and the normal derivative of the Green'"'"'s function at streamer locations {right arrow over (x)}_{r }based on Equation (22a) is summarized as follows. First, the parameter β in Equation (22a) is set to one. Next, letting {right arrow over (x)}′_{vs }approach the free surface and using the boundary condition G({right arrow over (x)}′_{vs}, ω)_{{right arrow over (x)}′}_{vs}_{={right arrow over (x)}}_{fs}=0, Equation (22a) becomes:
 f(x, y) approximates the profile or shape of the free surface S_{fs }above the streamer (i.e., f(x, y) is the approximate frozen freesurface profile).
Next, the computed normal derivative of the Green'"'"'s function at the free surface given by Equation (23a) is replaced in Equation (22a) and then setting β=1 and {right arrow over (x)}_{fs}={right arrow over (x)}_{r}, the Green'"'"'s function at the streamer location {right arrow over (x)}_{r }is given by:
Finally, taking the normal derivative of the Green'"'"'s function given by Equation (23b) at the receiver locations gives
The Green'"'"'s function G({right arrow over (x)}_{r}, ω) in Equation (23b) represents the reflectivity at receiver locations along the streamer and {right arrow over (n)}·∇G({right arrow over (x)}_{r}, ω) in Equation (23c) represents the normal derivative of the reflectivity at receiver locations along the streamer.
The Green'"'"'s function G({right arrow over (x)}_{r}, ω) and the normal derivative if {right arrow over (n)}·∇G({right arrow over (x)}_{r}, ω) may be computed numerically and substituted into Equation (18). The surface S_{fs }is broken into K small area elements dxdy. The Green'"'"'s function and the normal derivative of the Green'"'"'s function are considered constant over each small area element and are equal to their respective values at the center of each small area element. By mapping curved small area elements onto a flat surface
dS_{fs}=δdxdy
Equation (23a) may be rewritten as
Equation (23b) may be rewritten as
and Equation (23c) may be rewritten as
Equations (24a)(24c) may be used to compute the Green'"'"'s function G, which is then substituted into Equation (18) to compute the normal derivative of the pressure wavefield.
The shape or profile of the free surface at the time the pressure wavefield was measured is assumed to be in a fixed or frozenintime state called a “frozen free surface.” The approximate frozen freesurface profile f(x, y) approximates the actual, unknown, frozenintime state or profile of the free surface above the pressure sensors when the pressure wavefield was measured, which is used to compute the unit normal vector {right arrow over (n)}′. Computing an approximate freesurface profile above a streamer is now described with reference to Equations (25)(38) and
The fixed crosssectional shape of an approximate frozenfree surface is computed by forward and backward extrapolation of the pressure wavefield at a series of trial depths above the streamer.
Extrapolation is carried out by first transforming the pressure data generated by each receiver from the st domain to the kf domain as follows:
The transformation may be executed with an FFT or a DFT. Next, at each trial depth, the pressure data generated at each receiver is forward and backward extrapolated to the trial depth level to obtain forward and backward extrapolated wavefields that correspond to the trial depth. For each receiver, the pressure data is forward extrapolated to a trial depth Z_{n }according to
P^{F}(k_{x}, k_{y}, ωZ_{n})=P(k_{x}, k_{y}, ωz_{r})e^{−ik}^{z}^{(z}^{r}^{−z}^{n}^{)} (26)
and backward extrapolated to the same trial depth Z_{n }according to
P^{B}(k_{x}, k_{y}, ωZ_{n})=P(k_{x}, k_{y}, ωz_{r})e^{ik}^{z}^{(z}^{r}^{−z}^{n}^{)} (27)
 where
 k_{x }is the horizontal wavenumber in the xdirection;
 k_{y }is the horizontal wavenumber in the ydirection; and
 k_{z}=√{square root over (k^{2}−k_{x}^{2}−k_{y}^{2})}.
For each trial depth Z_{n}, the forward and backward extrapolated pressure data associated with each receiver is then transformed from the kf domain back to the st domain:
 where

 where
 superscript “F” represents forward extrapolated; and
 superscript “B” represents backward extrapolated.
The transformation may be executed with an IFFT or an IDFT.
 where
For each trial depth Z_{n}, the forward extrapolated pressure data computed for each receiver are collected to form a forward extrapolated gather
p^{F}(x, y, tZ_{n})={p^{F}(x_{r}, y_{r}, tZ_{n})}_{r=1}^{R} (29a)
 where R is the number of receivers.
Backward extrapolated pressure data computed for each receiver are also collected to form a backward extrapolated gather
p^{B}(x, y, tZ_{n})={p^{B}(x_{r}, y_{r}, tZ_{n})}_{r=1}^{R} (29b)
 where R is the number of receivers.
A difference gather may be computed for each pair of forward and backward extrapolated gathers. In other words, for each trial depth Z_{n}, a difference gather is computed as follows:
d(x, y, tZ_{n})=p^{F}(x, y, tZ_{n})−p^{B}(x, y, tZ _{z}) (30)
Amplitudes of the forward extrapolated gather p^{F}(x, y, tZ_{n}) are represented by a_{r}^{F}(t_{j}Z_{n}) and amplitudes of the backward extrapolated gather p^{B}(x, y, tZ_{n}) are represented by a_{r}^{B}(t_{j}Z_{n}), where r is a trace index, t_{j }is the time sample index, and Z_{n }identifies the trial depth. Amplitude differences, A_{r}(t_{j}Z_{n}), for each difference gather d(x, y, tZ_{n}) may be executed according to the following pseudo code:
A series of time windows are applied to each difference gather and a maximum energy difference is calculated for each time window.
E(x_{r}, y_{r}, t_{j}∈W_{m}Z_{n})=[A_{r}(t_{j}Z_{n})]^{2} (31)
Rectangle 1708 represents the time window W_{m }with energy differences calculated according to Equation (31) for each of the amplitude differences in the time window W_{m}. For example, solid circle 1710 represents the energy difference E(x_{9}, y_{9}, t_{j}Z_{n}) for the 9th trace at the jth time sample computed according to Equation (31). A maximum energy difference is determined for each trace with time samples in the time window W_{m }according to:
The quantity E_{max}(x_{r}, y_{r}, W_{m}Z_{n}) represents the maximum energy difference for the rth trace in the time window W_{m}. In
{right arrow over (E)}_{max}(W_{m}Z_{n})=[E_{max}(x_{1}, y_{1}, W_{m}Z_{n}) . . . E_{max}(x_{R}, y_{R}, W_{m}Z_{n})]^{T} (33)
For example, in
Next, for each time window applied to the N difference gathers, a peak energy is computed from the N vectors of maximum energy differences calculated for the time window.

 where z_{peak,r }equals the trial depth Z_{n }of the maximum energy difference in the set {E_{max}(x_{r}, y_{r}, W_{m}Z_{n})}.
Each maximum energy difference E_{max}(x_{r}, y_{r}, W_{m}Z_{n}) in the set {E_{max}(x_{r}, y_{r}, W_{m}Z_{n})} corresponds to a receiver as described above with reference toFIG. 17B . As a result, receiver coordinates (x_{r}, y_{r}) associated with the peak energy E_{peak}(x_{r}, y_{r}, z_{peak,r}) are the receiver coordinates associated with z_{peak,r}, which includes the subscript r to identify the receiver. The peaks z_{peak,r }and associated receiver coordinates (x_{r}, y_{r}) are collected to form a set of points {(x_{r}, y_{r}, Z_{peak,r})} that approximate the shape of the frozen free surface above the streamer.
 where z_{peak,r }equals the trial depth Z_{n }of the maximum energy difference in the set {E_{max}(x_{r}, y_{r}, W_{m}Z_{n})}.
The sideelevation shape or profile of the frozen free surface above the streamer may be approximated by interpolating over the set of points {(x_{r}, y_{r}, z_{peak,r})}. For example, polynomial interpolation, spline interpolation, and Gaussian interpolation may be used to compute a function or polynomial that approximates the sideelevation shape or profile of the frozen free surface above the streamer based on the set of points {(x_{r}, y_{r}, Z_{peak,r})}.

 where for the integer index q≥0,

 and for q<0, F(K_{q})=F(K_{−q})*.
The parameter W(K_{q}) is the PiersonMoskowitz spatial roughness spectrum, which for a fully developed sea surface in onedimension (e.g., xdirection) is given by:
 and for q<0, F(K_{q})=F(K_{−q})*.

 where
 K_{q }is the spatial wavenumber;
 U_{w }is the wind speed measured at a height of about 19 meters;
 α is 8.0×10^{−3};
 β is 0.74; and
 g is the acceleration due to gravity.
In Equations (36) and (37), the spatial wavenumber for component q is given by K_{q}=2πq/L, where L is the length of free surface. The random number N(0,1) may be generated from a Gaussian distribution having zero mean and a unit variance. As a result, the free surface is formed by adding each wavenumber component imposing random phase shifts. A frozenintime PiersonMoskowitz free surface may be computed from Equation (36) using a FFT for computational efficiency.
 where
The frozen freesurface extension f_{ext}(χ) may be combined with the approximate frozen freesurface profile f_{int}(xχ) to form an approximate frozen freesurface profile given by
In alternative implementations, the approximate frozen freesurface extension may be expanded to include a time parameter that characterizes the frozen free surface at different times. Freesurface waves are generally dispersive and in deep water, and the frequency and wavenumber are related by a dispersion relation given by:
Ω(K_{q})=√{square root over (gK_{q})} (39)
Equation (39) implies that each spatial harmonic component of the freesurface wavefield may move with a definite phase velocity. As a result, in general, freesurface waves of longer wavelengths travel faster relative to waves with shorter wavelengths. Combining Equations (36) and (39) gives a timevarying frozen free surface:

 where t is instantaneous time.
Equation (40) characterizes a onedimensional rough free surface moving in the positive xdirection and may be used to compute the frozen freesurface extension 2102 at earlier or later times.
 where t is instantaneous time.
Consider a freesurface shape at an instant in time t with wave heights given by Equation (40). The wavenumber spectrum F(K_{q}) of the free surface is computed according to Equation (36) and an arbitrary known dispersion relation Ω(K_{q}) is calculated according to Equation (39) may be used to calculate the frozen free surface at an earlier (t−Δt) or a later (t+Δt) time by:
The mathematical equations and formulas presented above are not, in any way, intended to mean or suggest an abstract idea or concept. Instead the mathematical equations described above provide a concise conceptual representation computational operations that may performed seismic data. The mathematical equations and methods described above are ultimately implemented on physical computer hardware, datastorage devices, and communications systems in order to obtain results that also represent physical properties, materials, and deposits within the earth'"'"'s interior. For example, as explained above, a pressure wavefield emanating from an actual subterranean formation after being illuminated with an acoustic signal is composed of actual physical pressure waves that are sampled using physical and concrete pressure sensors. The pressure sensors in turn produce physical electrical or optical signals that encode pressure data that is recorded on physical datastorage devices and undergoes computational processing using hardware as describe above to compute an approximate verticalvelocity wavefield that may in turn be used to compute upgoing wavefield data that represents physical and concrete upgoing pressure and verticalvelocity wavefields. The approximate verticalvelocity wavefield or the upgoing wavefield data may be displayed, or subjected to further seismic data processing, in order to interpret the physical structure and composition of the subterranean formation, such as in monitoring production of an actual hydrocarbon deposit within the subterranean formation.
Although the above disclosure has been described in terms of particular implementations, it is not intended that the disclosure be limited to these implementations. Modifications within the spirit of this disclosure will be apparent to those skilled in the art. For example, any of a variety of different implementations may be obtained by varying any of many different design and development parameters, including programming language, underlying operating system, modular organization, control structures, data structures, and other such design and development parameters.
The method described above may be implemented in real time while a marine survey is being conducted or subsequent to completion of the marine survey. The measured pressure wavefield and approximate normalvelocity wavefield may form geophysical data products indicative of certain properties of a subterranean formation. The geophysical data products may include processed seismic data and may be stored on a computerreadable medium as described above. The geophysical data products may be produced offshore (i.e. by equipment on the survey vessel 102) or onshore (i.e. at a computing facility on land) either within the United States or in another country. When the geophysical data product is produced offshore or in another country, it may be imported onshore to a datastorage facility in the United States. Once onshore in the United States, geophysical analysis may be performed on the data product.
It is appreciated that the previous description of the disclosed embodiments is provided to enable any person skilled in the art to make or use the present disclosure. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the generic principles defined herein may be applied to other embodiments without departing from the spirit or scope of the disclosure. Thus, the present disclosure is not intended to be limited to the embodiments shown herein but is to be accorded the widest scope consistent with the principles and novel features disclosed herein.