GEOPHYSICS, VOL. 75, NO. 6 共NOVEMBER-DECEMBER 2010兲; P. WB121–WB136, 18 FIGS. 10.1190/1.3478375

Data reconstruction with shot-profile least-squares migration

Sam T. Kaplan1, Mostafa Naghizadeh2, and Mauricio D. Sacchi1

ABSTRACT We introduce shot-profile migration data reconstruction 共SPDR兲. SPDR constructs a least-squares migrated shot gather using shot-profile migration and demigration operators. Both operators are constructed with a constant migration velocity model for efficiency and so that SPDR requires minimal information about the underlying geology.Applying the demigration operator to the least-squares migrated shot gather gives the reconstructed data gather. SPDR can reconstruct a shot gather from observed data that are spatially aliased. Given a constraint on the geological dips in an approximate model of the earth’s reflector, signal and aliased energy that interfere in the common shot data gather are disjoint in the migrated shot gather. In the least-squares migration algorithm, we construct weights to take advantage of this separation, suppressing the aliased energy while retaining the signal, and allowing SPDR to reconstruct a shot gather from aliased data. SPDR is illustrated with synthetic data examples and one real data example.

INTRODUCTION We introduce a new data reconstruction algorithm, shot-profile migration data reconstruction 共SPDR兲. SPDR uses a constant migration velocity model for efficiency and so that minimal assumptions are made about earth structure.At the core of least-squares migration are forward 共demigration兲 and adjoint 共migration兲 operators, the former mapping from model space to data space, and the latter mapping from data space to model space. SPDR uses least-squares migration to find an optimal instance of model space which, in turn, is mapped to data space using demigration, providing a reconstructed shot gather. SPDR applies shot-profile migration and demigration operators one shot gather at a time. Moreover, we assume a 2D ex-

-periment where SPDR is a 1D data reconstruction algorithm interpolating traces in one geophone dimension for any given shot gather. In seismic data reconstruction, algorithms tend to fall into two categories, being rooted in either signal processing or the wave equation. Examples of the former include f ⳮ x trace interpolation 共Spitz, 1991兲, f ⳮ k trace interpolation 共Gülünay, 2003兲, minimum weighted norm interpolation 共Liu and Sacchi, 2004兲, reconstruction using the curvelet transform 共Herrmann and Hennenfent, 2008兲, and multistep autoregressive reconstruction 共Naghizadeh and Sacchi, 2007兲; whereas examples of the latter include the use the Born approximation and asymptotic methods 共Chiu and Stolt, 2002; Stolt, 2002; Miranda, 2004; Ramírez et al., 2006兲, Green’s theorem 共Ramírez and Weglein, 2009兲, poststack Stolt migration operators and datafitting methods 共Trad, 2003兲, and dip moveout operators 共e.g., Baumstein and Hadidi, 2006兲. SPDR belongs to the family of waveequation based methods for data reconstruction. It differs from previous efforts in its parameterization of model space, being based on shot-profile migration operators; whereas Stolt 共2002兲, Chiu and Stolt 共2002兲, and Ramírez et al. 共2006兲 use a source-receiver parameterization and Trad 共2003兲 uses a poststack geometry. Additionally, SPDR relies on data-fitting methods such as those used in Trad 共2003兲, rather than direct inversion and asymptotic approximation which are used, for example, in Stolt 共2002兲. A challenge in data reconstruction is aliased signal. In particular, when aliased energy is present and interferes with signal, their separation becomes challenging 共but not impossible兲. A recent example of data reconstruction is Naghizadeh and Sacchi 共2007兲, which uses the nonaliased part of data to aid in the reconstruction of the aliased part. An alternative approach is to transform data via some operator that maps from data space to some model space, and such that in that model space, the corresponding representation of signal and alias are separable. This is a common approach in signal processing 共e.g., Trad et al., 2003兲 and also is the approach that we take in SPDR. In particular, the SPDR model space is the constant velocity shotprofile migration of a single shot gather 共i.e., a migrated shot gather兲. This means that the model space is a representation of the earth’s reflectors parameterized by pseudodepth 共i.e., depth under the assumption of a constant migration velocity model兲 and lateral posi-

Manuscript received by the Editor 10 March 2010; revised manuscript received 1 June 2010; published online 22 December 2010. 1 University of Alberta, Department of Physics, Edmonton, Alberta, Canada. E-mail: [email protected]; [email protected].. 2 University of Calgary, Calgary, CREWES, Department of Geoscience, Alberta, Canada. E-mail: [email protected], [email protected]. © 2010 Society of Exploration Geophysicists. All rights reserved.

WB121

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB122

Kaplan et al.

tion. We will show that under the assumption of limited dips in the earth’s reflectors, the SPDR model space allows for the suppression of alias while preserving signal, thus allowing for the reconstruction of aliased data. We begin with a description of shot-profile migration and demigration built from the Born approximation to the acoustic wavefield and constant velocity Green’s functions. We apply shot-profile migration to an analytic example to illustrate its mapping of signal and alias from the data space 共shot gather兲 to the model space 共migrated shot gather兲. With this mapping, we infer that with constrained dip in the earth’s reflectors, the signal and alias map to disjoint regions of the model space. We note that if the dips in the earth’s reflectors are large, then the model-space representation of signal and alias are not necessarily disjoint, and the SPDR algorithm will fail unless survey parameters are adjusted to increase Nyquist wavenumbers. We reproduce the analytic result with a similar numerical result. Given the demigration 共forward兲 and migration 共adjoint兲 operators, we construct a set of weighted least-squares normal equations for leastsquares migration. The normal equations are built such that 共1兲 to some prescribed noise tolerance, the reconstructed shot gather fits the observed shot gather, and 共2兲 the portion of the model-space representative of aliased data is suppressed. Solving the normal equations gives an optimal migrated shot gather, and, in turn, the demigration of the optimal migrated shot gather is the reconstructed data 共shot gather兲. This summarizes the procedure followed in the SPDR method. We apply SPDR to synthetic data generated by a finitedifference implementation of the acoustic wave equation, and to a real data example from the Gulf of Mexico.

SHOT-PROFILE MIGRATION AND DEMIGRATION OPERATORS SPDR relies on least-squares migration 共e.g., Nemeth et al., 1999兲. In turn, least-squares migration depends on demigration and migration operators, the former providing a mapping from model space to data space 共the forward operator兲, and the latter being its adjoint, mapping from data space to model space. Within the context of the Born approximation, the SPDR model space is a first-order approximation to the scattering potential 共e.g., Weglein et al., 2003兲, and under the shot-profile parameterization of the operators, this approximation is a migrated shot gather. Here, we describe constant velocity shot-profile wave-equation migration and demigration operators. The use of shot-profile migration rather than source-receiver migration 共Biondi, 2003兲 means that we cannot use the fast Stolt migration algorithm 共Stolt, 1978兲. However, it leads to a shot-profile parameterization of model space which is integral to the SPDR algorithm; namely, that the model is parameterized by pseudodepth and lateral position. The term pseudodepth describes a depth axis that honors the migration velocity model, but does not necessarily reflect the true depth of the reflectors. The forward operator 共demigration兲 for SPDR models the scattered seismic wavefield using the Born approximation under the assumption of an acoustic and constant velocity Green’s function so that 共e.g., Weglein et al., 2003兲

⬁

冕冕 冉冊

␺ s共xg,zg兩xs,zs;␻ 兲 ⳱ f共␻ 兲

G0共xg,zg兩x⬘,z⬘;␻ 兲

ⳮ⬁

⫻

␻

2

␣ 共x⬘,z⬘兲G0共x⬘,z⬘兩xs,zs;␻ 兲dx⬘dz⬘,

c0

共1兲 where G0 is a Green’s function for constant acoustic wavespeed c0, and such that 共e.g., DeSanto, 1992兲 ⬁

1 G0共xg,zg兩x⬘,z⬘;␻ 兲 ⳱ 2␲

冕冉

ⳮ

ⳮ⬁

1 i4kgz

冊

⫻eⳮikgx·共x⬘ⳮxg兲eikgz兩zgⳮz⬘兩dkgx . 共2兲 In equation 1, ␣ is the first-order approximation to the scattering potential. Within the context of least-squares migration and SPDR, ␣ is the model 共a migrated shot gather兲. The forward operator in equation 1 describes the mapping from the approximate scattering potential ␣ 共the model space兲 to the scattered wavefield ␺ s 共the data space兲 recorded at geophone positions 共xg,zg兲 where xg ⳱ 共xg,y g兲, and due to the seismic source f共 ␻ 兲 located at 共xs,zs兲 where xs ⳱ 共xs,y s兲. Equation 1 integrates over all possible scattering points 共x⬘,z⬘兲 where x⬘ ⳱ 共x⬘,y ⬘兲. The vertical wavenumber kgz in equation 2 is given by the dispersion relation

kgz ⳱ sgn共␻ 兲

冑

␻2 ⳮ kgx · kgx, c20

共3兲

where kgx ⳱ 共kgx,kgy兲 are the lateral wavenumbers 共Fourier conjugate variables of xg ⳱ 共xg,y g兲兲. The dispersion relation will play a role in understanding the null space of the demigration operator which, in turn, will help with our understanding of where alias in the data maps to in the model. Substituting equation 2 into equation 1 with zg ⳱ zs ⳱ z0, and taking the support of ␣ to be below z0 共see Appendix A兲 gives ⬁

冉 冊冕

␻ ␺ s共xg,␻ ;xs兲 ⳱ c0

2

F*u p共kgx,z⬘,␻ 兲F

z0

⫻关F*u p共kgx,z⬘,␻ 兲g共kgx,xs,␻ 兲兴␣ 共xg,z⬘兲dz⬘, 共4兲 where

u p共kgx,z⬘,␻ 兲 ⳱ ⳮ

