Methods and systems to separate wavefields using pressure wavefield data

This disclosure is directed to wavefield separation methods and systems. In one aspect, methods and systems compute an approximate vertical particle velocity wavefield based on a measured pressure wavefield and knowledge of free-surface when the pressure wavefield was measured. The measured pressure wavefield is used to compute an approximate frozen free-surface profile. The approximate frozen free-surface profile and the measured pressure wavefield are used to compute an approximate vertical particle velocity wavefield. The approximate vertical particle velocity wavefield and measured pressure wavefield may be used to compute separate up-going and down-going pressure, or vertical particle velocity, wavefields.

BACKGROUND

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 cable-like structures called streamers. The streamers may be equipped with a number of collated, dual pressure and particle motion sensors that detect pressure and vertical particle motion wavefields, respectively, associated with the acoustic signals reflected back into the water from the subterranean formation. The pressure sensors generate seismic data that represents the pressure wavefield and the particle motion sensors generate seismic data that represents the vertical particle motion wavefield. The survey vessel receives and records the seismic data generated by the sensors.

A wavefield that travels upward from the subterranean formation and is detected by the pressure or particle motion sensors is called an up-going wavefield, which 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 down-going wavefield created by reflection of the up-going wavefield from the water surface. The down-going wavefield is essentially the up-going 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 down-going wavefield combines with the up-going wavefield, resulting in recorded seismic data contaminated with, unwanted down-going wavefield energy that creates “ghost” effects in seismic images of the subterranean formation computed from the seismic data. Typical seismic data processing techniques use both the pressure wavefield and vertical particle motion wavefield to separate the pressure and vertical particle motion wavefields into up-going and down-going wavefields. The up-going wavefield may be used to compute an image of a subterranean formation without the ghost effects caused by the down-going wavefield.

DETAILED DESCRIPTION

This disclosure is directed to methods and systems that perform wavefield separation. The methods and systems compute an approximate vertical particle velocity wavefield based on a measured pressure wavefield and knowledge of a free-surface when the pressure wavefield was measured. The measured pressure wavefield is used to compute an approximate frozen free-surface profile of the free-surface. The approximate frozen free-surface profile and the measured pressure wavefield are used to compute an approximate vertical particle velocity wavefield. The approximate vertical particle velocity wavefield and measured pressure wavefield may be used to compute separate up-going and down-going pressure wavefields or vertical particle velocity wavefields. The up-going pressure wavefield or vertical particle velocity wavefield may, in turn, be used to compute seismic images with significantly higher resolution and deeper signal penetration than seismic images computed with seismic data contaminated with the down-going wavefield. Removal of the effects of the down-going wavefield may reduce unwanted noise during marine surveying or reservoir production monitoring.

FIGS. 1A-1Bshow side-elevation and top views, respectively, of an example seismic data acquisition system composed of a survey vessel102towing a source104and six separate streamers106-111beneath a free-surface112of a body of water. The body of water may be, for example, an ocean, a sea, a lake, or a river, or any portion thereof. In this example, each streamer is attached at one end to the survey vessel102via a streamer-data-transmission cable. The illustrated streamers106-111form a planar horizontal data acquisition surface with respect to the free-surface112. However, in practice, the data acquisition surface may be smoothly varying, for example, due to active sea currents and weather conditions. In other words, although the streamers106-111are illustrated inFIGS. 1A and 1Band subsequent figures as straight and substantially parallel to the free-surface112, in practice, the towed streamers may undulate as a result of dynamic conditions of the body of water in which the streamers are submerged. A data acquisition surface is not limited to having a planar horizontal orientation with respect to the free-surface112. The streamers may be towed at depths that angle the data acquisition surface with respect to the free-surface112or one or more of the streamers may be towed at different depths. A data acquisition surface is not limited to six streamers as shown inFIG. 1B. In practice, the number of streamers used to form a data acquisition surface can range from as few as one streamer to as many as 20 or more streamers. It should also be noted that the number of sources is not limited to a single source. In practice, the number of sources selected to generate acoustic energy may range from as few as one source to three or more sources and the sources may be towed in groups by one or more survey vessels.

FIG. 1Aincludes an xz-plane114andFIG. 1Bincludes an by-plane116of the same Cartesian coordinate system having three orthogonal, spatial coordinate axes labeled x, y and z. The coordinate system is used to specify orientations and coordinate locations within the body of water. The x-direction specifies the position of a point in a direction parallel to the length of the streamers (or a specified portion thereof when the length of the streamers are curved) and is referred to as the “in-line” direction. The y-direction specifies the position of a point in a direction perpendicular to the x-axis and substantially parallel to the free-surface112and is referred to as the “cross-line” direction. The z-direction specifies the position of a point perpendicular to the by-plane (i.e., perpendicular to the free-surface112) with the positive z-direction pointing downward away from the free-surface112. The streamers106-111are generally long cables containing power and data-transmission lines that connect receivers represented by shaded rectangles118spaced-apart along the length of each streamer to seismic data acquisition system and data-storage devices located on board the survey vessel102.

Streamer depth below the free-surface112may be estimated at various locations along the streamers using depth measuring devices attached to the streamers. For example, the depth measuring devices can measure hydrostatic pressure or utilize acoustic distance measurements. The depth measuring devices may be integrated with depth controllers, such as para vanes or water kites that control and maintain the depth and position of the streamers as the streamers are towed through the body of water. The depth measuring 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 help maintain the orientation and depth of the streamers below the free-surface112.