eikgz共z⬘ⳮz0兲 i4kgz

共5兲

and

g共kgx,xs,␻ 兲 ⳱ 2␲ f共␻ 兲eⳮikgx·xs .

共6兲

Equations 4–6 describe shot-profile demigration for a constant migration velocity model 共i.e., the forward operator兲. In equation 4, F is the 2D Fourier transform, mapping from lateral space xg to lateral wavenumbers kgx. We use symmetric 共1 / 冑2␲ 兲 normalization in our Fourier transforms so that F*, the adjoint operator of F, also is its corresponding inverse Fourier transform. The function u p in equa-

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR tion 5 is a phase-shift operator, propagating the wavefield toward, or away from, potential scattering points. In equation 5, the vertical wavenumber kgz is given by the dispersion relation in equation 3. Finally, equation 6 is the Fourier representation of the seismic point source located laterally at x ⳱ xs. For the construction of the forward operator in equation 4, we have made some important choices. First, we derived the forward operator such that the approximate scattering potential ␣ has a shotprofile parameterization. In other words, in ␣ 共xg,z⬘兲, xg are the lateral coordinates of the earth and z⬘ is pseudodepth. This means that the model space is characterized by the subsurface reflectors that fall within the aperture of the shot. We will use this parameterization in SPDR; in fact, it is integral. The shot-profile parameterization of the scattering potential does not allow us to use Stolt migration, thus producing a more expensive data reconstruction algorithm than the one derived, for instance, in Trad 共2003兲. The adjoint of the demigration operator 共equation 4兲 is migration. We find it by recognizing a Fredholm integral equation of the first kind in the demigration operator, so that equation 4 can be written as

WB123

The adjoint indicates how energy distributes in the model space given its computation from some instance of data. Thus, before using the forward and adjoint operators in least-squares migration, we will analyze the adjoint to illustrate how signal and alias are mapped from the data space to the model space. This analysis will, in turn, motivate the construction of the weights in the weighted leastsquares migration algorithm.

OBSERVED DATA ON A NOMINAL GRID SPDR uses least-squares migration for data reconstruction. The input to the least-squares migration algorithm is observed seismic data 共shot gather兲, and is sampled from the scattered wavefield ␺ s共xg,␻ ;xs兲. In practice, we have some finite set of geophones to represent ␺ s in its lateral dimensions xg. In this paper, we assume 2D earth models and a single shot gather so that the data spaces and model spaces are sufficiently described by a single lateral dimension xg. Then, we define a regular nominal grid xng for that same lateral dimension so that

⬁

␺ s共xg,␻ ;xs兲 ⳱

冕

u共␻ ,z⬘;xg,xs兲␣ 共xg,z⬘兲dz⬘,

共7兲

ⳮ⬁

where

u共␻ ,z⬘;xg,xs兲 ⳱

冉冊

␻ 2 F*u p共kgx,z⬘,␻ 兲F c0

⫻关F*u p共kgx,z⬘,␻ 兲g共kgx,xs,␻ 兲兴.

共8兲

u*共␻ ,z;xg,xs兲␺ s共xg,␻ ;xs兲d␻ ,

共9兲

Then, ⬁

␣ †共xg,z兲 ⳱

冕

xng 苸 兵xg兩xg ⳱ k⌬xng,k ⳱ 1 . . . nng其.

共12兲

Further, we define an observed grid xog so that xog 苸 xng. This means that each sample on the observed grid also is a sample on the nominal grid. Further, in the data space, each sample of the observed grid corresponds to an observation 共i.e., a seismic trace兲. This is a common idea used in data reconstruction algorithms 共e.g., Liu and Sacchi, 2004兲 and is illustrated in Figure 1. In this paper, we assume that the observed grid is regular with sample spacing ⌬xog, although the case of an irregular xog grid also is interesting, but perhaps less challenging because it tends to increase the Nyquist wavenumber 共e.g.,Babu and Stoica, 2010兲 共We note, as an aside, that field data are, by the complications of acquisition, necessarily irregular兲. SPDR is equally applicable to regular and irregular observed grids.

ⳮ⬁

where

u*共␻ ,z;xg,xs兲 ⳱

ADJOINT MAPPING OF SIGNAL AND ALIAS DUE TO THE OBSERVED GRID

冉冊

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲F.

共10兲

In other words, substituting equation 10 into equation 9, we find ⬁

␣ †共xg,z;xs兲 ⳱

冕冉 冊

ⳮ⬁

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲F␺ s共xg,␻ ;xs兲d␻ ,

共11兲

where 共 · 兲* denotes the adjoint when applied to the Fourier operator, and the complex conjugate otherwise. Equation 11 describes shotprofile migration for an acoustic and constant migration velocity model. With equations 4 and 11, we have the necessary tools to construct a constant velocity shot-profile least-squares migration algorithm. An efficient numerical implementation of the forward and adjoint operators in equations 4 and 11 is given in Appendix B. For the purpose of least-squares migration, ␺ s is called the data and ␣ † is the migrated shot gather found by applying the adjoint operator to the data.

In this section, we analyze the adjoint 共shot-profile migration兲 operator with an analytic example. We let data ␺ s be the sampling 共Dirac-comb兲 function corresponding to a regular observed grid xog. Our goal is to provide some understanding of where signal and alias in the data space map to in the model space via the adjoint operator. We begin with a data-space description of the signal and alias using the sampling theorem. Then, we find the adjoint 共model-space representation兲 of the sampling function by applying the shot-profile migration operator. This provides a basis for categorizing the energy in the model space as corresponding to either the data-space signal or the data-space alias; that is, as either model-space signal or modelspace artifacts resulting from aliased data. ∆xn

∆xo

xg

Figure 1. Schematic example of nominal and observed grids. The boxes 共open and filled兲 denote the nominal grid, and the observed grid is denoted by the filled boxes. ⌬xng denotes the grid spacing for the nominal grid, and ⌬x°g denotes the grid spacing for a regular observed grid.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB124

Kaplan et al.

¯s In the sampling theorem, we represent ␺ s共xg,␻ ;xs兲 with ␺ ⳱ IIIo共xg兲␺ s where IIIo共xg兲 is a sampling function such that

⬁

¯ †共xg,z;xs兲 ⳱ ␣

nog

冕冉 冊

ⳮ⬁

IIIo共xg兲 ⳱ 兺 ␦ 共xg Ⳮ k⌬xog兲.

共13兲

⫻F*u*p 共kgx,z,␻ 兲␺ s共kgx,␻ ;xs兲d␻

k⳱0

⬁

Then, a well known result of the sampling theorem states that 共e.g., Lathi, 1998, p. 319兲

兺

ⳮ⬁

␺ s共kgx Ⳮ 2␲ p/⌬xog,␻ ;xs兲.

共14兲

⬁

The series in equation 14 contains signal for p ⳱ 0 and alias for p ⫽ 0. We notice that equation 14 represents signal and alias in the data space. We are interested in their manifestation in the model space. For this purpose, we consider the migration operator in equation 11 ¯ s in equation 14, so that using the Fourier representation applied to ␺ ¯ s共kgx,␻ ;xs兲, of the scattered wavefield ␺

冕冉 冊

Ⳮ

ⳮ⬁

⬁

冕冉 冊

Ⳮ

冕冉 冊

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

¯ 共k ,␻ ;x 兲d␻ . ⫻F*u*p 共kgx,z,␻ 兲␺ s gx s

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲␺ s共kgx ⳮ 2␲ /⌬xog,␻ ;xs兲d␻

⬁

ⳮ⬁

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲␺ s共kgx Ⳮ 2␲ /⌬xog,␻ ;xs兲d␻

p⳱ⳮ⬁

¯ †共xg,z;xs兲 ⳱ ␣

冕冉 冊

Ⳮ

⬁

¯ 共k ,␻ ;x 兲 ⳱ ␺ s gx s

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

ⳮ⬁

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲␺ s共kgx Ⳮ 4␲ /⌬xog,␻ ;xs兲d␻ Ⳮ¯

共15兲

⬁

Then, substitution of equation 14 into equation 15 gives

⳱

a)

⬁

冕

␣ †0共xg,z,␻ ;xs兲d␻ Ⳮ

ⳮ⬁

冕

␣ †1共xg,z,␻ ;xs兲d␻

ⳮ⬁ ⬁

Ⳮ

冕

† ␣ ⳮ1 共xg,z,␻ ;xs兲d␻

ⳮ⬁

0

⬁

k = 0 k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 k = 7 k = 8 k = 9 k = 10 Lateral position

Ⳮ

b)

冕

␣ †2共xg,z,␻ ;xs兲d␻ Ⳮ ¯ ,

共16兲

ⳮ⬁ − cω0

ω c0

where

0 p = –2

Depth wavenumber

c)

p = –1 p=0 p=1 Lateral wavenumber

␣ †p共xg,z,␻ ;xs兲 ⳱

p=2

冉冊

␻ 2 关F*u*p 共kgx,z,␻ 兲g*共kgx,xs,␻ 兲兴 c0

⫻F*u*p 共kgx,z,␻ 兲␺ s共kgx Ⳮ 2␲ p/⌬xog,␻ ;xs兲.

0.0

− cω0 −

共17兲

2π ∆xg

− cω0

ω c0

ω c0

+

2π ∆xg

–0.1 p = –1 p=0 p=1 Lateral wavenumber

Figure 2. Analytic example for a single frequency: 共a兲 the sampling function for k ⳱ 0 . . . 10, 共b兲 the Fourier representation of the sampling function, and 共c兲 the adjoint of the sampling function in lateral kgx and depth kz wavenumber coordinates.

The term 兰␣ †0共xg,z,␻ ;xs兲d␻ from equation 16 is the hypothetical image that would be obtained from a continuously sampled spatial dimension xg in data 共i.e., ␣ †0 ⳱ ␣ †兲, and ␣ †p, p ⫽ 0 is the expression of the alias in the model space 共i.e., the model-space artifacts兲. We illustrate the model-space artifacts for a single frequency ¯ s共xg,␻ ;xs兲 ⳱ IIIo共xg兲␺ 0共 ␻ 兲, ¯ †. In particular, we let ␺ component of ␣ which is shown in Figure 2a for k ⳱ 0 . . . 10. The discrete Fourier transform of IIIo共xg兲 is 共e.g., Lathi, 1998, p. 320兲

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR

IIIo共kgx兲 ⳱

1

冉

⬁

兺

⌬xog p⳱ⳮ⬁

␦ kgx Ⳮ

2␲ p ⌬xog

冊

共18兲

¯ s共xg,␻ ;xs兲 ⳱ IIIo共xg兲␺ 0共 ␻ 兲, and is shown in Figure 2b. Next, we let ␺ so that from equation 14

␺ s共kgx Ⳮ 2␲ p/⌬xog,␻ ;xs兲 ⳱

1

␦ ⌬xog

冉

kgx Ⳮ

2␲ p ⌬xog

冊

␺ 0共␻ 兲. 共19兲

Upon substitution of equation 19 into equation 17, we find for any given frequency ␻ 共see Appendix C兲

␣ †p共kgx,kz,␻ ;xs ⳱ 0兲 ⳱

冉冊

f *共 ␻ 兲 ␺ 0共 ␻ 兲 ␻ 冑2␲ ⌬xog c0

2

␦ 共kz Ⳮ 共kgz共p兲

Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲, where k p ⳱ 2␲ p / ⌬xog,

kgz共p兲 ⳱ sgn共␻ 兲

冑

再

o
共x 兲 ⳱ 1, xg ⳱ k⌬xg,k ⳱ 0 . . . 9 III o g 0, otherwise,

␻2 2 2 ⳮ k p, c0

共20兲

共21兲

共22兲

Equation 22 is shown in Figure 2c for p ⳱ 0,Ⳳ1. In SPDR, we let the null space of the migration operator include the evanescent portion of the wavefield which, in turn, is defined for when kgz in the dispersion relation 共equation 3兲 is imaginary. This is common practice in wave-equation migration algorithms, although there are exceptions 共Huang et al., 1999兲. The nonevanescent portion of the wavefield is defined for when

ⳮ

␻ ␻ ⱕ kgx ⱕ . c0 c0

Lateral position (m)

共23兲

In Figure 2b, we plot the bounds of the nonevanescent region under some assumed frequency ␻ and migration wavespeed c0. Notice that energy for p ⳱ 0,Ⳳ1 falls within these bounds, and all other energy falls within the evanescent portion of the wavefield. Hence, for this example, the adjoint operator will annihilate energy mapped from the data space for where 兩p兩 ⱖ 2. This is exactly why in Figure 2c we show ␣ †p for only p ⳱ 0,Ⳳ1. In the model space, these bounds are translated by k p in equation 20 so that

ⳮ

␻ ␻ Ⳮ k p ⱕ kgx ⱕ Ⳮ k p . c0 c0

In the previous section, we explored the adjoint mapping of the sampling function. The example showed how the sampling function in data space maps to model space. Moreover, the resulting pictures show obvious overlap between the model-space representation of signal and alias 共e.g., Figure 4兲. This overlap would make it difficult to construct an algorithm that suppresses the model-space artifacts while preserving the model-space signal, and this is exactly how we intend SPDR to function. However, we remind ourselves that the SPDR model space provides a representation of the earth’s reflectors. This is a direct consequence of using shot-profile migration for

共24兲

In particular, we show in Figure 2c the lower and upper bounds for ␣ †p共kgx,kz,␻ ;xs ⳱ 0兲 for p ⳱ 0, the lower bound for p ⳱ ⳮ1, and the upper bound for p ⳱ 1. The synthetic experiment shown in Figure 3 confirms the analytic example. In particular, we define a nominal grid xng with ⌬xng ⳱ 5 m. ˆ o共xg兲 for xg 苸 xn and where In Figure 3a, we plot ␺ s共xg,␻ ;xs兲 ⳱ III g

a)

0

200

400

600

800

1000 1200 1400 1600 1800 2000

0 Time (s)

␻2 ⳮ 共kgx ⳮ k p兲2 . c20

REFLECTOR DIP CONSTRAINTS ON THE MODEL SPACE

1

2

b)

Lateral wavenumber (1/m) –0.06 0.00

–0.04

–0.02

0

0.02

0.04

0.06

Depth wavenumber (1/m)

冑

共25兲

and where ⌬xog ⳱ 300 m. In other words, the observed grid is regular ˆ o共xg兲 with samples spaced every 300 m. We show the adjoint of III evaluated for a temporal frequency of 6 Hz in Figure 3b. The resulting pattern mimics what is predicted by the analytic example. The examples presented in Figures 2 and 3 were for a single frequency. In Figure 4, we show results for when five frequencies, equally spaced between 6 and 10 Hz, contribute to the integral in the adjoint operator 共equation 11兲. In particular, Figure 4b shows the analytic result in the model space, whereas Figure 4a shows a corresponding synthetic example. The analytic example in this section gives us an understanding of where signal and alias in the data-space map in the model space. This understanding is useful in the development of SPDR. However, we point out that the application of the adjoint operator to the sampling function is equivalent to propagating a receiver side horizontal plane wave into the earth and then to applying the downward continued source, and the shot-profile migration imaging condition. In practice, the data will contain some continuum of plane-wave components, and each will contribute to the distribution of energy in the model space.

kz is the Fourier conjugate variable of z, and kgz was given by the dispersion relation in equation 3. Equation 20 means that ␣ †p is nonzero only when

kz ⳱ ⳮ kgz共p兲 ⳮ sgn共␻ 兲

冎

WB125

–0.1

Figure 3. Synthetic example for a single frequency: 共a兲 the sampling function, 共b兲 the adjoint of the sampling function in lateral kgx and depth kz wavenumber coordinates.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB126

Kaplan et al.

the adjoint operator. Thus, it follows that for any given frequency, the adjoint will be band-limited in kgx according to the maximum dip of the reflectors 共parameterized by lateral position and pseudodepth兲. We postulate that this band limitation allows for a disjoint representation of the model-space signal and artifacts. In this section, we illustrate this idea with two examples. First, we apply a band-limitation operator to the analytic example studied in the previous section. Second, we use a synthetic example in which the earth model consists of a single reflector. In the adjoint of the sampling function, we assume that ␣ †p is bandlimited so that equation 20 becomes

␣ †p共kgx,kz,␻ ;xs ⳱ 0兲 ⳱ h共kgx,kz兲

冉冊

f *共 ␻ 兲 ␺ 0共 ␻ 兲 ␻ 冑2␲ ⌬xog c0

2

␦ 共kz

Ⳮ 共kgz共p兲 Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲, 共26兲 where

h共kgx,kz兲 ⳱

再

1, 兩kgx兩 ⱕ kb共kz兲 0, 兩kgx兩 ⬎ kb共kz兲,

冎

共27兲

and where kb共kz兲 is built from the expected maximum dip of the earth’s reflectors in pseudodepth. In particular, if a reflector has dip ␰ such that

␣ 共xg,z兲 ⳱ ␦ 共z ⳮ ␰ xg兲,

共28兲

then it can be shown that

␣ 共kgx,kz兲 ⳱ 2␲ ␦ 共kgx Ⳮ ␰ kz兲.

共29兲

Hence, it must be that kb共kz兲 ⳱ ␰ kz. For example, when the dip of the reflector is null, so is kb, and likewise, kb increases with increasing ␰ . Figure 5 illustrates equation 26 for small dip 共 ␰ ⳱ 0.1兲. It is a bandLateral wavenumber (1/m) –0.06

–0.02

0

0.02

0.04

0.06

0.0

Thus far, we have identified the adjoint operator 共shot-profile migration兲 as a useful tool for mapping data 共common shot data gather兲 with interfering data-space signal and alias to a model-space representation 共migrated shot gather兲 where signal and the artifacts caused by the alias are separated. To move toward a complete description of SPDR, our next task is to find a method that suppresses the artifacts from the model space while preserving the signal. For this purpose, we use weighted least-squares inversion, rather than

–0.1

Depth wavenumber (1/m)

b)

LEAST-SQUARES DATA RECONSTRUCTION

0.0

–0.1 –0.06

–0.04

–0.02

0

0.02

0.04

0.06

Lateral wavenumber (1/m)

Figure 4. Multifrequency example: the adjoint of the sampling function for 共a兲 the synthetic example, and 共b兲 the analytic example as expressed by equation 22. Five frequencies spaced evenly between 6 and 10 Hz are used in the integral over frequency 共equation 11兲. In 共b兲 we show energy for p ⳱ 0 共solid line兲, p ⳱ 1 共long dashed line兲, p ⳱ ⳮ1 共dash-dot line兲, p ⳱ Ⳳ2 共short dashed line兲.

Depth wavenumber (1/m)

Depth wavenumber (1/m)

a)

–0.04