FIG. 1Ashows a cross-sectional view of the survey vessel102towing the source104above a subterranean formation120. Curve122represents a top surface of the subterranean formation120located at the bottom of the body of water. The subterranean formation120is composed of a number of subterranean layers of sediment and rock. Curves124,126, and128represent interfaces between subterranean layers of different compositions. A shaded region130, bounded at the top by a curve132and at the bottom by a curve134, represents a subterranean hydrocarbon deposit, the depth and positionally coordinates of which may be determined, at least in part, by analysis of seismic data collected during a marine survey. As the survey vessel102moves over the subterranean formation120, the source104is activated to produce an acoustic signal (often referred to as a “shot”) at spatial and/or temporal intervals. In other implementations, the source104may be towed by one survey vessel and the streamers may be towed by a different survey vessel. The source104may be an air gun, marine vibrator, or composed of an array of air guns and/or marine vibrators.FIG. 1Aillustrates an acoustic signal expanding outward from the source104as a pressure wavefield136represented by semicircles of increasing radius centered at the source104. The outwardly expanding wavefronts from the sources may be three-dimensional (e.g., spherical) but are shown in vertical plane cross section inFIG. 1A. The outward and downward expanding portion of the pressure wavefield136is called the “primary wavefield,” which eventually reaches the formation surface122, at which point the primary wavefield is partially reflected from the formation surface122and partially refracted downward into the subterranean formation120, becoming elastic waves within the subterranean formation120. In other words, in the body of water, the acoustic signal is composed of compression al pressure waves, or P-waves, while in the subterranean formation120, the waves include both P-waves and transverse waves, or S-waves. Within the subterranean formation120, at interfaces between different types of materials or at discontinuities in density or in one or more of various other physical characteristics or parameters, downward propagating waves may be partially reflected and partially refracted. As a result, each point of the formation surface122and each point of the interfaces124,126, and128may be considered a reflector that becomes a potential secondary point source from which acoustic and elastic wave energy, respectively, may emanate upward toward the receivers118in response to the acoustic signal generated by the source104. As shown inFIG. 1A, secondary waves of significant amplitude may be generally emitted from points on or close to the formation surface122, such as point138, and from points on or very close to interfaces in the subterranean formation120, such as points140and142. The upward expanding secondary waves emitted from the subterranean formation120are collectively called the “secondary wavefield.”

The secondary waves that compose the secondary wavefield may be generally emitted at different times within a range of times following the initial acoustic signal. A point on the formation surface122, such as the point138, may receive a pressure disturbance from the primary wavefield more quickly than a point within the subterranean formation120, such as points140and142. Similarly, a point on the formation surface122directly beneath the source104may receive the pressure disturbance sooner than a more distant-lying point on the formation surface122. Thus, the times at which secondary and higher-order waves are emitted from various points within the subterranean formation120may be related to the distance, in three-dimensional 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 source104as 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 waves. The superposition of waves emitted from within the subterranean formation120in response to the primary wavefield may be a generally complicated wavefield that includes information about the shapes, sizes, and material characteristics of the subterranean formation120, including information about the shapes, sizes, and locations of the various reflecting features within the subterranean formation120of interest to exploration geophysicists.

Each receiver118may comprise a pressure sensor that detects variations in water pressure over time.FIG. 2shows a side-elevation view of the seismic data acquisition system with a magnified view202of the receiver118. In this example, the magnified view202reveals that the receiver118comprises a pressure sensor204. The pressure sensors204may be, for example, hydrophones. Each pressure sensor204may measure non-directional, hydrostatic pressure changes over time and may produces pressure wavefield data denoted by p({right arrow over (x)}r,t), where {right arrow over (x)}r=(xr,yr,zr) represent the receiver Cartesian coordinates, and t represents time. The depth zrof each receiver118may be estimated from the depth measurements obtained from the depth measuring devices located along the streamers.

Seismic data includes data generated by the receivers118when detecting hydrostatic pressure changes over time. The streamers106-111, receivers118, and the survey vessel102may include sensing electronics and data-processing facilities that allow seismic data generated by each receiver118to be correlated with the time and location of each source activation, absolute positions on the free-surface112, and absolute three-dimensional positions with respect to an arbitrary three-dimensional coordinate system. The seismic data may be stored at the receivers118and/or may be sent along the streamers and data transmission cables to the survey vessel102, where the data may be stored electronically or magnetically on data-storage devices located onboard the survey vessel102. The seismic data generated by the receivers118may represent pressure changes in the secondary wavefield emitted from the subterranean formation120.

At least a portion of the secondary wavefield emitted from the subterranean formation120may propagate upward toward the free-surface112, forming an up-going wavefield. InFIG. 2, directional arrow206represents the direction of an up-going wavefield at the location of receiver118, and dashed arrows210and212represents a down-going wavefield produced by the up-going wavefield reflection from the free-surface112before reaching the receiver118. In other words, the pressure wavefield data p({right arrow over (x)}r,t) measured by the receivers118is composed of an up-going pressure wavefield component and a down-going pressure wavefield component. The down-going wavefield contaminates seismic data and creates notches in the seismic data spectral domain.

Each receiver118may generate an electrical or optical signal that encodes the seismic data, which may be recorded on a physical data-storage device that may be located at the receiver, at components along the streamer, or onboard the survey vessel. The seismic data is generally a time series of consecutively measured values, called amplitudes, separated in time by a sample rate. The time-series seismic data measured by a receiver is called a “trace,” which may consist of thousands of time samples of amplitudes collected at a sample rate of about 1 to 5 ms. A trace is generally a record of a subterranean formation response to acoustic energy that passes from an activated source, into the subterranean formation where the reflected acoustic energy is detected by a receiver as described above. A trace generated by a pressure sensor is pressure wavefield data that may be represented as a set of time-dependent pressure amplitudes denoted by:
p({right arrow over (x)}r,t)={ar(tj)}j=1J(1)where r is a positive integer trace, receiver, or channel index;j is a time sample index;J is the number of time samples; andar(tj) is the pressure amplitude of the r-th trace at time sample tj.
Each trace also includes 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.

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 “source-receiver 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. Various sets of traces are collected to form seismic data structures called “gathers” that may be further processed using various seismic data processing techniques in order to extract information about the structure of the subterranean formation.

FIG. 3shows example ray paths that represent an acoustic signal300that travels from the source104into the subterranean formation120. Dashed-line rays, such as rays302, represent acoustic energy reflected from the formation surface122to the receivers118located along the streamer108, and solid-line rays, such as rays304, represent acoustic energy reflected from an interface124to the receivers118located along the streamer108. Note that for simplicity of illustration only a handful of ray paths are represented, and ray paths that extend to deeper interfaces are not shown. Each receiver118measures pressure changes resulting from the acoustic energy reflected from the subterranean formation120. The pressure wavefield data generated at each receiver118, p({right arrow over (x)}r,t), where the receiver subscript r equals 1, 2, 3, 4, and 5, are time recorded as separate traces in one or more data-storage devices as described above with reference to Equation (1). In the example ofFIG. 3, the collection of traces generated by the five receivers118located along the streamer108for a single activation of the source104may be collected to form a seismic data structure called a “common-shot gather.” The traces generated by the receivers118located along each of the six streamers108, shown inFIG. 1B, for the same activation of the source104may be collected to form six separate common-shot gathers.