limited version of Figure 4b. Notice that with the band-limiting function h共kgx,kz兲 applied, the model-space signal and artifacts of the sampling function are separable. This will make it possible to build a data reconstruction algorithm that filters out the artifacts but preserves the signal. To further illustrate the band limitation, we construct a synthetic example using the model in Figure 6a, and the finite-difference approximation to the acoustic wave equation, constructing data ␺ s in Figure 6b. The source 共f共 ␻ 兲 in equation 1兲 is a minimum-phase Ricker wavelet with a peak frequency of 15 Hz. For the data, we use a nominal grid spacing of 5 m, and a regular observation grid, with observations 共geophones兲 spaced every 300 m. The data are severely aliased as illustrated by their f ⳮ k spectrum in Figure 6c. Because the model consists of a single reflector with null dip, we expect the adjoint operator to separate signal and alias. In Figure 6d, we plot the adjoint of the data, using the shot-profile migration operator. A quick comparison to the band-limited analytic example in Figure 5 allows us to label the clusters of energy corresponding to the model-space signal and artifacts. In particular, we identify energy corresponding to p ⳱ 0,Ⳳ1,Ⳳ2. We can construct a lower bound for the maximum reflector dip ␰ resulting in separable model-space signal and artifacts, and allowing for the successful application of SPDR. This dip will be a function of the observed grid spacing 共i.e., the lateral Nyquist wavenumber兲. In particular, we note that equation 29 gives the relation kgx ⳱ ⳮ␰ kz so that max兩kgx兩 ⳱ ␰ max兩kz兩. At the very least, it must be that max兩kgx兩 ⬍ k1, where k1 ⳱ 2␲ / ⌬xog was defined in the previous section. Hence, it follows that ␰ max兩kz兩 ⬍ 2␲ / ⌬xog, or, ␰ ⬍ 2␲ / 共⌬xog max兩kz兩兲. This shows, for example, an inverse relation between the observed grid sampling interval and the maximum allowed reflector dip in the model space. In other words, it shows a direct relation between the lateral Nyquist wavenumber and the maximum allowed dip.

0.0

–0.1 –0.06

–0.04

–0.02

0

0.02

0.04

0.06

Lateral wavenumber (1/m)

Figure 5. Analytic and multifrequency example: We represent the adjoint of the sampling function under the band-limitation constraint introduced in equation 26 using five frequencies evenly spaced between 6 and 10 Hz.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR

WB127

Key to the construction of m*, and thus also to the reconstructed the adjoint, to construct a model. The consequences are twofold. First, it will ensure that the forward operator 共demigration兲 applied data d*, is the choice of the model and data weight matrices, Wm and to the model produces modeled data that fit the observed data. SecWd. The data weight matrix Wd is a diagonal matrix such that its ith ond, it allows us to penalize those regions of the model space that diagonal element is contain an expression of the data alias rather than signal. In the previous two sections, we identified disjoint regions of the model space that correspond to either the model-space signal or artifacts. We be1, i 苸 Id 关Wd兴ii ⳱ 共33兲 gin with a description of the weighted least-squares normal equa0, i 苸 Id, tions. Then, we construct data-weights that honor the observed data, and model weights that penalize the energy from the aliased data. where Id is the set of indices corresponding to the observed grid xog. To construct the weighted least-squares inversion, we begin with This ensures that the cost function in equation 30 works to fit only the some definitions. We let d be a data vector of length M realized from n n observed data, and not the missing data. The construction of model a single shot gather ␺ s共xg,␻ ;xs兲, and where we recall that xg deweights is slightly more involved, and draws on the discussion in the scribes the nominal grid for the lateral dimension xg. We let dobs also previous two sections. In particular, we define be a vector of length M. It is sampled by the nominal grid xng, and will be nonzero only for the intersection of the nominal xng and observed xog grids. For this single shot gather, M ⳱ n␻ nng where n␻ is the number of nonnegative frequency samples, and as before, nng is the number of grid points in the nominal grid Lateral position (m) Offset (m) xng. Likewise, we let m be a model vector of length a) b) 0 1000 2000 3000 4000 5000 –1200 –800 –400 0 400 800 1200 N realized from ␣ 共xng,z;xs兲, where N ⳱ nznng, and 0 0 nz are the number of samples in depth for the shotprofile migrated gather. Last, we let A be the M 250 ⫻ N matrix representation of the forward opera1500 2000 m/s tor 共shot-profile demigration in equation 4兲, and 1 500 AH its corresponding adjoint operator 共shot-profile migration in equation 11兲. Then, A maps from m to d. 750 To find an optimal m that honors the observed 2 data dobs, we find the minimum of the cost function,

␾ 共m兲 ⳱ 储Wd共dobs ⳮ Am兲储22 Ⳮ ␮储Wmm储22

H 共AHWH d WdA Ⳮ ␮WmWm兲m obs ⳱ A HW H d W dd ,

共32兲

which we solve to find an optimal scattering potential m ⳱ m*. Then, the data can be reconstructed using the wavefield modeling 共forward兲 operator such that d* ⳱ Am*. That is, d* is the SPDR reconstructed data. In practice, equation 32 is solved using the least-squares conjugate gradient method 共Paige and Saunders, 1982兲 and implicit construction of the matrices: A, AH, Wd, and Wm in equation 32.

Frequency (Hz)

where Wd are data weights and Wm are model weights. It is convenient to partition the cost function into two components ␾ d and ␾ m, with ␾ d being the data misfit function and ␾ m being the model-norm function. The parameter ␮ allows for a trade-off between fitting the observed data and honoring the model norm. Finding the minimum of equation 30 results in the normal equations

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

0

10

20

30

Lateral wavenumber (1/m)

d)

–0.06

–0.04

–0.02

0

0.02

0.04

0.06

0.0

Depth wavenumber (1/m)

共31兲

Lateral wavenumber (1/m)

c)

共30兲 ⳱ ␾ d共m兲 Ⳮ ␮␾ m共m兲,

冎

Time (s)

Depth (m)

再

–0.1

p = −1 –0.2

p = −2

p=1 p=0

p=2

Figure 6. Synthetic example for a zero-dip reflector: 共a兲 the earth model is a half-space over a half-space, 共b兲 synthetic data generated by the finite-difference method, 共c兲 f ⳮ k spectrum of the data, and 共d兲 kgx ⳮ kz spectrum of the adjoint. We identify energy corresponding to p ⳱ 0,Ⳳ1,Ⳳ2 in equation 17.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB128

Kaplan et al.

Wm ⳱ Fⳮ1WF,

共34兲

k1共kz兲 and k2共kz兲 are defined according to the model dip constraints in equation 29. Namely,

where F is the 2D Fourier transform over lateral position xg and depth z. Then, W penalizes the model-space artifacts according to the windowing function

k1共kz兲 ⳱ ␰ˆ 兩kz兩 Ⳮ ␩

共36兲

w共kgx,kz兲

k2共kz兲 ⳱ ␰ˆ 兩kz兩 Ⳮ ␩ Ⳮ ␶ ,

共37兲

冦

冧

冉 冉

kgx Ⳮ k2共kz兲

冉 冉冉

冊冊

冊冊冊

Model weights

Time (s)

1

kgx ⬍ ⳮ k2共kz兲

Frequency (Hz)

Depth wavenumber (1/m)

Time (s)

so that ␶ defines the length of a cosine taper, and ␩ is the expected extent of the model-space signal along the kgx dimension. The pa1 Ⳮ cos ␲ , 1Ⳮ ⳮ k2共kz兲 ⬍ kgx ⬍ ⳮ k1共kz兲 rameter ␰ˆ is the estimated maximum dip of the earth’s reflectors as ␶ 2⑀ parameterized by pseudodepth and lateral position in the migrated ⳱ 1, ⳮ k1共kz兲 ⬍ kgx ⬍ k1共kz兲 shot gather. For example, we could set ␰ˆ ⳱ 0 and ␩ ⳱ 2␲ / ⌬xog, the kgx ⳮ k1共kz兲 1 former for simplicity and the latter motivated by the sampling theo1 Ⳮ cos ␲ 1 Ⳮ , k1共kz兲 ⬍ kgx ⬍ k2共kz兲 1Ⳮ rem 共see equation 14兲. The parameter ␶ defines the length of the co2⑀ ␶ sine taper. The function kb共kz兲 was defined in equation 29, and dek2共kz兲 ⬍ kgx, 1 Ⳮ 1/⑀ , pends on some user-selected maximum dip of the earth’s reflectors. 共35兲 We note that equation 35 can be generalized to allow for specification of separate maximum positive and negative dips. where ⑀
1 is some small dimensionless constant. The functions To illustrate, we continue with the example shown in Figure 6. In particular, we use the decimated data from Figure 6b for dobs in equation 32. Then, we solve equaOffset (m) Offset (m) a) –1200 –800 –400 0 400 800 1200 b)–1200 –800 –400 0 400 800 1200 tion 32 for m*. We recall that m* is a discrete rep0 0 resentation of the approximate scattering potential ␣ ⳱ ␣ *共xg,z;xs兲 共migrated shot gather兲, and in Figure 7c we plot its amplitude in the wavenumber domain, 兩 ␣ *共kgx,kz;xs兲兩. This should be com1 1 pared to Figure 6d, where we plot the adjoint of ? ? ¯ †共kgx,kz;xs兲兩. Finally, we use the decimated data 兩 ␣ m* to compute the reconstructed data d* ⳱ Am*. 2 2 We show the reconstructed data in Figure 7b. For comparison, we show the same data computed using the finite-difference approximation to the Lateral wavenumber (1/m) acoustic wave equation in Figure 7a. In Figure 7d, c) –0.06 –0.04 –0.02 0 0.02 0.04 0.06 0.0 we plot the f ⳮ k spectrum of the reconstructed data which should be compared to the f ⳮ k spectrum of the decimated data which was shown in Figure 6c. –0.1 We make some comments about the results in Figure 7. First, in the recovered data, we notice 1/
some missing energy at the far-offsets 共indicated 1 by the arrows in Figure 7a兲. This energy corre–0.2 sponds to the seismic head wave, and belongs to the null space of our forward operator. Second, in Lateral wavenumber (1/m) d) Figure 7c, we plot the model weights. In this ex–0.6 –0.4 –0.2 0 0.2 0.4 0.6 0 ample, we set ␰ˆ ⳱ 0 so that the model weights are invariant to the vertical wavenumber kz. Finally, 10 we plot the data misfit ␾ d as a function of conjugate gradient iteration in Figure 8; it shows that the conjugate gradient algorithm used to solve 20 equation 32 converges after 60 iterations. One could use, for example, an ␹ 2 statistic built from 30 an estimate of the noise in the data to determine a stopping criterion for the conjugate gradient iterations 共e.g., Hansen, 1998兲. Numerical tests show Figure 7. Synthetic example for a zero-dip reflector: 共a兲 data computed using the acoustic that the algorithm tends to converge sufficiently finite-difference approximation to the acoustic wave equation on the nominal grid. The after tens of iterations 共of course, this will depend head wave is indicated by two vertical arrows. 共b兲 The reconstructed data computed from the observed data in Figure 6b, and using SPDR. 共c兲 The wavenumber spectrum of the on the size of the model space兲. least-squares inverse m*. The model weights used in the least-squares migration normal Our description of SPDR now is complete. The equations are also shown in 共c兲. 共d兲 The f ⳮ k spectrum of the SPDR reconstructed data data weights in the least-squares inversion were d . 1 Ⳮ 1/⑀ ,

*

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR

increased, and the separation between the signal and artifacts degrades, in turn, reducing the effectiveness of SPDR. The location of the artifacts is, in part, dictated by the observed 0.01

Data misfit

built to honor the known data on the observed grid xog, while giving no weight to the remaining samples on the nominal grid xng 共i.e., the missing traces in the shot gather兲. The model weights were built to penalize the portion of the model space that corresponds to artifacts generated by aliased data rather than signal. The synthetic example given in this section provided an illustration of SPDR. However, the example in this section is limited to a model with a single zero-dip reflector which is, exactly, the least challenging scenario for SPDR. In the following two sections, we present more challenging examples for SPDR, where the reflectors have nonzero dip and the actual earth’s velocity is not constant.

WB129

0.05

0.00 0

10

20

30

40

50

60

Conjugate gradient iteration

SYNTHETIC DATA EXAMPLES FOR MODELS WITH VARIABLE DIP REFLECTORS

Figure 8. The least-squares normal equations are solved using leastsquares conjugate gradients. We show the data-misfit as a function of conjugate gradient iteration, showing the convergence of the algorithm.

Depth wavenumber (1/m)

Time (s)

Time (s)

Depth (m)

Thus far, we have limited ourselves to a synthetic example where the earth model consists of a single zero-dip reflector. The earth model at and above the reflector Lateral position (m) Lateral position (m) Lateral position (m) consisted of a single velocity of 1500 m / s, prea) 2000 4000 2000 4000 2000 4000 cisely what we use for c0 in the migration and 0 demigration operators. Moreover, the data set 1500 2000 1500 2000 1500 2000 m/s m/s m/s consisted of a single primary event, and neglect250 ed, for example, free-surface multiples. In this section, we apply SPDR to more involved syn500 thetic examples, where the models contain a 750 range of reflector dips. In a first example, we use a range of models, each containing a different dip b) 0.5 for its interface. This illustrates the dip limitation of the algorithm, and its connection to the Nyquist wavenumber. In a second example, we consider a model with multiple reflectors and a range of re1.0 flector dips 共to a maximum absolute dip of approximately ␰ ⳱ 0.7, or a slope of approximately 45° as measured from the horizontal兲.Additional1.5 ly, we allow the data to contain free-surface multiples, violating the single scattering approximac) 0.5 tion made in the migration/demigration operators, and used in our data reconstruction algorithm. 共The data also contain internal multiples, 1.0 etc.兲 For the first example, we illustrate the dip limitation of SPDR, and its relation to the lateral Ny1.5 quist wavenumber. Figure 9 shows the application of SPDR to three velocity models. From left to right, we show models with increasing reflecLateral wavenumber (1/m) Lateral wavenumber (1/m) Lateral wavenumber (1/m) tor dips. In particular, Figure 9a plots the velocity d) –0.1 0 0.1 –0.1 0 0.1 –0.1 0 0.1 models for these three cases; Figure 9b plots the 0.0 observed data built from the finite-difference approximation to the acoustic wave equation, and a –0.1 shot located at 共xs,zs兲 ⳱ 共2500 m,50 m兲; Figure 9c plots the SPDR reconstructed data; and Figure –0.2 9d plots the wavenumber spectra of the migrated gathers computed from the data in 9b. We notice that in the shallowest dip case, shown in the first –0.3 column of Figure 9, SPDR is successful. This success is due to the clear separation between the Figure 9. Synthetic data example showing the effect of dips: 共a兲 the velocity models with adjoint representation of data signal and alias dips 共angle to the horizontal兲 of, from left-to-right, 7.6°, 12.6°, and 17.3°, 共b兲 the obshown in the wavenumber spectra 共Figure 9d, served data, 共c兲 the SPDR reconstructed data, and d兲 the wavenumber spectra of the confirst column兲. On the other hand, in the second stant wavespeed migrated observed data. Each column plots a different experiment, with and third columns of Figure 9, the reflector dip is the dips of the reflectors increasing from left to right.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB130

Kaplan et al.

grid. In this example, the observed grid is regular with ⌬xog ⳱ 150 m. From the band-limited adjoint representation of data 共equation 26兲, we predict that the model-space artifacts are band-limited and centered at the lateral wavenumbers kgx ⳱ Ⳳ2␲ / 150⳱ Ⳳ0.04, and this is observed in Figure 9d. Therefore, increasing the Nyquist wavenumber by decreasing ⌬xog will, in turn, improve the separation between the model-space signal and artifacts 共regardless of reflector dip兲, and subsequently improve the quality of the SPDR result. This, of course, is not a surprising conclusion, and simply states that data reconstruction becomes less challenging as ⌬xog decreases. For the second example, Figures 10–15 show the experiment and the SPDR results. In Figure 10a, we plot the acoustic velocity and density models for the experiment. We use the finite-difference approximation to the acoustic wave equation to generate the synthetic shot gather in Figure 10b. The synthetic shot gather is built with a source located in the acoustic model at 共xs,zs兲 ⳱ 共2500 m, 50 m兲 and with a Ricker wavelet with a peak frequency of 15 Hz. The observations are on the regular grid xog such that geophones are placed every 100 m, and run from 0 m to 5000 m, each at a depth of 50 m. The observations are placed on the nominal grid xng with ⌬xng ⳱ 5 m spacing. In Figure 10b, the data are plotted using a subset of the nominal grid. In particular, we plot the recorded data for every fourth nominal grid sample and near offsets. This is for illustrative purposes only, and allows us to show in detail a representative portion of the data. Note that with this subset of data used, the plotted traces are

spaced every 20 m; hence, in our plots, we represent the data with four nominal grid points falling between each observed trace. This is applicable to the displays of the input data, SPDR reconstructed data, and the finite-difference constructed data 共the last two of which are shown in Figure 11兲. In Figure 10c, we plot the f ⳮ k spectrum of the recorded data. We solve equation 32 for the optimal model 共migrated shot gather兲 m*, subsequently computing the SPDR reconstructed data d* ⳱ Am . In finding m , we set the model weights such that ␰ˆ ⳱ 0, ␩ *

*

⳱ .09kN, ␶ ⳱ .01kN, and ␮ ⳱ 10ⳮ2, and where kN ⳱ 0.063 is the Nyquist wavenumber for the observed grid xog. The result, again shown for near offsets and every fourth trace on the nominal grid, is shown in Figure 11b. For comparison, we compute the data on the nominal grid using the finite-difference approximation to the acoustic wave equation, plotting the resulting data in Figure 11a. Figure 11c plots the difference between the SPDR reconstructed data in Figure 11b, and the finite-difference data in Figure 11a. Finally, in Figure 12, we plot the f ⳮ k spectra of the finite-difference and SPDR reconstructed data. In particular, Figure 12a plots the f ⳮ k spectrum of the finite-difference data, and Figure 12b plots the f ⳮ k spectrum of the SPDR reconstructed data. To illustrate the sensitivity of the method to the selection of ␩ , we ran SPDR several times for different values of ␩ , and plotted, in FigOffset (m)

a)

–400

–200

–400

–200

–400

–200

0

200

400

200

400

200

400

Lateral position (m) 1000

2000

3000

4000

5000

Depth (m)

0 250

1500

2000

ρ = 1000kg/m3

m/s

2500

3000

Time (s)

0.5

a)

1.0

500 ρ = 1500kg/m3

750

1.5

ρ = 2000kg/m3 ρ = 2500kg/m3

Offset (m)

b) Offset (m)

b)

–200

0

200

400

0.5

Time (s)

–400 0.5

Time (s)

0

1.0

1.0 1.5 1.5

Offset (m)

c)

0

Lateral wavenumber (1/m)

c)

–0.6

–0.4

–0.2

0

0.2

0.4

Time (s)

Frequency (Hz)

0.5

0.6

0 10

1.0

20 30

1.5

40

Figure 10. Synthetic data example: 共a兲 the acoustic velocity and density 共 ␳ 兲 model, 共b兲 The observed data shown for its near-offset traces, and using every fourth trace in the nominal grid, 共c兲 the f ⳮ k spectrum of the observed data.

Figure 11. Synthetic data example: 共a兲 data synthesized using the finite-difference approximation and the known velocity/density models, 共b兲 the SPDR reconstructed data, 共c兲 the difference between the finite-difference data in 共a兲 and the SPDR reconstructed data in 共b兲. In all plots, every fourth trace in the nominal grid is shown for the first 500 m of offsets.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR ure 13, the resulting signal-to-noise ratios 共S/N兲 as defined by

S/N⳱ 10 log10

储d1储22 储d1 ⳮ d*储22

共38兲

,