FIG. 4shows a plot of a synthetic common-shot gather composed of example traces401-405of the pressure wavefield data recorded by the five receives located along the streamer108shown inFIG. 3. Vertical axis430represents time and horizontal axis432represents trace numbers. The traces are arranged so that trace401represents the seismic data generated by the receiver118located closest to the source104and trace405represents the seismic data generated by the receiver118located farthest (along the length of the streamer) from the source104. In this example, the traces401-405represent variation in the amplitude of the seismic data recorded by the five receivers118. The example traces include wavelets or pulses406-410and411-415represented by peaks colored black that represent the up-going wavefield measured by the pressure sensors204. The time intervals along the traces401-405from the trace number axis432(i.e., time zero) to the wavelets406-410represents two-way travel time of the acoustic energy output from the source104to the formation surface122and to the receivers118located along the streamer108, and the wavelets411-415represent longer two-way travel time of the acoustic energy output from the source104to the interface124and to the same receivers118located along the streamer108. The amplitudes of the peaks or troughs of the wavelets406-410and411-415indicate the magnitude of the reflected acoustic energy recorded by the receivers118.

The arrival times of acoustic energy at the receivers increases with increasing source-receiver offset. The wavelets generated by a formation surface and an interface reflection of acoustic energy are collectively called a “reflected wave” or simply “reflection” that, in this example, tracks a parabolic-shaped curve. For example, curve416represents the distribution of the wavelets406-410reflected from the formation surface122, which are called a “formation surface reflected wave”, and curve418represents the distribution of the wavelets411-415from the interface124, which are called an “interface reflected wave” or “interface reflection.”

FIG. 5shows an example expanded view of a synthetic gather composed of 38 traces. Each trace, such as trace502, varies in amplitude over time and represents acoustic energy reflected from the formation surface and five different interfaces within a subterranean formation as measured by a pressure sensor. In the expanded view, wavelets that correspond to reflection from the same formation surface or interface of the subterranean formation appear chained together or to overlap. For example, wavelets504with the shortest transit time represent a formation surface reflection, and wavelets506represent an interface reflected wave emanating from an interface just below the formation surface. Reflected waves508-511represent reflections from interfaces located deeper within the subterranean formation.

The synthetic gathers shown inFIGS. 4 and 5represent ideal cases in which acoustic reflections from features of a subterranean formation are measured directly by the receivers. But, in practice, seismic data typically collected in actual marine surveys record other types of acoustic energy reflections that contaminate the seismic data. For example, seismic data obtained from a marine survey records the up-going wavefield scattered directly from a subterranean formation and the down-going wavefield reflected from the free-surface described above with reference toFIG. 2. The down-going wavefield is essentially a time-delayed up-going wavefield. The down-going wavefield interferes with the up-going wavefield by cancelling frequencies (i.e., creating notches) in the frequency spectrum of the seismic data and creates “ghost” effects in seismic images generated from the seismic data. The up-going wavefield may be separated from the down-going wavefield when the vertical particle velocity wavefield of the acoustic energy reflected from features of the subterranean formation is known. Methods and systems now described compute an approximate vertical particle velocity wavefield based on a measured pressure wavefield and the shape of the free-surface at the time the pressure wavefield is measured. The approximate vertical particle velocity wavefield and the measured pressure wavefield may then be used to compute separate up-going and down-going wavefields.

FIGS. 6A-6Cprovide an overview of how a measured pressure wavefield may be separated into up-going and down-going pressure wavefields without a measured vertical particle velocity wavefield. InFIG. 6A, synthetic common-shot gather600represents a portion of a pressure wavefield measured by a number of receivers, such as receivers located along a streamer, after activation of a source. For the sake of simplicity, the gather600shows only up-going and down-going pressure wavefield components of the pressure wavefield. The solid curves represent the up-going pressure wavefield and dashed curves represent the down-going pressure wavefield. For example, solid curve606represents pressure variations created by a formation surface reflection, and dashed curve608represents pressure variations created by the same formation surface reflection with a time delay610resulting from the time it takes for acoustic energy to travel up past the streamer to the free-surface and back down to the streamer, as described above with reference toFIG. 2. The methods and systems described herein compute602an approximate vertical particle velocity wavefield represented by a gather604based on the pressure wavefield data recorded in the gather600and knowledge of the free-surface shape above the streamers at the time the pressure wavefield is measured. In other words, for each trace in the gather600that corresponds to the pressure wavefield measured at a receiver, an approximate vertical particle velocity trace is computed for the same receiver:
p({right arrow over (x)}r,t)→{tilde over (v)}z({right arrow over (x)}r,t)  (2)

where {tilde over (v)}z({right arrow over (x)}r,t) represents a vertical particle velocity trace at the receiver (xr,yr,zr).

The approximate vertical particle velocity {tilde over (v)}z({right arrow over (x)}r,t) is an approximation of the vertical particle velocity vz({right arrow over (x)}r,t) that may have been obtained from a particle motion sensor collated with the receiver used to measure the pressure wavefield data p({right arrow over (x)}r,t). The approximate vertical particle velocity wavefield represented by the gather604is similarly composed of up-going and down-going vertical particle velocity wavefield components identified by solid and dashed curves, respectively.

Once the approximate vertical particle velocity wavefield is computed, the measured pressure wavefield and the approximate vertical particle velocity wavefield may be used to compute separate up-going and down-going pressure wavefield components of the measured pressure wavefield. InFIG. 6B, the pressure wavefield, denoted by p, is transformed612from the space-time (“s-t”) domain using a fast Fourier transform (“FFT”), or a discrete Fourier transform (“DFT”), to obtain pressure wavefield P in the frequency-wavenumber (“f-k”) domain. In the frequency-wavenumber domain, the pressure wafefield P may be decomposed into a sum of the up-going pressure wavefield and the down-going pressure wavefield as follows:
P=Pup=Pdown(3)

wherePuprepresents the up-going pressure wavefield in the f-k domain; andPdownrepresents the down-going pressure wavefield in the f-k domain.
The approximate vertical particle velocity wavefield, denoted by {tilde over (v)}z, is also transformed614from the s-t domain using an FFT, or a DFT, to obtain an approximate vertical particle velocity wavefield {tilde over (V)}zin the f-k domain. The up-going and down-going pressure wavefields in the f-k domain are computed616according to

whereρ is the density of water;ω is angular frequency; andkzis the z-direction or vertical wavenumber.
In other words, once the approximate vertical particle velocity wavefield {tilde over (V)}zis computed, the pressure wavefield P may be separated, or decomposed, into up-going and down-going pressure wavefields according to Equations (4a) and (4b). The separate up-going and down-going pressure wavefields, Pupand Pdown, may be transformed618and620from the f-k domain back to the s-t domain using an inverse FFT (“IFFT”), or inverse (“IDFT”), to obtain separate up-going and down-going pressure wavefields, pupand pdown, in the s-t domain. Alternatively, the measured pressure wavefield and the approximate vertical particle velocity wavefield may be used to compute approximate up-going and down-going vertical particle velocity wavefield components of the vertical particle velocity wavefield.