where d1 is the original data before decimation 共Figure 16a兲, and 储 · 储 2 is the L2 norm computed for an xg grid common to both d1 and d*. For this example, the method is sensitive to the choice of ␩ , resulting in the narrow width of the peak in Figure 13. This is due to the separation 共or lack thereof兲 between the model-space signal and artifacts. If the separation were larger due to either smaller reflector dips or finer sampling, then the method would become less sensitive to ␩ . In the extreme case where the aliased data is in the null space of the adjoint operator, SPDR is insensitive to ␩ , provided that it is set to a value that is larger than the bandwidth of the model-space signal. We use the example to consider the data reconstruction of freesurface multiple events. The single scattering approximation used in the demigration operator does not model free-surface multiples correctly. Indeed, a first-order approximation to a free-surface multiple requires a second-order scattering interaction 共e.g. Weglein et al., 2003兲. However, the single scattering approximation does see the multiple event as a false primary, and maps it to a copy of the true reflector at an erroneous depth, allowing for the extra traveltime of the wavefield through the water column. This can be seen in Figure 14. Lateral wavenumber (1/m)

a)

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

Frequency (Hz)

0 10

WB131

In Figure 14a, we allow the model space to have a maximum depth of 1000 m. This allows for the model space to contain the pseudodepth location of the reflectors, but does not allow for the first-order copy that is due to the free-surface multiples. Meanwhile, in Figure 14b, we allow the model space to extend to a depth of 2000 m. Then, the copy of the reflector that corresponds to the multiple energy falls within the bounds of the allotted model space. To clearly see the multiple energy, we applied a clip to the image. In turn, this shows various low-amplitude artifacts. The effect in the reconstructed data of enlarging the model space to allow for the multiple of the reflector is shown in Figure 15. In Figure 15a, we plot finite-difference data, and in Figure 15b and c, we plot the data reconstructions using, respectively, the larger and smaller model spaces. The corresponding dataspace effect of the extended model space is shown in Figure 15b, where the multiples are reconstructed. In contrast, the limited model-space result in Figure 15c excludes the multiples.

REAL DATA EXAMPLE For a real data example, we use a shot gather from a marine data set from the Gulf of Mexico. The shot gather has receivers spaced every 26.7 m, with the near-offset receiver placed at a distance of 20.7 m from the source and the far-offset receiver placed at a distance of 4874.7 m from the source. In total, there are 183 recorded geophones. In our experiment, we use a nominal grid with 365 points spaced every 13.3 m. Finally, we decimate the data so that the geophones are spaced every 80.0-m. In other words, there are five nominal grid points that need to be interpolated between each point on the observed grid.

20

Lateral position (m)

30

a)

200

400

600

800

1000

800

1000

0

Lateral wavenumber (1/m)

b)

–0.6

–0.4

–0.2

0

0.2

0.4

0.6

Frequency (Hz)

0

Depth (m)

40

500 750

10 20

Lateral position (m)

b)

30 40

400

600

0

500

Depth (m)

13 11 9 7 5 0

200

250

Figure 12. Synthetic data example: f ⳮ k spectrum of 共a兲 the finitedifference data in Figure 11a and 共b兲 the SPDR reconstructed data in Figure 11b. For comparison, consider the f ⳮ k spectrum of the observed data in Figure 10c.

SNR

250

750 1000 1250 1500 1750

0.05

0.1

0.15

0.2

0.25

0.3

Normalized lateral wavenumber

Figure 13. Synthetic data example: signal-to-noise ratios 共S/N兲 of the SPDR reconstructed data as a function of ␩ . The abscissa is normalized by the lateral Nyquist wavenumber.

Figure 14. Synthetic data example: 共a兲 The inverse m* computed for when the model space has a maximum depth of 1000 m; 共b兲 the inverse m* computed for when the model space has a maximum depth of 2000 m. The multiple energy in 共b兲 is indicated by the dashed box.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

Kaplan et al.

1.0

800

900

1000

1100

1200

1.8

–2.0 –1.5 –1.0 –0.5

1.8

2.8 3.3 3.8

–2.0 –1.5 –1.0 –0.5

1.8

2.3

2.3

2.8 3.3

Offset (m)

Time (s)

800

900

1000

1100

1200

a)

1.5 -

Offset (m)

Time (s)

800

900

1000

1100

1200

1.5

Frequency (Hz)

2.0

1.0

2.8 3.3 3.8

Lateral wavenumber (1/m) –0.01

0.01

b)

Lateral wavenumber (1/m) –0.01

Lateral wavenumber (1/m)

c)

0.01

0

0

0

10

10

10

20

20

20

30

30

30

40

50

Frequency (Hz)

Time (s)

-

2.0

c)

Offset (km) –2.0 –1.5 –1.0 –0.5

Figure 16. Real data example: 共a兲 the original data, 共b兲 the decimated data, 共c兲 the reconstructed data, 共d兲 the residual, or difference between 共c兲 and 共a兲. An automatic gain control algorithm is applied prior to plotting. In addition, we clip the data at 33% of the largest magnitude. All four plots use the same scale.

1.5

1.0

3.3

d)

3.8

b)

2.8

3.8

Offset (km) 1.8

–2.0 –1.5 –1.0 –0.5

2.3

Time (s)

2.3

c)

Offset (km)

b)

40

50

Frequency (Hz)

a)

Offset (km)

a)

Time (s)

Offset (m)

plots the adjoint of the decimated data 共from the 80.0 m observational grid兲. The artifacts resulting from the decimated data are evident. Finally, the third column 共Figure 18c and f兲 plots the weighted leastsquares inverse of the decimated data, which is used in producing the SPDR result in Figure 16c. In applying SPDR, we set ␩ ⳱ 0.09, ␶ ⳱ 0.01, and ␰ˆ ⳱ 0, the effect of which is seen in the comparison of Figure 18e and f. Additionally, we used 60 least-squares conjugate gradient iterations to find the solution of the least-squares normal

Time (s)

In Figure 16, we plot results from the SPDR algorithm. In Figure 16a, we plot the original real data on its original observational grid with 26.7 m spacing. In Figure 16b, we plot the data after decimation, and on its nominal grid with 13.3 m spacing. Finally, we show the SPDR result in Figure 16c. The result is plotted using the 26.7-m grid spacing of the original data. This allows us to compare the SPDR result to the observed data. In particular, Figure 16d plots the difference between the data in Figure 16a and the SPDR result in Figure 16c. Note that before plotting, we apply an automatic gain control algorithm to emphasize the small amplitude events. In addition, the gained data is clipped to 33% of the largest magnitude. It is clear from Figure 16d that SPDR has failed to reconstruct some portion of the back-scattered energy. This is likely because the energy falls within the null space of the demigration operator. In Figure 17, we plot the corresponding f ⳮ k spectra. In particular, Figure 17a plots the f ⳮ k spectrum of the original data on the 26.7-m grid, Figure 17b plots the f ⳮ k spectrum of the decimated data 共80.0-m grid兲, and Figure 17c plots the f ⳮ k spectrum of the SPDR reconstructed data. In Figure 18a-c, we plot various realizations from the model space, and in Figure 18d-f, we plot their wavenumber spectra. The first column 共Figure 18a and d兲 plots the adjoint of the original data 共on its 20.7-m grid兲, whereas the second column 共Figure 18b and e兲

Time (s)

WB132

–0.01

0.01

40

50

60

60

60

70

70

70

80

80

80

2.0

Figure 15. Synthetic data example, illustrating the data reconstruction of free-surface multiples: 共a兲 data generated for the nominal grid using the finite-difference approximation to the wave equation, 共b兲 the SPDR reconstructed data using a model space to 2000-m depth, 共c兲 the SPDR reconstructed data using a model space to 1000-m depth. The multiples are indicated by arrows in 共a兲 and 共b兲 and are absent from 共c兲.

Figure 17. Real data example: f ⳮ k spectra of 共a兲 the original data, 共b兲 the decimated data, and 共c兲 the reconstructed data.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR equations. The signal-to-noise ratio for the reconstructed data is 11 dB. For comparison, the signal-to-noise ratio for the first synthetic example, presented in Figure 7, is 16 dB.

DISCUSSION In this paper, the SPDR algorithm is limited to a single 2D shot gather. We imagine that the extension to a 3D shot gather would be straightforward. The extension to data reconstruction between shot gathers is more involved. However, we point out that the model space can be made independent of the acquisition geometry, including the shot location. Therefore, it is possible to reconstruct missing shot gathers using some extension of the SPDR algorithm. Whereas this paper explored the model-space separation of signal and adjoint that overlap in the data space, it did not give definite conditions for when the model-space signal and adjoint will be disjoint. Rather, we chose to make broader statements about the need for limited dip reflectors, giving a lower bound for the maximum dip. This presents a possible topic for future research: Find, given some acquisition geometry, rigorous bounds on the model space that would ensure the success of the algorithm. The forward and adjoint operators for the examples in this paper were implemented with the seismic source f共 ␻ 兲 in equations 4 and 11 set to 1. This assumption does not match the 2D source modeled in the synthetic data, or the 3D source inherent to the real data, and

–1

b)

0

Depth (km)

Depth (km)

2

Lateral wavenumber (1/m) –0.1

0

0.1

–0.6

2

e)

–2

–1

0

2

–0.1

0

0.1

–0.2

–0.4

–0.6

that is processed by a 2D algorithm. The least-squares component of SPDR makes this incomplete implementation sufficient. In short, the least-squares data fitting absorbs the effects of the source into the least-squares approximation of the scattering potential. As with many signal processing algorithms, SPDR includes a parameter selection problem. Although there are many methods for parameter selection in inverse theory, we do not extensively explore them in this paper. We note that for the examples presented, parameter selection was not onerous. We found that the most important parameter to select, not surprisingly, was ␩ , and a limited analysis of its effect was presented. It is used to determine the windowing function used to construct the model weights. In general, we chose it to be a value slightly smaller than the location of the alias in kgx, and which corresponds to k p for p ⳱ Ⳳ1, a property of the observed grid.