InFIG. 6C, the up-going and down-going pressure wavefields combined in the pressure wavefield represented by the gather600ofFIG. 6Aare shown in separate up-going pressure wavefield gather622and down-going pressure wavefield gather624. The pressure wavefield represented by the up-going pressure wavefield in the gather622may be subjected to further seismic data processing to remove noise and other effects and serve as input to imaging methods that generate seismic images of the subterranean formation. The resulting seismic images may have significantly higher resolution and deeper signal penetration into the subterranean formation than seismic images computed with the unseparated seismic data represented in the gather600.

Methods and systems compute an approximate vertical particle velocity wavefield based on the measured pressure wavefield and on knowledge of the shape of the free-surface above the data acquisition surface when the pressure wavefield is measured. Methods and systems described herein compute the shape of the free-surface above the data acquisition surface based on the measured pressure wavefield.FIGS. 7-15illustrate calculation of the shape of a free-surface above a streamer based on a pressure wavefield measured by receivers located in the streamer. The shape of the free-surface at the time the pressure wavefield is measured is assumed to be in a fixed or frozen-in-time state called a frozen free-surface.

FIG. 7shows a side-elevation view of a streamer702located in a body of water below a free-surface704. InFIG. 7, horizontal axis706represents the in-line, or x-direction, vertical axis708represents depth, and circle710represents the cross-line or, y-direction, orthogonal to the xz-plane. The depth zrand the elevation of the free-surface are estimated with respect to the geoid, which corresponds to the x-axis706. The geoid is the hypothetical surface of the earth that coincides with the mean sea level and is used to define zero elevation (i.e., z=0). The streamer702includes 14 spaced apart receivers, such as receiver712, that each measure a different portion of a pressure wavefield. The streamer702may be part of larger data acquisition surface composed of any number of streamers towed by a survey vessel not shown. Each receiver may generate pressure wavefield data p({right arrow over (x)}r,t) as described with reference to Equation (1), such as pressure wavefield data p({right arrow over (x)}s,t) generated for the fifth receiver712of the streamer702.

FIGS. 7-9 and 16included the frozen free-surface704which represents a snapshot of the free-surface above the streamer702when the pressure wavefield data p({right arrow over (x)}r,t) is recorded. The frozen free-surface704has a fixed cross-sectional shape above the streamer702when the pressure wavefield data p({right arrow over (x)}r,t) is generated for each activation of a source, consequently referred to as a “frozen free-surface.” In practice, the actual shape of the frozen free-surface704is not known but the frozen free-surface704is included in the follow illustrations in order to depict computation of an approximate frozen free-surface as now described with reference toFIGS. 8-16.

The fixed cross-sectional shape of the frozen free-surface704is computed by forward and backward extrapolation of the pressure wavefield at a series of trial depths above the streamer.FIG. 8shows the streamer702and frozen free-surface704′ and an example series of trial depths Z1through ZN, where N is a positive integer. Dots, such as dot802, identify the series of regularly spaced trial depths above the receiver712. In this example, the trial depths extend beyond the frozen free-surface704as represented by final trial depth804ZN.

Extrapolation is carried out by first transforming the pressure wavefield data generated by each receiver from the s-t domain to the f-k domain as follows:
p({right arrow over (x)}r,t)→P(kx,ky,ω|zr)  (5)
The transformation may be executed with an FFT or a DFT. Next, at each trial depth, the pressure wavefield 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 wavefield data is forward extrapolated to a trial depth Znaccording to
PF(kx,ky,ω|Zn)=P(kx,ky,ω|zr)e−ikz(zr−zn)(6)
and backward extrapolated to the same trial depth Znaccording to
PB(kx,ky,ω|Zn)=P(kx,ky,ω|zr)eikz(zr−zn)(7)

wherei is the imaginary unit √{square root over (−1)};kxis the horizontal wavenumber in the x-direction;kyis the horizontal wavenumber in the y-direction; andkz=√{square root over (k2−kx2−ky2)}.
For each trial depth Zn, the forward and backward extrapolated pressure wavefield data associated with each receiver is then transformed from the f-k domain back to the s-t domain:
PF(kx,kx,ω|Zn)→pF(xr,yr,t|Zn)  (8a)
PB(kx,ky,ω|Zn)→pB(xr,yr,t|Zn)  (8b)

wheresuperscript “F” represents forward extrapolated; andsuperscript “B” represents backward extrapolated.
The transformation may be executed with an IFFT or an IDFT.

FIG. 9shows an example of forward and backward extrapolated pressure wavefield data at the trial depth Zn902for the fifth receiver712. The forward and backward extrapolated pressure wavefield data at the trial depth Znare represented by pF(x5,y5,t|Zn) and pB(x5,y5,t|Z7), respectively. For each trial depth Znwith n=1, . . . , N, forward and backward extrapolated pressure wavefield data pF(xr,yr,t|Zn) and pB(xr,yr,t|Zn) are computed for each receiver (i.e., r=1, . . . , 14) of the streamer702.

For each trial depth Zn, the forward extrapolated pressure wavefield data computed for each receiver are collected to form a forward extrapolated gather
pF(x,y,t|Zn)={pF(xr,yr,t|Zn)}r=1R(9a)

where R is the number of receivers.

Backward extrapolated pressure wavefield data computed for each receiver are also collected to form a backward extrapolated gather
pB(x,y,t|Zn)={pB(xr,yr,t|Zn)}r=1R(9b)

FIG. 10shows a series of forward extrapolated gathers and corresponding backward extrapolated gathers computed for each trial depth Znwith n=1, . . . , N. Rectangle1002represents a common-shot gather of pressure wavefield data p generated by a number of receivers located along a streamer. Pair of rectangles represent the forward and backward extrapolated gathers computed for each of the trial depths Zn. For example, pair of rectangles1004and1006represent forward and backward extrapolated gathers for the first trial depth Z1and pair of rectangles1008and1010represent forward and backward extrapolated gathers for the final trial depth ZN.

A difference gather may be computed for each pair of forward and backward extrapolated gathers. In other words, for each trial depth Zn, a difference gather is computed as follows:
d(x,y,t|Zn)=pF(x,y,t|Zn)−pB(x,y,t|Zz)  (10)
Amplitudes of the forward extrapolated gather pF(x,y,t|Zn) are represented by arF(tj|Zn) and amplitudes of the backward extrapolated gather pB(x,y,t|Zn) are represented by arB(tj|Zn), where r is a trace index, tjis the time sample index, and Znidentifies the trial depth. Amplitude differences, Ar(tj|Zn), for each difference gather d(x,y,t|Zn) may be executed according to the following pseudo code:

FIG. 11shows difference gathers computed from the pairs of forward and backward extrapolated gathers shown inFIG. 10. Rectangle1102represents a difference gather obtained by computing the difference between the pair of forward and backward extrapolated gathers1004and1006shown inFIG. 10. Rectangle1104represents a difference gather obtained by computing the difference between the pair of forward and backward extrapolated gathers1008and1010shown inFIG. 10.

A series of time windows are applied to each difference gather and a maximum energy difference is calculated for each time window.FIGS. 12A-12Cshow an example calculation of a maximum energy difference for each time window.FIG. 12Ashows an example of a series of time windows applied to a time region1200of a difference gather d(x,y,t|Zn)1202. In this example, the time windows are represented by rectangles that span a time subinterval and the full set of traces of the difference gather1202. The time windows are denoted by Wm, where m=1, . . . , M is the time window index and M is the total number of time windows in the time window series. Initially, time window W1is located over the earliest time interval of the difference gather1202, time window Wmis located over an intermediate time interval, and time window WMis located over a later time interval. In certain implementations the time windows may intersect while in other implementations the time windows may not intersect. In alternative implementations, the time windows and the region of the difference gather the series of time windows are applied to may be hyperbolic in order to approximate the curved shaped of reflections created by source receiver offset in common-shot gathers described above with reference toFIG. 4. In still other implementations, the time series windows may be applied to the full difference gather.

FIG. 12Bshows a rectangle1204that represents an enlargement of the time window Wmlocated over a time interval of the difference gather1202shown inFIG. 12. Solid circles located within the time window Wmrepresent amplitude differences Ar(tj|Zn) calculated according to Equation (10). For example, solid circle1206represents an amplitude difference A9(tj|Zn) calculated for a trace r=9. For each amplitude difference in the time window Wm, a corresponding energy difference is calculated according to
E(xr,yr,tjεWm|Zn)=[Ar(tj|Zn)]2(11)
Rectangle1208represents the time window Wmwith energy differences calculated according to Equation (11) for each of the amplitude differences in the time window Wm. For example, E(x9,y9,tj|Zn) represents the energy difference for the 9-th trace at the j-th time sample determined by computing the square of the amplitude difference A9(tj|Zn) represented by circle1206in rectangle1204. A maximum energy difference is determined for each trace with time samples in the rectangle1208according to:

FIG. 12Cshows the time windows shown inFIG. 12Awith corresponding vectors of maximum energy differences computed for each time window. For example, time window W1has associated vector of maximum energy differences, {right arrow over (E)}max(W1|Zn), and time window W2has associated vector of maximum energy differences, {right arrow over (E)}max(W2|Zn).

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.FIG. 13shows the M time windows applied to the N difference gathers, each of which is identified by a trial depth Zn. For example, time window Wmis applied to each of the N difference gathers over the same time interval and a vector of maximum energy difference {right arrow over (E)}max(Wm|Zn) is computed for each of the N difference gathers. The maximum energy differences of the vectors of maximum energy differences formed for the time window Wmapplied to each of the N difference gathers are collected to form a set of maximum energy differences {Emax(xr,yr,Wm|Zn)} for r=1, . . . , R and n=1, . . . , N. For each time window Wm, a peak energy is then identified from the set of maximum energy differences

where zpeak,requals the trial depth Znof the maximum energy difference in the set {Emax(xr,yr,Wm|Zn)}.

Each maximum energy difference Emax(xr,yr,Wm|Zn) in the set {Emax(xr,yr,Wm|Zn)} corresponds to a receiver as described above with reference toFIG. 12B. As a result, receiver coordinates (xr,yr) associated with the peak energy Epeak(xr,yr,zpeak,r) are the receiver coordinates associated with zpeak,r, which includes the subscript r to identify the receiver. The peaks zpeak,rand associated receiver coordinates (xr,yr) are collected to form a set of points {(xr,yr,zpeak,r)} that approximate the shape of the frozen free-surface above the streamer.

FIG. 14shows a side-elevation view of the streamer702and frozen free-surface704described above with reference toFIG. 7. Open circles located along the frozen free-surface704′ above the receivers represent points in the set {(xr,yr,zpeak,r)}. For example, open circle1402represent a point (x5,y5,zpeak,5) along the frozen free-surface704′ above the fifth receiver712with receiver coordinates (x5,y5,z5).

Implementations are not limited to generating a difference gather as described above with reference toFIG. 11. In alternative implementations, two separate time windows may be applied to the same time intervals of corresponding forward and backward extrapolated gathers and calculation of amplitude differences is limited to the amplitudes within the two time windows.