CONCLUSIONS We introduce a new data reconstruction algorithm, called shotprofile migration data reconstruction 共SPDR兲. The algorithm shares similarities with other wave-equation based data-reconstruction methods. It differs in its choice of model space for which we use a shot-profile migrated shot gather parameterized by lateral position and pseudodepth. An analysis of SPDR shows that it is applicable to earth models consisting of reflectors with limited dip. We illustrated the effectiveness of the algorithm using synthetic and real data examples. Whereas the single scattering approximation used does not explicitly model multiple reflections, a careful parameterization of the SPDR model space maps data multiples to false copies of the modelspace reflectors. In turn, this allows SPDR to interpolate multiple events, as well as primary events. On the other hand, portions of data corresponding to the evanescent wavefield are not reconstructed by the SPDR algorithm, namely because the evanescent wavefield is in the null space of the forward operator.

ACKNOWLEDGMENTS 3

f)

–0.1

0.0

Depth wavenumber (1/m)

Depth wavenumber (1/m)

–0.4

c)

0

Lateral wavenumber Lateral wavenumber (1/m) (1/m)

0.0

–0.2

–1

3

3

d)

–2

Depth (km)

–2

0

We thank the Signal Analysis and Imaging Group 共SAIG兲 consortium sponsors for their continued support.Additionally, we thank the associate editor and reviewers Adriana Citlali Ramírez and two anonymous reviewers for their careful reading, helpful suggestions, and criticisms.

0.1

0.0

Depth wavenumber (1/m)

a)

Offset (km)

Offset (km)

Offset (km)

WB133

–0.2

–0.4

–0.6

APPENDIX A DERIVATION OF THE FORWARD OPERATOR We find ␺ s in equation 4 using the Born approximation in equation 1, and the Green’s function in equation 2. Taking zg ⳱ zs ⳱ z0, and the support of ␣ to be below z0, we find by substitution of equation 2 into equation 1,

␺ s共xg,␻ ;xs兲 ⳱ Figure 18. Real data example: the model-space representation of 共a兲 the original data and 共b兲 the decimated data, both using the adjoint operator, and 共c兲 the decimated data under the weighted leastsquares inversion. The corresponding wavenumber spectra are shown in 共d兲 for the original data and 共e兲 for the decimated data, both using the adjoint operator, and 共f兲 the decimated data using the weighted least-squares inversion.

冉 冊冉 冊 冕冕 冕冉 1

3

2␲ ⬁

␻

2

c0

⬁

⬁

⫻

ⳮ

z0 ⳮ⬁ ⳮ⬁

1 i4kgz

冊

eⳮikgx·共x⬘ⳮxg兲

⫻eikgz共z⬘ⳮz0兲dkgx ⫻ ␣ 共x⬘,z⬘兲

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB134

Kaplan et al.

冕冉 冊

冉冊

⬁

⫻

ⳮ

ⳮ⬁

1

i4kz⬘

eikx⬘·x⬘eikz⬘共z⬘ⳮz0兲

In equation A-1, kx⬘ is the Fourier conjugate variable of x⬘, and we have recognized an expression for the synthetic source used in shotprofile migration algorithms,

g共kx⬘,␻ ,xs兲 ⳱ 2␲ f共␻ 兲eⳮikx⬘·xs .

␻2 ⳮ kx⬘ · kx⬘ . c20

共A-3兲

Rearranging terms and regrouping integrals in equation A-1 gives

␺ s共xg,␻ ;xs兲 ⳱

冉 冊冉 冊 1 2␲

3

␻ c0

冕冕

eikgx·xg

z0 ⳮ⬁ ⬁

⫻

冕

ⳮ⬁

u p共l兲共kgx,␻ 兲 ⳱ ⳮ

eikgz共zlⳮz0兲 , i4kgz

共B-2兲

and the earth model is partitioned into nz layers of constant thickness ⌬z. We implement equations B-1 and B-2 using two iterations. First, we define vs共l兲 for l ⳱ 1 . . . nz such that

vs共1兲共kgx,␻ ;xs兲 ⳱ u p共1兲共kgx,␻ 兲g共kgx,xs,␻ 兲 vs共l兲共kgx,␻ ;xs兲 ⳱ ⌬u p共kgx,␻ 兲vs共lⳮ1兲共kgx,␻ ;xs兲,

共B-3兲

冉

vr共1兲共kgx,␻ 兲 ⳱ u p共1兲共kgx,␻ 兲

冊

1 ⳮ eikgz共z⬘ⳮz0兲 i4kgz

vr共l兲共kgx,␻ 兲 ⳱ ⌬u p共kgx,␻ 兲vr共lⳮ1兲共kgx,␻ 兲.

共B-4兲

Then, the forward operator in equation B-1 becomes

冕 冉 ⬁

eⳮikgx·x⬘

where

where u p共1兲 ⳱ exp共ikgz共z1 ⳮ z0兲兲 / 共i4kgz兲 and ⌬u p ⳱ exp共ikgz⌬z兲. Second, we define vr共l兲 so that

2

⬁ ⬁

⫻

共B-1兲

共A-2兲

Moreover, we have used the dispersion relation

冑

⫻F 关F*u p共l兲共kgx,␻ 兲g共kgx,xs,␻ 兲兴␣ 共xg,zl兲,

共A-1兲

⫻g共k⬘x ,␻ ,xs兲dk⬘x dx⬘dz⬘ .

kz⬘ ⳱ sgn共␻ 兲

n

z ␻ 2 ␺ s共kgx,␻ ;xs兲 ⳱ ⌬z 兺 u p共l兲共kgx,␻ 兲 c0 l⳱1

eikx⬘·x⬘ ⳮ

ⳮ⬁

1 i4kz⬘

冊

␺ s共kgx,␻ ;xs兲 ⳱

n

z ␻ 2 ⌬z 兺 vr共l兲共kgx,␻ 兲 c0 l⳱1

⫻F 关F*vs共l兲共kgx,␻ ;xs兲兴␣ 共xg,zl;xs兲.

⫻eikz⬘共z⬘ⳮz0兲g共kx⬘,␻ ,xs兲dkx⬘ ⫻␣ 共x⬘,z⬘兲dx⬘dkgxdz⬘ .

冉冊

共A-4兲

Finally, we make a one-to-one change of notation for variables of integration, from kx⬘ to kgx and x⬘ to xg, and recognize three 2D Fourier transforms 共using the 1 / 共2␲ 兲 symmetric normalization factors兲 to give equation 4. Note that the change of integration variables to lateral geophone 共space and wavenumber兲 is somewhat arbitrary. In practice, we extend xg beyond the aperture of the geophones by adding zero traces to either side of the shot gather. Moreover, we can let xg correspond to the nominal grid used for data reconstruction, rather than a grid dictated by the geophone locations. This, in effect, makes the aperture of ␣ for shot xs 共the migrated shot gather兲 independent of the recording aperture.

APPENDIX B NUMERICAL IMPLEMENTATION OF THE FORWARD AND ADJOINT OPERATORS To produce efficient numerical implementations of the forward and adjoint operators in equations 4 and 11, we find iterative algorithms that implement the operators. The derivation of the iterations is straightforward, and the resulting algorithms should be familiar to those versed in shot-profile migration. We begin with the forward operator 共equation 4兲, rearranging terms and approximating the integral with a Riemann sum so that 共using a Fourier representation of the scattered wavefield兲

共B-5兲

Equations B-3–B-5 constitute an algorithm that implements the demigration operator in equation 4, and the bulk of computation is in the two 2D Fourier transforms required per depth and frequency. Next, we consider the implementation of the adjoint operator 共equation 11兲, approximating the integral over angular frequency with a Riemann sum, n␻

␣ †共xg,zl;xs兲 ⳱ ⌬␻ 兺

j⳱1

冉 冊 ␻j c0

2

关F*u*p共l兲共kgx,␻ 兲g*共kgx,xs,␻ 兲兴

⫻F*u*p共l兲共kgx,␻ 兲␺ s共kgx,␻ j;xs兲.

共B-6兲

Analogous to the forward operator, we implement equation B-6 us* for l ⳱ 1 . . . n so that ing two iterations. First, we define vs共l兲 z * 共k ,␻ ;x 兲 ⳱ u* 共k ,␻ 兲g*共k ,x ,␻ 兲 vs共1兲 gx s gx s p共1兲 gx * 共k ,␻ ;x 兲 ⳱ ⌬u*共k ,␻ 兲v* 共k ,␻ ;xs兲. vs共l兲 gx s p gx s共lⳮ1兲 gx

共B-7兲 * for l ⳱ 1 . . . n so that Second, we define vr共l兲 z

* 共k ,␻ ;x 兲 ⳱ u* 共k ,␻ 兲␺ 共k ,␻ ;x 兲 vr共1兲 gx s s gx s p共1兲 gx * 共k ,␻ ;x 兲 ⳱ ⌬u*共k ,␻ 兲v* 共k ,␻ ;xs兲. vr共l兲 gx s p gx r共lⳮ1兲 gx

共B-8兲 Then, the adjoint operator in equation B-6 becomes

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

SPDR

␣ †共xg,zl;xs兲 ⳱ ⌬␻ 兺 j

冉 冊 ␻j c0

2

Changing the order of integration and rearranging terms gives

* 共k ,␻ ;x 兲兴 关F*vs共l兲 gx s

* 共k ,␻ ;x 兲, ⫻F*vr共l兲 gx s

WB135

共B-9兲

␣ †p共kgx,z,␻ ;xs ⳱ 0兲 ⳱

1 1 f *共 ␻ 兲 ␺ 0共 ␻ 兲 2␲ ⌬xog ⬁

again requiring two 2D Fourier transforms per depth and frequency. Note that equation B-9 is shot-profile migration for a constant velocity reference medium.

⫻e

ⳮ⬁ ⬁