FIG. 15shows an example calculation of amplitude differences for two separate time windows1502and1504applied to the same time intervals of forward and backward extrapolated gathers1506and1508. In this example, the forward and backward extrapolated gathers1506and1508represent the extrapolated gathers used to compute the difference gather1202described above with reference toFIGS. 12A-12C. Instead of computing the entire difference gather1202, a time window1510containing amplitude differences, arF(tj|Zn)−arB(tj|Zn, of corresponding amplitudes in the two time windows1502and1504is computed for every time sample in the time window1510. The amplitude differences in the time window1510correspond to the amplitude differences in the time window1204inFIG. 12B. The energy differences are computed for corresponding amplitude differences in the time windows1502and1504as follows:
E(xr,yr,tjεWm|Zn)=[arF(tj|Zn)−arB(tj|Zn)]2(14)
The energy differences computed from the amplitude differences in the time window1510correspond to the energy differences in the time window1208shown inFIG. 12Band a maximum energy difference is determined for each trace as described above with reference to Equation (12b). Next, for each time window applied to the N corresponding forward and backward extrapolated gathers, a peak energy is computed as described above with reference to Equation (13) to obtain a set of points {(xr,yr,zpeak,r)} that approximate the shape of the frozen free-surface above the streamer.

The shape of the frozen free-surface above the streamer may be approximated by applying interpolation to the set of points {(xr,yr,zpeak,r)}. For example, polynomial interpolation, spline interpolation, and Gaussian interpolation may be used to compute an approximate frozen free-surface profile above the streamer based on the set of points {(xr,yr,zpeak,r)}.

FIG. 16Ashows a side-elevation view of the streamer702and frozen free-surface704and open circles that represent the set of points [(xr,yr,zpeak,r)]. Dashed curve1602represents an approximate frozen free-surface profile, fint(x), of the frozen free-surface704above the streamer702between the first receiver coordinate x1and the last receiver coordinate xlidentified by dashed lines1606and1608, respectively. Note that the cross-line receiver coordinate yris suppressed, because the approximate frozen free-surface profile is determined in the in-line direction.FIG. 16Aalso shows a source1604(represented by a shaded circle) towed by a survey vessel (not shown). But the approximate frozen free-surface profile1602does not approximate the frozen free-surface above the source1604. The approximate frozen free-surface profile may be extended to approximate the frozen free-surface above the source1604by computing a frozen free-surface extension.

FIG. 16Bshows a plot of an extended portion of the approximate frozen free-surface profile. Dotted curve1610represents an approximate frozen free-surface extension above the source1604. The frozen free-surface extension1610may be calculated from a free-surface model based on parameters associated with the weather conditions measured at the time of the marine survey. For example, a Pierson-Moskowitz model of the free-surface may be used to calculate the frozen free-surface extension1610. The Pierson-Moskowitz model of a free-surface is based on the wind blowing steadily for a long period of time over a large free-surface area to produce waves that eventually reach a state of equilibrium. This condition is referred to as a “fully developed sea.” The Pierson-Moskowitz model used to calculate an extension to the approximate frozen free-surface profile at a point x in the x-direction is given by:

where for the integer index q≧0,

⁢F⁡(Kq)=2⁢π⁢⁢LW⁡(Kq)⁢{[N⁡(0,1)+j⁢⁢N⁡(0,1)]⁢/⁢2for⁢⁢i≠0,Q⁢/⁢2N⁡(0,1)for⁢⁢i=0,Q⁢/⁢2⁢⁢⁢and⁢⁢for⁢⁢q<0,F⁡(Kq)=F⁡(K-q)*.(16)
The parameter W(Kq) is the Pierson-Moskowitz spatial roughness spectrum, which for a fully developed sea surface in one-dimension (e.g., x-direction) is given by:

whereKqis the spatial wavenumber;Uwis the wind speed measured at a height of about 19 meters;α is 8.0×10−3;β is 0.74; andg is the acceleration due to gravity.
In Equations (15) and (16), the spatial wavenumber for component q is given by Kq=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 frozen in time Pierson-Moskowitz free-surface may be computed from Equation (15) using a FFT for computational efficiency.

The frozen free-surface extension fext(x) may be combined with the approximate frozen free-surface profile fint(x) to represent the frozen free-surface above the source and receivers by

f⁡(x)={fext⁡(x)for⁢⁢x0≤x<x1fint⁡(x)for⁢⁢x1≤x≤xl(18)
When extending the approximate frozen free-surface profile to approximate the frozen free-surface above the source1604, the approximate frozen free-surface profile f(x) reduces to fint(x).

In alternative implementations, the approximate frozen free-surface extension may be expanded to include a time parameter that characterizes the frozen free-surface at different times. Free-surface waves are generally dispersive and in deep water, and the frequency and wavenumber are related by a dispersion relation given by:
Ω(Kq)=√{square root over (gKq)}  (19)
Equation (19) implies that each spatial harmonic component of the free-surface wavefield may move with a definite phase velocity. As a result, in general, free-surface waves of longer wavelengths travel faster relative to waves with shorter wavelengths. Combining Equations (15) and (19) gives a time-varying frozen free-surface:

fext⁡(x,t)=1L⁢∑q=0Q-1⁢⁢F⁡(Kq)⁢ⅇⅈ⁢⁢(Kq⁢x-Ω⁡(Kq)⁢⁢t)(20)
where t is instantaneous time. Equation (20) characterizes a one-dimensional rough free-surface moving in the positive x-direction and may be used to compute the frozen free-surface extension1610at earlier or later times.

Consider a free-surface shape at an instant in time t with wave heights given by Equation (20). The wavenumber spectrum F(Kq) of the free-surface is computed according to Equation (16) and an arbitrary known dispersion relation Ω(Kq) is calculated according to Equation (19) may be used to calculate the frozen free-surface at an earlier (t−Δt) or a later (t+Δt) time by:

The frozen free-surface wavefield reflectivity (i.e., the response of a unit point source at the receiver position {right arrow over (r)}r) may be computed for a source at a location {right arrow over (x)}s=(xs,ys,zs) and a receiver at a location {right arrow over (x)}r=(xr,yr,zr) using:

where

Hn(1)⁡(x)≅2π⁢⁢x⁢exp⁡(ⅈ⁡(x-π4-n⁢⁢π2))
is the asymptotic form of the first-order Hankel function with n=0 and 1. The parameters of Equation (22) are represented inFIG. 17as follows:

k is the wavenumber of the propagating wavefield;

f(x) is the approximate frozen free-surface profile represented by dashed curve1702;

[x0,f(x0)] is a coordinate position1704of a running scattering point on the approximate frozen free-surface profile;

{right arrow over (x)}0is a vector1706from the origin of the Cartesian coordinate system to the running scattering point1704;

{right arrow over (r)}ris a vector1708from the origin to a receiver1710;

{right arrow over (x)}0−{right arrow over (x)}ris a vector1712from the running scattering point1704to the receiver1710;

{right arrow over (x)}dis a vector1714from the origin to a source1716;

{right arrow over (x)}dis a vector1718from the source1716to the receiver1710;

{right arrow over (ρ)} is a vector1720from the source1716to the scattering point1704;

Sris the path of the streamer which may be interpolated from the depth measurements obtained from the depth measuring devices location along the streamer; and

η(x′)={circumflex over (n)}·{circumflex over (ρ)}=cos θ is the obliquity factor with normal vectors {circumflex over (n)}1722and {circumflex over (ρ)}1724corresponding to the frozen free-surface normal and the unit vector direction of the incident field at [x0,f(x0)]1704and θ is the angle between the vectors {circumflex over (n)}1722and {circumflex over (ρ)}1724.

Next, the gradient of the pressure wavefield in the frequency domain, denoted by ∇rP, may be computed. Computation of the gradient ∇rP generally depends on the depth of the source in relation to the depth of the streamer.FIG. 18Ashows an example side-elevation view of the source1604located shallower than the streamer702, andFIG. 18Bshows an example side-elevation view of the source1604located deeper than the streamer702. InFIGS. 18A-18B, vectors, such as vector {right arrow over (x)}1802, represent coordinates (x,y,z). Subscript-r vectors, such as vector {right arrow over (r)}r1804, represent receiver coordinates (xr,yr,zr), and subscript-s vectors, such as vector {right arrow over (r)}sinFIG. 18Aand vector {right arrow over (x)}s1808inFIG. 18Brepresent source coordinates (xs,ys,zs). When the depth of the source1604is less than the depth of the streamer, as shown inFIG. 18A, the gradient of the pressure wavefield may be calculate by solving the following integral equation for ∇rP({right arrow over (x)}r,ω) in the frequency domain:
∫SrdSr{right arrow over (n)}·R({right arrow over (x)}r,{right arrow over (x)})∇rP({right arrow over (r)}r,ω)=a(ω)R({right arrow over (x)}s,{right arrow over (x)})+∫SrdSr{right arrow over (n)}·P({right arrow over (x)}r,ω)∇rR({right arrow over (x)}r,{right arrow over (x)})  (23)

where

a(ω) is the Fourier transform of the source-time function for the source at the coordinate location {right arrow over (x)}s;

P({right arrow over (x)}r,ω) is the pressure wavefield data in the frequency domain obtained from transforming the pressure wavefield data p({right arrow over (x)}r,t) using an FFT or a DFT.

R({right arrow over (x)},{right arrow over (x)}s) is the frozen free-surface reflectivity calculated according to Equation (22);

∇rR({right arrow over (x)}r,{right arrow over (x)}) is the radient of the reflectivity; and ∇rR({right arrow over (x)}r,ω) is the gradient of the pressure wavefield at the receiver.

Equation (23) is a Fredholm integral equation of the first kind for the gradient of the pressure wavefield ∇rP({right arrow over (x)}r,ω) where the right-hand side of Equation (23) contains only known parameters such as the pressure wavefield P({right arrow over (x)}r,ω) and the free-surface wavefield reflectivity R({right arrow over (x)}r,{right arrow over (x)}). The gradient of the reflectivity, ∇rR({right arrow over (x)}r,{right arrow over (x)}), may be computed using numerical gradient techniques applied to the free-surface wavefield reflectivity R({right arrow over (x)}r,{right arrow over (x)}) in Equation (22). On the other hand, when the source is located at a depth below the streamer Sr, as shown inFIG. 18B, the expression used to calculate the gradient of the pressure wavefield ∇rP({right arrow over (x)}r,ω) over the frequency range is give by:
∫SrdSr{right arrow over (n)}·R({right arrow over (x)}r,{right arrow over (x)})∇rP({right arrow over (x)}r,ω)=∫Sr{right arrow over (n)}·P({right arrow over (x)}r,ω)∇rR({right arrow over (x)}r,{right arrow over (x)})  (24)
In Equation (24), the source function a(ω) is zero. Note that the solutions of Equations (23) and (24) become unstable when the spectrum of the pressure wavefield has very small values (e.g., close to receiver ghost notches).

Depending on whether the source is at a depth above or below the streamer, as shown inFIGS. 18A and 18B, respectively, Equations (23) and (24) may be solved numerically for the gradient of the pressure wavefield, ∇rP({right arrow over (x)}r,ω) at receiver locations along the streamer using quadrature or expansion methods. For quadrature methods, the integrals may be approximated by quadrature formulas and the resulting system of algebraic equations is solved. For expansion methods, the solution may be approximated by an expansion in terms of basis functions.

An approximate normal particle velocity at each receiver location along the streamer may be calculated according to:

where

{right arrow over (n)} is a normal vector at a receiver; and

{right arrow over (n)}·∇rP({right arrow over (x)}r,ω) is the normal derivative of the pressure wavefield P at the receiver.

FIG. 19shows a segment of a streamer1902located beneath an approximate frozen free-surface profile1904in the xz-plane. A normal vector1906to the streamer1902at the receiver1908is given by:

n⇀=(nx,nz)=(-sin⁢⁢ϕ,cos⁢⁢ϕ)=(-ⅆzrⅆx,ⅆxrⅆq)(26)
The resulting approximate vertical particle velocity for each receiver 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,ω)  (27)
The approximate vertical particle velocity {tilde over (V)}z({right arrow over (x)}r,ω) may be transformed from the space-frequency domain to the f-k domain using an FFT or DFT to obtain {tilde over (V)}z({right arrow over (k)},ω), where {right arrow over (k)} is the wavevector (i.e., {right arrow over (k)}=(kx,ky,kz)), which may be used to compute separate approximate up-going and down-going pressure wavefields according to Equations (4a) and (4b).

FIG. 20shows a flow-control diagram of a method that computes separated wavefields from measured pressure wavefields. In block2001, a gather of pressure wavefield data is received. The gather may be a common-shot gather of pressure wavefield data obtained from receivers located along a streamer of a data acquisition surface, as described above with reference toFIG. 3. In block2002, the gather may be transformed from the s-t domain to the f-k domain as described above with reference to Equation (5). A for-loop beginning with block2003repeats the computational operations of blocks2004-2014for N trial depths. In block2004, a forward extrapolated wavefield computed for a trial depth zn, as described above with reference to Equation (6). In block2005, the forward extrapolated wavefield is transformed from f-k domain back to the s-t domain as described above with reference to Equation (8a) using an FFT or a DFT to obtain a forward extrapolated wavefield in the s-t domain. In block2006, a backward extrapolated wavefield computed for a trial depth zn,as described above with reference to Equation (6). In block2007, the backward extrapolated wavefield is transformed from f-k domain back to the s-t domain as described above with reference to Equation (8a) using an FFT or a DFT to obtain a backward extrapolated wavefield in the s-t domain. In block2008, a difference gather is computed from the forward and back extrapolated wavefields as described above with reference to Equation (10) and with reference toFIG. 11. A for-loop beginning with block2009repeats the computational operations of blocks2010-2012for a series of time windows as described above with reference toFIGS. 12A-12C. In block2010, a routine “compute maximum energy difference” is called as described below with reference toFIG. 21. In decision block2011, when window index m equals the number M of windows in the time window series, control flows to decision block2012. Otherwise, control flows to block2012and the window index is incremented and the operation represented by block2010is repeated. In decision block2013, when trial depth index n equals the number N of trial depth in the series of trial depths, control flows to block2015. Otherwise, the method increments the trial depth index n in block2014and repeats the operations represented by blocks2004-2012. In block2015, a routine “wavefield separation” is called as described below with reference toFIG. 22.

In an alternative implementation to the method represented inFIG. 20, the computational operation represented by block2008may be omitted. Rather than computing a difference gather from the forward and backward extrapolated wavefields, difference amplitudes may be calculate for each time of the series of time windows after block2009, as described above with reference toFIG. 15.

FIG. 21shows a flow-control diagram of the routine “compute maximum energy difference” called in block2010ofFIG. 20. In block2101, amplitude differences in a time window Wmare identified. In block2102, energy differences are calculated for each amplitude difference in the time window Wm, as described above with reference toFIG. 12BandFIG. 15. In block2103, a maximum energy difference for the time window Wmis identified as described above with reference to Equation (12).

FIG. 22shows a flow-control diagram of the routine “wavefield separation” called in block2015ofFIG. 20. A for-loop beginning with block2201repeats the computational operations represented by blocks2202-2209for each time window. A for-loop beginning with block2202repeats the computational operations represented by blocks2203-2205for each trial depth. In block2203, a set of maximum energy differences is formed as described above with reference toFIG. 13. In block2204, when the trial depth index n does not equal the number of trial depths, control flows to block2205in which the trial depth index is incremented and the operation represented by block2203is repeated. Otherwise, control flows to block2206and a peak trial depth is determined as described above with reference to Equation (13). In block2207, receiver coordinates associated with the peak trial depth are determined as described above with reference toFIG. 14. In decision block2208, when the time window index in does not equal the number of time windows in the time window series M, control flows to block2209and the time window index is incremented and the operations represented by blocks2202-2208are repeated. In block2209, a routine “compute vertical particle velocity wavefield” is called to compute an approximate vertical particle velocity wavefield. In block2210, separate wavefields are computed from the measured pressure wavefield received in block2001ofFIG. 20and the approximate vertical particle velocity wavefield.

FIG. 23shows a flow-control diagram of the routine “compute vertical particle velocity wavefield” called in block2209ofFIG. 23. In block2301, an approximate frozen free-surface profile that represents a frozen in time profile of the free-surface above a streamer is computed as described above with reference toFIG. 16A. In block2302, an approximate frozen free-surface extension is computed as described with reference toFIG. 16Band Equations (15)-(21). A for-loop beginning with block2303repeats the computational operations represented by blocks2304-2307for each receiver. In block2304, frozen free-surface wavefield reflectivity is computed as described above with reference to Equation (22) and shown inFIG. 17. In block2305, a pressure gradient is computed as described above with reference to Equation (23) when the source is located shallower than the streamer or computed as described above with reference to Equation (24) when the source is located deeper than the streamer. In block2306, a vertical particle velocity is computed as describe above with reference to Equation (25). In decision block2307, when a vertical particle velocity has not been computed for the full set of receivers that generated the measured pressure wavefield received in block2001ofFIG. 20, the operations represented by blocks2304-2306are repeated.

The mathematical equations and gathers presented above are not, in any way, intended to mean or suggest an abstract idea or concept. Instead the mathematical equations and gathers described above represent actual physical and concrete concepts and properties of materials that are in existence. The mathematical equations and methods described above are ultimately implemented on physical computer hardware, data-storage devices, and communications systems in order to obtain results that also represent physical and concrete concepts of materials that are in existence. 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 wavefield data that is physically recorded on physical data storage devices and undergoes computational processing using hardware as describe above to obtain up-going wavefield data that represents physical and concrete up-going pressure and vertical particle velocity wavefields. The up-going 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.

Note that in certain implementations, the computational operations represented blocks2004-2007may be executed in parallel as shown inFIG. 20. In other implementations, the computational operation represented by block2004may be executed before the computational operation represented by block2006followed a transformation of the gathers obtained in blocks2004and2006from the f-k domain to the s-t domain in block2401as shown inFIG. 24A. In still of the implementations, the computational operation represented by block2004may be executed after the computational operation represented by block2006followed a transformation of the gathers obtained in blocks2004and2006from the f-k domain to the s-t domain in block2401as shown inFIG. 24B.

FIG. 25shows an example of a computer system programmed to execute efficient methods of computing approximate vertical particle velocity wavefields and separate up-going and down-going pressure wavefields based on measured pressure wavefields and therefore represents a geophysical-analysis data-processing system. The internal components of many small, mid-sized, and large computer systems as well as specialized processor-based storage systems may be described with respect to this generalized architecture, although each particular system may feature many additional components, subsystems, and similar, parallel systems with architectures similar to this generalized architecture. The computer system contains one or multiple central processing units (“CPUs”)2502-2505, one or more electronic memories2508interconnected with the CPUs by a CPU/memory-subsystem bus2510or multiple busses, a first bridge2512that interconnects the CPU/memory-subsystem bus2510with additional busses2514and2516, or other types of high-speed interconnection media, including multiple, high-speed serial interconnects. The busses or serial interconnections, in turn, connect the CPUs and memory with specialized processors, such as a graphics processor2518, and with one or more additional bridges2520, which are interconnected with high-speed serial links or with multiple controllers2522-2527, such as controller2527, that provide access to various different types of computer-readable media, such as computer-readable medium2528, electronic displays, input devices, and other such components, subcomponents, and computational resources. The electronic displays, including visual display screen, audio speakers, and other output interfaces, and the input devices, including mice, keyboards, touch screens, and other such input interfaces, together constitute input and output interfaces that allow the computer system to interact with human users. Computer-readable medium2528is a data-storage device, including electronic memory, optical or magnetic disk drive, USB drive, flash memory and other such data-storage device. The computer-readable medium2528can be used to store machine-readable instructions that encode the computational methods described above and can be used to store encoded data, during store operations, and from which encoded data can be retrieved, during read operations, by computer systems, data-storage systems, and peripheral devices.

FIGS. 26A-26Bshow simulation results in calculating an approximation to a frozen free-surface.FIG. 26Ashows a total pressure wavefield obtained from an actual numerical simulation. Horizontal axis2602represents channel (i.e., trace) numbers and vertical axis2604represents time. Dashed hyperbolic-shaped curves2606and2608identify lower and upper the hyperbolic-shaped time series boundaries. A series of time windows was applied to pressure wavefield data within the boundaries2606and2608to compute an approximate free-surface represented inFIG. 26B.FIG. 26Bshows as a true free-surface represented by solid-line curve2610and an approximate frozen free-surface profile represented by dashed-line curve2612. Horizontal axis2614represents distance along a streamer and vertical axis2616represents depth. Methods described above were applied to the synthetic pressure wavefield located between the boundaries2606and2608to compute the approximate frozen free-surface profile2612, which shows a close approximation to the true free-surface distances below about 4500 meters but still tracks the overall profile of the true free-surface for distances greater than about 4500 meters.

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 up-going and down-going wavefield gathers computed as described above form a geophysical data product indicative of certain properties of a subterranean formation. The geophysical data product may include processed seismic data and may be stored on a computer-readable medium as described above. The geophysical data product may be produced offshore (i.e. by equipment on the survey vessel102) 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 data-storage facility in the United States. Once onshore in the United States, geophysical analysis may be performed on the data product.