APPENDIX C ⫻e

DERIVATION OF ADJOINT SAMPLING OPERATOR

冤

⬁

1

冑2␲

冕

eikgxxgeⳮikgzzdkgx

ⳮ⬁

1

冑2␲

冕

eikgxxgeⳮikgzz

ⳮ⬁

⬘ Ⳮkp兲xgdx dk⬘ . eⳮikgxxgei共kgx g gx

To proceed, we recognize that

冥

⬘ Ⳮ k p兲兲其 Fⳮ1兵␦ 共kgx ⳮ 共kgx ⬁

⳱

1

1

冑2␲

␺ 0共 ␻ 兲 ⌬xog

⫻␦ 共kgx ⳮ k p兲dkgx,

冕

共C-6兲

⬁

⫻

ⳮikgz ⬘z

ⳮ⬁

We find ␣ †p in equation 20. We substitute equation 19 into equation 17 while letting xs ⳱ 0, so that

␣ †p共xg,z,␻ ;xs ⳱ 0兲 ⳱ 冑2␲ f *共␻ 兲

冕

ⳮikgz共p兲z

⳱

共C-1兲

where

冕

⬘ Ⳮ k p兲兲dkgx 共C-7兲 eikgxxg␦ 共kgx ⳮ 共kgx

ⳮ⬁

1

冑2␲ e

i共kgx ⬘ Ⳮkp兲xg

共C-8兲

.

Hence,

k p ⳱ 2␲ p/⌬xog .

共C-2兲

In equation C-1, we use the sifting property of the Dirac delta function so that

␣ †p共xg,z,␻ ;xs ⳱ 0兲 ⳱

1

1

冑2␲ ⌬xog f

*共 ␻ 兲 ␺ 0共 ␻ 兲

冕

eikgxxgeikpxgeⳮikgzzdkgx,

ⳮ⬁

共C-3兲 where

kgz共p兲 ⳱ sgn共␻ 兲冑␻ 2 /c20 ⳮ k2p .

共C-4兲

⬘ , letting kgz ⬘ Making a change of variables from kgx to kgx ⬘ 兲2, and taking the Fourier transform of equa⳱ sgn共 ␻ 兲冑␻ 2 / c20 ⳮ 共kgx tion C-3 with respect to xg, gives

1 1 ␣ †p共kgx,z,␻ ;xs ⳱ 0兲 ⳱ f *共 ␻ 兲 ␺ 0共 ␻ 兲 2␲ ⌬xog

1

冕

eⳮikgxxg

ⳮ⬁ ⬁

⫻

冕

ⳮ⬁

1

冑2␲ ⌬xog f

*共 ␻ 兲 ␺ 0共 ␻ 兲 ⬁

⫻eⳮikgz共p兲z

冕

⬘z eⳮikgz

ⳮ⬁

⬘ Ⳮ k p兲兲dkgx ⬘. ⫻␦ 共kgx ⳮ 共kgx

共C-10兲

As before, we use the sifting property of the Dirac delta function so ⬘ ⳱ kgx ⳮ k p in the integrand. In that the integral evaluates to where kgx particular, if we recall that

⬘ ⳱ sgn共␻ 兲 kgz

冑

␻2 ⬘, ⳮ kgx c20

共C-11兲

then, equation C-10 becomes

⬁

⫻eⳮikgz共p兲z

共C-9兲

Substituting equation C-9 into equation C-6 gives

␣ †p共kgx,z,␻ ;xs ⳱ 0兲 ⳱

⬁

⫻eⳮikgz共p兲z

⬘ Ⳮkp兲xg其 ⳱ 冑2␲ ␦ 共k ⳮ 共k⬘ Ⳮ k 兲兲. F兵ei共kgx gx p gx

␣ †p共kgx,z,␻ ;xs ⳱ 0兲 ⳱

1

1

冑2␲ ⌬xog f

*共␻ 兲␺ 0共␻ 兲eⳮikgz共p兲z

⫻eⳮi sgn共␻ 兲

冑␻2/c20ⳮ共kgx ⳮ kp兲2z . 共C-12兲

⬘ xgeikpxgeⳮikgz ⬘ zdk⬘ dx . 共C-5兲 eikgx gx g

Taking the Fourier transform of equation C-12 with respect to z gives

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/

WB136

Kaplan et al. ⬁

1 1 ␣ †p共kgx,kz,␻ ;xs ⳱ 0兲 ⳱ f *共 ␻ 兲 ␺ 0共 ␻ 兲 2␲ ⌬xog ⫻eⳮi共kgz共p兲Ⳮsgn共␻ 兲

冕

REFERENCES e

ⳮikzz

ⳮ⬁

冑␻2/c20ⳮ共kgx ⳮ kp兲2兲zdz. 共C-13兲

To simplify equation C-13, we let Fz denote the Fourier transform over z, and note that

Fzⳮ1兵␦ 共kz Ⳮ 共kgz共p兲 Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲其 ⳱

1

冑2␲

冕

eikzz␦ 共kz Ⳮ 共kgz共p兲

Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲dkz,

共C-14兲

which simplifies to

Fzⳮ1兵␦ 共kz Ⳮ 共kgz共p兲 Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲其 ⳱

1

冑2␲ e

冑

ⳮi共kgz共p兲Ⳮsgn共␻ 兲 ␻ 2/c20ⳮ共kgx ⳮ k p兲2兲z

.

共C-15兲

Hence,

Fz兵eⳮi共kgz共p兲Ⳮsgn共␻ 兲

冑␻2/c20ⳮ共kgx ⳮ kp兲2兲z其

⳱ 冑2␲ ␦ 共kz Ⳮ 共kgz共p兲 Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲. 共C-16兲 With equation C-16, we can simplify equation C-13, giving

␣ †p共kgx,kz,␻ ;xs ⳱ 0兲 ⳱

1

1

冑2␲ ⌬xog f

*共␻ 兲␺ 0共␻ 兲␦ 共kz Ⳮ 共kgz共p兲

Ⳮ sgn共␻ 兲冑␻ 2 /c20 ⳮ 共kgx ⳮ k p兲2兲兲, 共C-17兲 which is the same as equation 20.

Babu, P., and P. Stoica, 2010, Spectral analysis of nonuniformly sampled data — A review: Digital Signal Processing, 20, no. 2, 359–378. Baumstein, A., and M. T. Hadidi, 2006, 3D surface-related multiple elimination: Data reconstruction and application to field data: Geophysics, 71, no. 3, E25–E33. Biondi, B., 2003, Equivalence of source-receiver migration and shot-profile migration: Geophysics, 68, Short Note: 1340–1347. Chiu, S. K., and R. Stolt, 2002, Applications of 3-D data mapping — Azimuth moveout and acquisition-footprint reduction: 72nd Annual International Meeting, SEG, Expanded Abstracts, 2134–2137. DeSanto, J. A., 1992, Scalar wave theory: Green’s functions and applications: Springer-Verlag, volume 12, of Springer Series on Wave Phenomena. Gülünay, N., 2003, Seismic trace interpolation in the Fourier transform domain: Geophysics, 68, 355–369. Hansen, P. C., 1998, Rank-deficient and discrete ill-posed problems: Society for Industrial and Applied Mathematics. Herrmann, F. J., and G. Hennenfent, 2008, Nonparametric seismic data recovery with curvelet frames: Geophysical Journal International, 173, no. 1, 233–248. Huang, L.-J., M. C. Fehler, and R.-S. Wu, 1999, Extended local Born Fourier migration method: Geophysics, 64, 1524–1534. Lathi, B. P., 1998, Signal processing and linear systems: Berkely Cambridge Press. Liu, B., and M. D. Sacchi, 2004, Minimum weighted norm interpolation of seismic records: Geophysics, 69, 1560–1568. Miranda, F., 2004, Seismic data reconstruction of primaries and multiples: 74th Annual International Meeting, SEG, Expanded Abstracts, 2100– 2103. Naghizadeh, M., and M. D. Sacchi, 2007, Multistep autoregressive reconstruction of seismic records: Geophysics, 72, no. 6, V111–V118. Nemeth, T., C. Wu, and G. T. Schuster, 1999, Least-squares migration of incomplete reflection data: Geophysics, 64, 208–221. Paige, C. C., and M. A. Saunders, 1982, LSQR: an algorithm for sparse linear equations and sparse least squares: ACM Transactions on Mathematical Software, 8, no. 1, 43–71. Ramírez, A. C., and A. B. Weglein, 2009, Green’s theorem as a comprehensive framework for data reconstruction, regularization, wavefield separation, seismic interferometry, and wavelet estimation: A tutorial: Geophysics, 74, no. 6, W35–W62. Ramírez, A. C., A. B. Weglein, and K. Hokstad, 2006, Near offset data extrapolation: 76th Annual International Meeting, SEG, Expanded Abstracts, 2554–2558. Spitz, S., 1991, Seismic trace interpolation in the F-X domain: Geophysics, 56, 785–794. Stolt, R. H., 1978, Migration by Fourier transform: Geophysics, 43, 23–48. ——–, 2002, Seismic data mapping and reconstruction: Geophysics, 67, 890–908. Trad, D., 2003, Interpolation and multiple attenuation with migration operators: Geophysics, 68, 2043–2054. Trad, D., T. J. Ulrych, and M. Sacchi, 2003, Latest views of the sparse Radon transform: Geophysics, 68, 386–399. Weglein, A. B., F. V. Araújo, P. M. Carvalho, R. H. Stolt, K. H. Matson, R. T. Coates, D. Corrigan, D. J. Foster, S. A. Shaw, and H. Zhang, 2003, Inverse scattering series and seismic exploration: Inverse Problems, 19, no. 6, R27–R83.

Downloaded 03 Jan 2011 to 99.53.249.50. Redistribution subject to SEG license or copyright; see Terms of Use at http://segdl.org/