Migration methods for imaging different-order multiples

Geophysical Prospecting, 2007, 55, 1–19 Migration methods for imaging different-order multiples Zhiyong Jiang1∗ , Jianming Sheng1 , Jianhua Yu2 , Gerard T. Schuster1 and Brian E. Hornby2 1 University of Utah, Department of Geology and Geophysics, 135S 1460 E, Salt Lake City, UT 84112, USA; 2 BP America E&P Technology Group, 501 Westlake Park Boulevard, Houston, TX 77079, USA Received July 2005, revision accepted July 2006 ABSTRACT Multiples contain valuable information about the subsurface, and if properly migrated can provide a wider illumination of the subsurface compared to imaging with VSP primary reflections. In this paper we review three different methods for migrating multiples. The first method is model-based, and it is more sensitive to velocity errors than primary migration; the second method uses a semi-natural Green’s function for migrating multiples, where part of the traveltimes are computed from the velocity model, and part of the traveltimes (i.e., natural traveltimes) are picked from the data to construct the imaging condition for multiples; the third method uses cross-correlation of traces. The last two methods are preferred in the sense that they are significantly less sensitive to velocity errors and statics because they use “natural data” to construct part of the migration imaging conditions. Compared with the interferometric (i.e., crosscorrelation) imaging method the semi-natural Green’s function method is more computationally efficient and is sometimes less prone to migration artifacts. Numerical tests with 2-D and 3-D VSP data show that a wider subsurface coverage, higher-fold and more balanced illumination of the subsurface can be achieved with multiple migration compared with migration of primary reflections only. However, there can be strong interference from multiples with different orders or primaries when multiples of high order are migrated. One possible solution is to filter primaries and different orders of multiples before migration, and another possible solution is least squares migration of all events. A limitation of multiple migration is encountered for subsalt imaging. Here, the multiples must pass through the salt body more than twice, which amplifies the distortion of the image. INTRODUCTION In the typical migration strategy, only primary reflections are considered as signal and are therefore migrated. Other events, including multiples, are attenuated as much as possible before migration. However, multiples contain valuable information about the subsurface, and if properly migrated, they can provide wider illumination, higher fold, and better vertical res- ∗ E-mail: zjiang@mines.utah.edu  C 2007 European Association of Geoscientists & Engineers olution of the subsurface than primaries, as shown in Fig. 1 (Schuster 2003). In recent years, there has been some effort in the fields of multiple migration and various methods have been explored (Reiter et al. 1991; Berkhout and Vershuur 1994; Yu and Schuster 2001; Schuster 2003; Yu and Schuster 2004; Sheng 2001; Brown and Guitton 2005; Jiang et al. 2005). A strategy for multiple migration is to use the velocity model to compute the trajectory of multiple reflection rays and their associated traveltimes. The raypaths for the multiples are typically much longer than those for the primary reflections so the multiple migration image will be more sensitive to velocity 1 2 Z. Jiang et al. Figure 1 Multiples can provide wider illumination, higher fold, and better vertical resolution of the subsurface than primaries. The raypaths shown are for the multiple (solid) and primary (dashed) reflections. errors than the primary migration image. Moreover, the longer raypaths promote lower signal-to-noise ratios in the recorded multiples because of greater geometrical spreading and attenuation losses. How can these problems be overcome so we can benefit from the merits of multiple migration described above? To overcome the increased sensitivity to velocity errors, Claerbout (1968), and Rickett and Claerbout (1996, 1999) proposed a strategy for transforming “daylight” seismic data to CMP-like (common midpoint) data. In their formulation, the “daylight” data are assumed to be generated by a random distribution of buried sources with unknown excitation times. The strategy is to crosscorrelate a master trace with the surrounding traces to get a virtual CDP shot gather, and to sum these shot gather traces for every subsurface shot. The master trace position is also the virtual shot location. In this way the ghost reflections from the free surface are kinematically transformed to reflections for virtual sources at the surface; in general the Nth-order ghost is transformed into the N − 1th-order ghost. They demonstrated the feasibility of their approach with helioseismic data and synthetic seismic data. The important merit of this daylight imaging approach is that the velocity of the subsurface medium does not need to be known so that the migration of these virtual CMP data has the same sensitivity to velocity errors as the real CMP data. This greatly mitigates the velocity sensitivity problem that afflicts standard migration of ghost reflections. Later, Schuster and Rickett (2000) and Schuster (2001) extended the daylight imaging strategy to the more general procedure they denoted as interferometric imaging. Here, the restriction to random distributions of sources was eliminated, imaging of reflectors in arbitrary 3-D media was allowed and new imaging conditions were formulated to image source distribu-  C tions (Schuster et al. 2004), and interface geometry from transmitted arrivals (Sheng et al. 2003), converted arrivals (Sheng et al. 2003), VSP (Vertical Seismic Profile) multiples (Yu and Schuster 2001; Yu and Schuster 2002; Yu and Schuster 2004) and primaries and multiples in CDP (Common Depth Point) data (Sheng 2001; Zhou et al. 2005). For VSP data (Yu and Schuster 2004), the additional benefit is that the well statics are eliminated by interferometric imaging. An important development was the proof by Wapenaar (2003) and Wapenaar et al. (2004) that the daylight imaging concept introduced by Claerbout could be mathematically justified for 3-D media using Green’s theorem. A problem with interferometric imaging is that it can be computationally expensive. The traces must be crosscorrelated with one another and, since each trace acts as a secondary source, the interferometric migration can be more costly than standard migration (Yu and Schuster 2004). Moreover, crosscorrelation of traces can broaden the source wavelet and introduce spurious artifacts denoted as virtual multiples (Schuster et al. 2004). To mitigate these problems, Schuster (2002 and 2003) proposed the use of semi-natural Green’s functions to migrate seismic data. The idea is to construct the Kirchhoff imaging condition for multiples by a combination of traveltimes τ model computed from the model, and traveltimes τ̃ data picked from the data. Using Fermat’s principle, the calculation model data + τ̃xg ) yields the Kirchhoff imaging condiτsxg = min(τsx tion for multiples, where s is the source position, x is a trial image point, g is the geophone, and τ xx′ is the traveltime for waves to propagate from x to x′ . The advantage of this method over crosscorrelation migration (i.e., interferometric imaging) is that the traces do not need to be correlated, the computational expense is reduced, and migration artifacts are reduced. The disadvantages are that primary traveltimes of the strong multiple generator needs to be picked and the full Fresnel zone along the free-surface is not fully utilized to migrate the ghost energy; but wavepath migration can be used to mitigate this last problem (Sun and Schuster 1999). This paper provides the theory and numerical examples for migrating multiples of any order using standard modelbased migration as well as the semi-natural methods that use semi-natural Green’s functions (also designated as specular interferometric imaging) and interferometric imaging. For semi-natural methods the key requirement (and potential disadvantage) is that the traveltimes of a lower-order event must be picked (or windowed in the prestack data) in order to image a higher-order event. Based on the numerical results shown in this paper, we conclude that migration of multiples can sometimes provide increased structural resolution, higher fold 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 3 and wider subsurface illumination compared with migration of primary reflections only. However, this benefit will not be practical if the multiples are much weaker than the primaries. The first part of this paper presents the theory, and the next part presents numerical examples for both synthetic and field data. The final part summarizes the conclusions of this paper. T H E O RY There will be three parts to the theory section: model-based multiple migration, multiple migration by a semi-natural Green’s function and multiple imaging by interferometry. The latter two methods are semi-natural imaging methods in the sense that part of the imaging condition is obtained from the data and the other part calculated from the assumed velocity model. We will use the standard diffraction stack imaging formula in the frequency domain, but it will be understood that more sophisticated imaging formulas can be used as well such as Kirchhoff migration or wave equation migration. The example we will use below is mainly that of a 1st-order freesurface ghost, but any order free-surface ghost and pegleg multiples (provided they can be identified) are amenable to the migration methods described below. The free surface is denoted as the B 0 interface and the underlying multiple generating interface will be denoted as the B 1 surface (e.g., see Fig. 2). This multiple generating layer will often be referred to as the water layer in this paper. Model-Based Multiple Imaging We will first present the imaging formulas for migrating ghostrelated multiples in VSP data, and then for CDP data. These formulas use traveltimes computed entirely from the assumed velocity model. VSP Data. The general imaging formula in the frequency domain for a VSP ghost reflection (Fig. 2a) is Specular Refl. m(x, ω) =
g,s A(s, g, x)d(s, g)e  −iω(τsx +τxg ′ 0   + τg0′ g ) , (1) where d(s, g) represents a trace for a source at s and a receiver in the well at g; τ sx is the time for energy to propagate from the surface source to the trial image point at x; τ xg′0 is the time for energy to propagate from the trial image point at x to the free surface specular reflection point g′0 ; and τ g0′ g is the time energy takes to propagate back down to the VSP geophone at g from the free-surface reflection point at g′0 . Here  C Figure 2 Rays for (a) a source-side VSP ghost reflection and (b) a receiver-side CDP ghost reflection. In both cases xo is the subsurface specular reflection point, and the position go is the specular reflection point on the free surface corresponding to xo; x is the trial image point, and g′o is the specular reflection point corresponding to x. When x → x o we will have g′o → g o . g′ is a point on the free surface. See the theory section for details. A(s, g, x) is a function that determines the type of imaging method to be used; e.g., if A(s, g, x) consists of the Beylkin determinant (along with the obliquity factor) and the source and geometric spreading terms (Bleistein et al. 2001) then this is known as the Born asymptotic inversion formula. Without loss of generality, we will conveniently set A(s, g, x) = 1 so the imaging method described in this paper is diffraction stack migration. Coefficients such as powers of the frequency ω are absorbed into the d(s, g) function for notational tidiness. A ray is specular in the sense that the xg′0 g portion of the ray honors Snell’s law everywhere in the medium, including the bounce point g′0 at the free surface. The times τxg0′ and τ g0′ g without tildes are computed by a raytracing method using an assumed velocity model, unlike the semi-natural methods that use the picked traveltime τ̃g0′ g with a tilde. If the free surface is horizontal then a convenient way to compute the VSP traveltime for τxg0′ + τg0′ g is to mirror both the buried geophone at g 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 4 Z. Jiang et al. leg. Here x is the trial image point, x′B is the specular reflection point at the B 1 interface of the first leg and g′0 is the specular reflection point along the free surface B 0 for the trial image point x. If the reflector geometry is known, then ray tracing can be used to compute the specular primary reflection times ′ ′ τ (0) sg′ = τ sxB + τ xB g′ for all positions of the sources and receivers (otherwise known as a zeroth-order receiver-side ghost). This specular primary reflection time can be used with Fermat’s principle to compute the ghost migration imaging condition in equation 3:   (0) τsx′B + τx′B g0′ + τg0′ x + τxg = min (4) τ ′ + τg ′ x + τ xg . sg ′ g ǫ B0 Figure 3 The geometry for the mirror imaging method when the surface is flat. gmirror is the mirror location of g. Note that the mirror imaging method is a special case of the model-based multiple imaging method. and the subsurface velocity model across the free surface (Fig. 3). Then, trace a specular ray from the mirrored geophone at gmirror to the trial image point x. If the free surface along B 0 is irregular, then a simple application of Fermat’s principle will determine the specular τxg0′ + τg0′ g time, rather than using a mirror velocity model: (τxg′ + τg′ g ), τxg0′ + τg0′ g = min ′ (2) g ǫ B0 where g′ ranges over all the positions on the irregular free surface, and each specular reflection point location g′0 depends on the trial image point x for fixed source s and receiver g positions (Fig. 2). If there is more than one local minima then multi-arrival Kirchhoff migration should be used. The merit of this Fermat minimization approach is that it is simple and inexpensive to implement; but the problem is that the modelbased raypaths can be considerably longer than those for a primary reflection, hence the multiple migration image is more sensitive to migration velocity errors. CDP Data. For CDP data with 1st-order receiver-side ghosts (see Fig. 2b), the migration imaging formula is Specular Refl. m(x, ω) =
g,s d(s, g)e  −iω(τsx′ B   + τx′B g0′ +τg′ x +τxg ) 0 , Migration with Semi-Natural Green’s Functions To avoid some of the errors in calculating the specular reflection times in either equations 2 or 4 we replace some of the model-based traveltimes by their natural traveltimes (denoted by a tilde in τ̃xx′ ), where the natural traveltimes are picked from the data (Schuster 2002; Schuster 2003). That is, the VSP multiple migration equation takes a form similar to equation 1: Semi−natural Refl. m(x, ω) =
d(s, g)e g,s  −iω(τsx + τxg ′ 0   + τ̃g0′ g ) (5) except that the semi-natural reflection traveltime τxg0′ + τ̃g0′ g is partly obtained from the data (picked direct arrival times τ̃g0′ g ) and partly from the model (calculated τxg0′ times) using Fermat’s principle (Matsuoka and Ezaka 1992; Asakawa and Matsuoka 2002): (3) where τsx′B + τx′B g0′ is the specular traveltime for the first leg of the multiple and τg0′ x + τxg represents the traveltime of the last  C We notice that for given s and g locations, the position of g′o depends on the location of x. When x → x0 we will have g′0 → g0 , so that the left-hand side of equation 4 becomes τsxB + τxB g0 + τg0 x0 + τx0 g , which is the traveltime of the specular ray from s to g (see Fig. 2b). Again, the problem with equation 3 is that the reflection traveltimes must be computed from the assumed velocity model, which always has errors; and so the migration image accuracy will be more sensitive to the associated timing errors than the primary reflection image. To overcome this problem, we now introduce migration of multiples with semi-natural Green’s functions. τ xg0′ + τ̃ g0′ g Traveltime  Model Direct  ′′ = min τxg + ′′ g ǫ B0 Picked Direct Traveltime  τ̃g′′ g  . (6) 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 5 The natural direct traveltimes τ̃g′′ g are picked from the direct arrivals in the VSP data for all receiver positions g in the well and source positions g′′ on the surface. These natural traveltimes overcome two problems with the model-based migration in equation 1: they eliminate the well statics and partly reduce timing errors for a downgoing leg of the multiple raypath. The well statics are eliminated because the picked direct arrival traveltime in the VSP data contains, naturally, time shifts due to delayed source excitation times or statics in the well. In equation 5 the exponential is denoted as the semi-natural asymptotic Green’s function (neglecting the geometrical spreading term), because it is constructed from both the natural and model-based traveltimes. Imaging data with a semi-natural Green’s function will be denoted as specular interferometric imaging, and the first successful application of it to VSP field data was by Luo (2005). Similarly, the CDP multiple migration equation becomes m(x, ω) =
(0) d(s, g)e −iω(τ̃ ′ +τg ′ x +τxg ) sg 0 (7) 0 g,s (0) where τ̃sg ′ = τ̃sx′ + τ̃ x′ g ′ is the natural traveltime (picked from B B 0 0 the data) for the primary reflection from the B 1 interface for all source and receiver positions on B 0 . It is straightforward to include an obliquity term cos θ based on the ray angles from the source and specular reflection point at the surface. To use this imaging condition with wave equation migration we backward extrapolate the data using the, e.g., phase-shift migration method, and correlate these data with the shifted (0) iω(τ̃ ′ +τg ′ x ) sg 0 0 source wavelet W(ω)e , where W(ω) is spectrum of the source wavelet. The extension of this multiple migration to any order multiple reflection is easily extended to the Nth-order ghost by m(x, ω) =
(N−1) d(s, g)e −iω(τ̃ ′ sg 0 +τg ′ x +τxg ) 0 (8) , g,s (N−1) where τ̃sg is the natural traveltime for the N − 1th-order ′ 0 reflection in the water-layer, which can be picked from the data or estimated from lower-order water-layer reflections. For example, if the 0th-order multiple reflections in the water Figure 4 (a) VSP upgoing first-order multiple; (b) VSP downgoing second-order multiple. x 0 and g 0 are specular reflection points. layer are picked then the 1st-order water layer reflection can be estimated by Asakawa and Matsuoka (2002) = min ′′ g ǫ B0 m(x, ω) =
g,s d(s, g)e −iω( 0   (0) τ̃sg ′′  τ̃g(0) ′′ g + (9) . and these 1st-order and 0th-order natural reflections can be used to estimate any higher-order water layer reflection using Fermat’s principle. Similar to the above analysis, equation 5 can be extended to semi-naturally migrate any order of VSP multiples related to the free surface (e.g., see Fig. 4). The above migration formula migrates data to the last bounce of the multiple. In fact, the reflection energy can be naturally migrated to any of the bounce points (Schuster 2003). For example, if the multiple consists of both receiver-side and source-side ghosts then the obvious CDP migration formula is Natural Receiver−Side Primary Natural Source−Side Primary  (0) τ̃sg ′ 0th−order Refl. 0th−order Refl. (1) τ̃sg +τg ′ x +τxg ′′ + 0 0  τ̃g(0) ′′ g 0 ) (10) . (0th) (0) (0) where τ̃sg ′ + τg ′ x = ming′ ǫ B0 (τ̃sg ′ + τg ′ x ), and τ xg ′′ + τ̃g ′′ g = 0 0 0 ming′′ ǫ B0 (τxg′′ + τ̃g(0) ′′ g ).  C 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 0 6 Z. Jiang et al. sitivity to migration velocity errors (due to longer raypaths). To alleviate both of these problems we multiply d(g, g′ )direct∗ with d(g, s)ghost to get: If the migration formula above is summed over all frequencies then the diffraction stack migration formula in the spacetime domain is Natural Receiver−Side Primary Natural Source−Side Primary M(x) = m(x, ω)dω =
d(s, g)e iω( g,s  (0) τ̃sg ′  τ̃g(0) ′′ g +τg ′ x +τxg ′′ + 0 0 0 0 ) dω, (11) 
 (0) (0) d s, g, τ̃sg = ′ + τg ′ x + τ xg ′′ + τ̃g ′′ g , 0 0 0 0 g,s where d(s, g, t) represents the weighted trace at position g, for a source at s and recording time t. Migration of pegleg multiples can be accomplished by a straightforward application of the above principles. Migration of Multiples by Interferometric Imaging The theory of interferometric imaging of multiples (Schuster 2001) was previously described for CDP (Sheng 2001) and VSP multiples (Yu and Schuster 2001; Yu and Schuster 2004). Interferometric imaging of multiples differs from the previous approach using semi-natural Green’s functions in that the traces are crosscorrelated prior to migration, and reflections do not need to be picked from the data. These crosscorrelations shift the multiple reflections so that they become kinematically equivalent to primary reflections (Schuster 2003), which can then be migrated by a standard Kirchhoff imaging method. It will be shown that, asymptotically, the stationaryphase approximation to the interferometric imaging formula will yield migration kernels equivalent to those used for seminatural migration. We now present the interferometric principle for imaging downgoing first-order multiples (Fig. 2a). Assume that the downgoing direct wave arrival d(g, g′ )direct (for a surface source at g′ and VSP receiver at g) and ghost arrival d(g, s)ghost are recorded in the well and expressed as: d(g, g′ )direct = W(ω)eiω(τ̃gg′ +τ̃ static ) , d(g, s)ghost = m(xo, ω)W(ω)eiω(τ̃sxo +τ̃xo go +τ̃go g +τ̃ static ) (12) , where W(ω) is the spectrum of the source wavelet, τ̃xx′ is the natural traveltime for waves to propagate from x to x′ , τ̃ static is the receiver-related static error for a VSP geometry (note, τ̃ static is the source-related static error for an RVSP geometry), geometrical spreading is ignored, and m(x o , ω) is the reflection coefficient-like term at x o . The position g o is the specular reflection point on the free surface for the ghost reflection. Two problems with VSP imaging of ghosts that lead to smeared images are well receiver statics and increased sen-  C ∗ (g, g′ , s) = d(g, g′ )direct d(g, s)ghost = m(xo, ω)|W(ω)|2 eiω(τ̃sxo +τ̃xo go +τ̃go g −τ̃gg′ ) , (13) which becomes the correlated data in the time domain after a temporal Fourier transform. Note, the well receiver static errors τ̃ static are eliminated, and the traveltime τ̃go g for most of the downgoing ghost wave is largely eliminated by the direct wave time −τ̃gg′ if g′ is selected to be near g o .1 It is shown in the appendix that stationary-phase integration automatically selects g′ → g o . These correlated data can be migrated by diffraction-stack migration:
M(x) = m(x, ω), ω =



ω = g′ s


ω g′ (g, g′ , s)e−iω(τg′ x +τxs ) , g γ (s, g′ )e−iω(τg′ x +τxs ) , (14) s where for convenience we assume a wideband source wavelet so |W(ω)| = 1. Here x is the trial image point and γ (s, g′ ) = ′ g (g, g , s) represents the data extrapolated to the surface with the virtual geophone position at g′ and source at s (Jiang et al. 2005). In other words the VSP data have been kinematically transformed to be (approximately) virtual CDP data with the potential for much greater subsurface illumination compared to standard VSP data with primary reflections. Note, the τ ’s without tildes denote traveltimes computed by a ray tracing code. A stationary-phase justification that equation (14) will image the reflector at x o is given in the appendix. NUMERICAL TESTS Numerical tests are carried out on both synthetic and field data to illustrate the benefits and limitations of imaging multiples. 1 The statics due to the weathering zone along the free surface are still present. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 7 Figure 5 The dipping layer velocity model for the 2-D synthetic test. The eight interfaces are numbered 1 to 8, with 1 the shallowest interface and 8 the deepest. Six data sets are used: one set of 2-D VSP synthetic data, three VSP field data sets, one 3-D VSP synthetic data set and one 2-D CDP synthetic data example. 2-D Synthetic Data Test This VSP data set is simulated by computing a finite-difference approximation to the solution of the 2-D acoustic wave equation. The horizontal length of the model is 925 meters and the vertical distance is 1300 meters. There are nine layers with nine different velocities in the model (Fig. 5), and the minimum and maximum velocities are 1900 m/s and 4000 m/s, respectively. Ninety-two sources are evenly deployed along the free surface with a spacing of 10 meters. The first source is located at (10m, 0m) and the last source at (920m, 0m); the origin (0m, 0m) is located at the upper left corner of the model. One vertical well is located at the left boundary of the model. There are 91 receivers deployed in the well with an even spacing of 10 meters, and the position of the first receiver is (0m, 50m) and that of the deepest receiver is (0m, 950m). The wavefield recorded at the geophones in the well consists of upgoing and downgoing components. The two wavefields can be separated by F-K filtering, and then they can be migrated separately with different imaging conditions. Migration of First-order Multiples with Interferometric Imaging. This data set consists of 91 receiver gathers. But due to the expense of VSP experiments in practice, a small number of  C receivers are usually deployed in a well, especially in the case of offshore environments. Fig. 6(a) shows the primary migration image, and Fig. 6(b) shows the interferometric migration image of downgoing first-order multiples with only eight receiver gathers. These eight receivers are evenly deployed in the deep portion of the well, with depths ranging from 510 to 650 meters. In interferometric migration, the traces were shifted by the observed direct arrival time and migrated. This served as an inexpensive proxy for correlation. Figs 6(a) and 6(b) show that the multiples provide a wider subsurface coverage than given by the primaries. Also, the primaries are only able to illuminate below the receivers, while the multiples are able to see structures both above and below the receivers. Our test shows that even with only 8 receivers the multiples provide an image that is comparable to the primary image obtained with 91 receivers (which is not shown here). However, coherent artifacts are stronger in the multiple image, which include interference from higher-order multiples, interbed multiples, etc. Migration of First-order Multiples with Semi-natural Green’s Functions. The downgoing first-order multiples are also migrated with semi-natural Green’s functions (i.e., specular interferometry), and the migration image is shown in Fig. 6(c). There are no obvious differences between the images obtained with interferometry (Fig. 6b) and specular interferometry (Fig. 6c), except that the specular interferometric migration image shows slightly higher resulotion. The reason might be that in specular interferometry only the energy from the specular points along the free surface is imaged (see equation 5). In comparison, interferometric migration images the energy from all the scattering points on the free surface (see equation 14), which can have a smoothing effect in the final migration image. Entire Trace Crosscorrelation Migration. The theory was presented for migration of first-order multiples by interferometric imaging, where only the direct wave is used to crosscorrelate with the trace (see equation 13). In contrast, entire trace crosscorrelation migration correlates an entire trace with another trace (Yu and Schuster 2006). In theory, entire trace crosscorrelation migration has the potential advantage of migrating many different types/orders of multiples at the same time (as illustrated in Fig. 7), which can lead to better resolution. Fig. 6(b) shows the migration image obtained by interferometric migration of downgoing first-order multiples, and 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 8 Z. Jiang et al. 0 200 X (m) 400 600 800 0 200 200 400 400 600 600 Z (m) Z (m) 0 800 800 1000 1000 1200 1200 (a) Migration of Primary Reflections 0 200 X (m) 400 600 800 0 200 200 400 400 600 600 800 800 1000 1000 1200 1200 (c) Specular Interferometric Mig. of 1st-order Ghost 200 X (m) 400 600 800 (b) Interferometric Migration of 1st-order Ghost Z (m) Z (m) 0 0 0 200 X (m) 400 600 800 (d) Entire Trace Crosscorrelation Migratio n Figure 6 (a) Migration image of primary, (b) interferometric migration image of downgoing first-order multiples, (c) migration image of downgoing first-order multiples with semi-natural Green’s functions (i.e., specular interferometry), and (d) entire trace crosscorrelation migration image. All these images are obtained with data recorded by only 8 receivers. See the numerical tests section for details. Fig. 6(d) shows the entire trace crosscorrelation migration image. It can be seen that there is no significant difference between the two images. One reason might be that in VSP data, the downgoing first-order multiples are the dominant events compared to other types/orders of multiples, so that imaging of other types/orders of multiples is negligible. That is, the downgoing first-order multiples have the main contribution to the entire trace crosscorrelation migration image. To facilitate comparison, the migration images corresponding to a subset of the image are plotted in Fig. 8, which includes the primary migration image (Fig. 8a), multiple migration image with interferometric imaging (Fig. 8b), multiple migration image with semi-natural Green’s functions (Fig. 8c), and en-  C tire trace crosscorrelation migration image (Fig. 8d). As expected, the entire trace crosscorrelation migration image contains moderately more noise than the other images. This is largely due to coherent cross-talk noise. However, correlating traces with additive random noise has the desirable property of noise reduction. So there is a tradeoff in entire trace crosscorrelation migration between reduction of additive random noise and increased coherent noise from cross-talk. Statics. Besides the illumination range and coverage, interferometric migration of multiples has the advantage over primary migration of being insensitive to well statics. Fig. 9 shows an example when there are statics associated with receivers. It is 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 9 Figure 7 Ray diagrams for three different upgoing rays from a reverse VSP experiment. Correlating the entire trace at A with that at B yields a trace at B that is kinematically equivalent to a primary reflection with a source at A and receiver at B. This assumes that the trace at A only contains the above upgoing events and the trace at B only contains the free-surface reflections shown above. Figure 8 Similar to Fig. 6, but plotted in wiggles for a subset of the model to facilitate comparison. Migration images of (a) primary, (b) downgoing first-order multiples with interferometric imaging, and (c) downgoing first-order multiples with semi-natural Green’s functions (i.e., specular interferometry); (d) Entire trace crosscorrelation migration image.  C 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 10 Z. Jiang et al. 0 200 X (m) 400 600 800 0 200 200 400 400 600 600 Z (m) Z (m) 0 800 800 1000 1000 1200 1200 0 200 (a) 0 200 600 800 600 800 (b) X (m) 400 600 800 0 200 200 400 400 600 600 Z (m) Z (m) 0 X (m) 400 800 800 1000 1000 1200 1200 (c) 0 200 X (m) 400 (d) Figure 9 Migration images of the primary reflections in the presence of a random receiver static error, the maximum value of which is (a) 2.5 ms, (b) 6.25 ms, and (c) 12.5 ms. (d) Interferometric migration image of the downgoing first-order multiple. shown that when the static reaches a certain magnitude the primary image is affected (Fig. 9b,c), while the interferometric migration image is unaffected (Fig. 9d). Migration of Higher-Order Multiples. Fig. 10 shows the images obtained by migrating upgoing first-order multiples and downgoing second-order multiples (see Fig. 4) with the seminatural Green’s functions (see Equation 5). Fig. 10(a) shows the migration images obtained from upgoing first-order multiples reflecting from interface 8 in the model. Interface 1 refers to the shallowest interface and interface 8 is the deepest inter-  C face in our model. By saying the multiple is from some interface we mean that the bounce point A in Fig. 4 is located at that interface. Fig. 10(b) shows the migration images obtained from downgoing second-order multiples reflecting from interface 7 in the model. If we denote the primary as an upgoing 0th-order multiple, we can see that in both images there is interference from lower-order multiples and inter-bed multiples, and the images are consequently blurred. Thus the lower-order multiples should be filtered out before the higher-order multiples are migrated. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 11 200 X (m) 400 600 800 0 200 200 400 400 Depth (m) Depth (m) 0 0 600 200 X (m) 400 600 800 600 800 800 1000 1000 1200 1200 (a) Upgoing 1st Multiple, interface8 0 (b) Downgoing 2nd Ghost, interface7 Figure 10 (a) Migration of upgoing first-order multiples from interface 8. See Fig. 4(a) for the raypath. (b) Migration of downgoing second-order multiples from interface 7. For the raypath refer to Fig. 4(b). The solid dark lines show the locations of the interfaces in the actual model. The images are blurred due to interference from low-order multiples and inter-bed multiples. 2-D IVSPWD Data An Inverse Vertical Seismic Profiling While Drilling (IVSPWD) data set donated by Union Pacific Resources Corporation (UPRC) is processed with the crosscorrelation migration algorithm. The UPRC data are recorded on the earth’s free surface while a tri-cone drill bit and down-hole motor were used to drill along a horizontal trajectory at the depth of around 2800 m in the Austin Chalk formation. The schematic of the acquisition survey is shown in Fig. 11. The traces in any “shot” gather have a recording length of 20 seconds with a sample interval of 2 ms. Each “shot” record is acquired by ten 3-component receivers on the surface with offsets ranging from about 821 m to 1915 m relative to the drill-rig position. The drill bit and the receiver line approximately lie in the same vertical plane. The detailed description of the UPRC data and the processing procedure are given in Yu and Schuster (2001 and 2002). Fig. 12 shows the crosscorrelograms with correlation window lengths of 8 s, 12 s and 16 s. The window overlap length is 10% of the correlation window length and the master trace is selected to be at the location of the 6th trace. A prestack time-migration algorithm with the ghost imaging condition is applied to the crosscorrelograms with an antialiasing migration operator (Lumley et al. 1994) and a wavelet deconvolution. The final surface-related ghost images by ghost crosscorrelogram migration with different window lengths are  C Figure 11 The map view of the IVSPWD acquisition survey. The circles denotes the receivers deployed on the surface. Each shot gather is recorded by ten 3-component receivers on the surface with offsets ranging from about 821 m to 1915 m relative to the drill-rig position. The drill bit is along the depth of 2800 m while recording these data. Drill bit and receiver line approximately lay in the same vertical plane. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 12 Z. Jiang et al. Traces 5 Traces 5 10 2 3 2 3 4 2 3 4 (a) 10 1 Time (s) 1 Time (s) Time (s) 1 Traces 5 10 4 (b) (c) Figure 12 Crosscorrelograms with different correlation window lengths: (a) 8s, (b) 12s, and (c) 16s. The window overlap length is 10 percent of the correlation window length. The master trace is at position 6. obtained. The final migration image is displayed with a trace interval of 3 m. Fig. 13 shows the comparisons of the crosscorrelogram migration images (inserted) with the near-by surfaceCDP section of Line AC4. The reflectors in the crosscorrelogram migration image roughly correlate with the main reflectors in the surface-CDP section, but some parts of the ghost image are polluted by extra events in the crosscorrelogram migration results. This may be caused by the unintended imaging of primary reflections, virtual multiples or the fact that surface CDP line is about a half-mile away from the location of the IVSPWD image. This example illustrates one of the key advantages of crosscorrelation migration. Namely, if first arrivals can not be easily identified (such as in the IVSPWD data example) then semi-natural migration is not a straightforward possibility.  C 2-D VSP Marine Data This VSP marine data set was recorded using 12 receivers between the depths of 4568 meters and 4736 meters. Each trace has a recording length of 14 seconds with a sample interval of 2 ms. The velocity model is shown in Fig. 14, where the positions of the well and the receiver array are also shown. The first-order downgoing multiple (Fig. 2a) is migrated using three methods: model-based imaging, migration with semi-natural Green’s functions (i.e., specular interferometry), and interferometric imaging. In the interferometric imaging, the traces were shifted by the observed direct arrival time and migrated. This served as a proxy for correlation. The corresponding migration images are shown in Figs 15(a) to 15(c). It can be seen that there is little difference between the images of 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 13 0 0 0.2 Distance (m) 0.4 0.6 x10 4 0.8 1.0 0.2 Depth (m) 0.4 0.6 0.8 1.0 x10 4 2000 Figure 13 Comparison of ghost crosscorrelogram image with surface stacked section of Line AC4. The correlation window length is 12 s, and the gray solid arrow indicates an approximate projection of the drill-bit trajectory on Line AC4. model-based imaging and semi-natural Green’s functions (i.e., specular interferometry), which means that the velocity errors and the receiver statics are small. Comparing the images of interferometry and specular interferometry, the former shows lower frequencies in the image. The reason might be that the energy from all the scattering points along the free surface is used in interferometric imaging (see equation 14), which has a smoothing effect in the final migration image. In comparison, specular interferometry only images the energy from the specular point on the free surface (see equation 5). To facilitate comparison, the migration images obtained with the three methods are plotted in wiggles for a subset of the model, which are shown in Fig. 16(a) to 16(c). 3-D VSP Marine Data A 3D VSP marine data set was recorded using 12 receivers between the depths of 4.69 km and 4.86 km. Data along one radial line of shots were used to migrate the ghost reflections to give the image at the top of Fig. 17. A section obtained by migrating primary reflections in the CDP data collected over  C 3000 Vel (m/s) 4000 Figure 14 The velocity model for the 2-D VSP marine example. The thin vertical line at the distance of 5.4 km shows the position of the well, and the thick vertical bar at the depth of about 4.6 km denotes the receiver array. the same area is shown at the bottom of Fig. 17. Both images show a similar reflectivity distribution, but the VSP image has about 1.5 times the temporal resolution of the CDP data. The ghost image was computed by ghost migration with a mirror imaging condition rather than interferometric migration. Nevertheless, this image compared well with the CDP image and appeared to be of higher frequency than the CDP data. The receivers were mostly below the image area so the reflectivity shown in Fig. 17 is unseen by traditional migration of VSP primary reflections. 3-D Synthetic VSP Data Test A 3-D VSP data set was generated by solving the acoustic wave-equation using a finite-difference method. The dimension of the 3-D velocity model used is 2 km by 2 km by 2 km. Altogether 1089 sources are evenly deployed on the surface, with 111 receivers evenly distributed between the depths of 100 meters and 1200 meters along a vertical well, whose position is in the center of the model. The recording length for the seismic traces is 6.0 seconds, with a temporal sampling interval of 0.001 second. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 14 Z. Jiang et al. Distance (m) 4000 5000 6000 7000 3000 2000 3000 3000 4000 4000 Depth (m) Depth (m) 3000 2000 5000 6000 7000 5000 6000 6000 7000 7000 8000 Distance (m) 4000 5000 8000 (a) Model-Based Migration 3000 Distance (m) 4000 5000 6000 (b) Migration with Specular Interferometry 7000 2000 3000 Depth (m) 4000 5000 6000 7000 8000 (c) Interferometric Imaging Figure 15 Multiple migration images with (a) model-based imaging (the mirror imaging condition is used here), (b) semi-natural Green’s functions (i.e., specular interferometry), and (c) interferometric imaging. Similar to Fig. 14, the thin vertical line and the thick vertical bar denote the positions of the well and the receiver array, respectively. Primaries and downgoing first-order multiples are migrated from this 3-D data set. The mirror imaging method is used to migrate the downgoing first-order multiples because the surface is perfectly flat in this case. The direct waves are muted, and upgoing and downgoing wavefields have been separated using f-k filtering before migration. Migration results for slices in the Y and X directions are shown in Figs 18 and 19. From Figs 18 and 19, we can see that the salt boundaries are better imaged with the downgoing first-order multiples than with the primary reflections. The primary reflections do not clearly illuminate the structures above the salt, which are clearly imaged with the multiple reflections. This example illustrates a greater subsurface coverage enjoyed by ghost imaging compared to imaging of primaries in VSP data.  C 2-D CDP Data Test Synthetic 2-D acoustic data are generated for the surface sources and receivers in the SEG/EAGE salt model. The freesurface multiples were migrated by the crosscorrelation migration method (Sheng 2001; Schuster 2001). Now, we test the effectiveness of migration ghost-reflections with the seminatural Green’s function (Schuster 2002; Schuster 2003). The model is modified by increasing the depth to the ocean bottom by 366 meters (see Fig. 20) such that the water-layer primary reflections can be easily identified and picked. A 4th-order finite-difference solution to the 2-D acoustic wave equation (Levander 1988), with perfectly matched layer (PML) boundary conditions (Festa and Nielsen 2003), is used to compute 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 15 3500 3500 4500 4500 5000 Distance (m) 5400 5600 4000 Depth (m) 4000 Depth (m) 4000 Distance (m) 5400 5600 Depth (m) 3500 Distance (m) 5400 5600 4500 5000 (a) 5000 (b) (c) Figure 16 Similar to Fig. 15, but plotted in wiggles for a subset of the model to facilitate comparison. Multiple migration images with (a) model-based imaging, (b) semi-natural Green’s functions (i.e., specular interferometry), and (c) interferometric imaging. Figure 17 Migration images obtained by migrating (top) downgoing first-order multiple in 3D VSP marine data and (bottom) primary reflections in 3D CDP data over the same area. The VSP image has higher resolution than the CDP image. of the primary reflections in which some artifacts (indicated by the white arrows) due to the water-bottom multiples appear at about 2 km at the left side, and at 1 and 3 km to the right side. Migration with the semi-natural Green’s functions (equation 11) for both source- and receiver-side ghosts gives the image shown in Fig. 21(b), in which ghost reflections are migrated to their correct position while primary reflections are not correctly migrated and result in the artifacts. In Figs 21(a) and 21(b), it can be seen that the reflectors are correctly imaged and the distribution and pattern of the artifacts in the two images are quite different. This means that, at the actual reflector positions, the two images are quite consistent while elsewhere they typically disagree. Thus, the artifacts can be suppressed by weighting the primary migration image with the local coherence (coherency window is about 10 traces wide and several wavelengths tall) between the primary and ghost migration images (Sheng et al. 2003), and the merged image is shown in Fig. 21(c) which shows fewer artifacts compared to the primary image. CONCLUSIONS the seismograms with a 8 Hz peak frequency Ricker wavelet as the source wavelet. A total of 325 shots (54 m shot spacing) and 176 traces per shot (27 m trace spacing) are migrated and Fig. 21(a) shows the Kirchhoff depth migration image  C Three methods for migrating different-order multiples in the VSP and CDP cases are presented. The first method is to use the velocity model to predict the traveltimes of the multiples. Due to the longer raypaths of the multiples, the multiple 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 16 Z. Jiang et al. Figure 18 (a) A vertical slice of the velocity model at X=1000 meters for the 3-D synthetic VSP example. The vertical line shows the location of the well; the circles denote the receivers in the well; (b) Primary migration image for the same slice as shown in (a); (c) Downgoing first-order multiple migration image for the same slice. The white lines in the migration images show the locations of the reflectors in the velocity model. migration image obtained with this method will be more sensitive to the velocity errors than the primary migration image. The second method of multiple migration uses semi-natural Green’s functions and Fermat’s principle to compute the traveltimes for low-order specular reflections. Practical applications of CDP multiple migration with semi-natural Green’s functions will be mostly limited to CDP data containing reflections from strongly reflecting multiple generators, e.g., sea  C Figure 19 (a) A vertical slice of the velocity model at Y=1000 meters. The vertical line shows the location of the well; the circles denote the receivers in the well; (b) Primary migration image for the same slice as shown in (a); (c) Downgoing first-order multiple migration image for the same slice. The white lines in the migration images show the locations of the reflectors in the velocity model. floor or salt bodies. In these cases the primary reflections can be identified in the common shot gathers and windowed, and the crosscorrelation of these windowed events can be used to estimate τ̃sg0′ (e.g., see equation 7). This procedure can be automated without time picking, just as time shifts are estimated for estimating residual static shifts from common shot gathers. These estimated traveltimes can then be used to estimate the specular traveltimes by Fermat’s principle for use in multiple migration with semi-natural Green’s functions. For 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 17 0 3 Distance (km) 6 9 12 0 15 3 Distance (km) 6 9 12 15 Depth (km) 1 Depth (km) 1 2 2 3 3 (a) Kirchhoff Depth Migration (Primary) SEG/EAGE Salt Model Figure 20 The 2-D SEG/EAGE salt model. 0  C Distance (km) 6 9 12 15 Depth (km) 1 2 3 (b) Semi-Natural Depth Migration (Ghost) 0 3 Distance (km) 6 9 12 15 1 Depth (km) VSP data, the direct arrivals can be estimated by a similar procedure. The third method images multiples by migrating crosscorrelograms. With the latter two methods multiple migration will have the same sensitivity to the velocity errors as primary migration, and in the VSP case, the well statics are eliminated. Compared with interferometric imaging, migration with a semi-natural Green’s function reduces the computational cost and migration artifacts because now the traces do not need to be crosscorrelated before migration. However, the key advantage of interferometric imaging over migration with semi-natural Green’s functions is when first arrivals can not be identified (such as in the IVSPWD data example) then semi-natural migration is not possible, while interferometric imaging is still applicable. Another potential advantage is the reduction of additive random noise by correlation of traces. Table 1 shows a comparison among the multiple and primary imaging methods. In VSP case, if there are ns shots on the surface and ng receivers in the well, for each common shot gather there are ng traces. In interferometric imaging the well receivers are redatumed to be at surface shot positions (see equation 14), and now each shot gather has ns traces. The increase in trace number is on the order of ns/ng, which is also approximately the increase in cost for interferometric imaging compared with primary migration (Table 1), if we neglect the cost in crosscorrelating traces. Note here we assume a dense distribution of sources on the free surface. If the source distribution is sparse, we will need to interpolate the ns sources to ns′ sources; thus the increase in cost for interferometric imaging compared with primary migration is on the order of ns′ /ng, instead of ns/ng. Migration results with VSP multiples show that, for our examples, a wider subsurface coverage and more balanced illumination can be achieved by migrating downgoing first- 3 2 3 (c) Kirchhoff Depth Migration (Primary+Ghost) Figure 21 (a) The Kirchhoff depth migration image for the primary reflections. The white arrows indicate the artifacts due to the waterbottom caused ghost-reflections. (b) The Kirchhoff depth migration image for the source- and receiver-side ghost reflections using seminatural Green’s function. (c) The coherence-weighted migration image obtained by weighting the primary image with the local coherence between the primary and ghost migration images. order ghosts. Even wider subsurface coverage and greater fold from the many bounce points can be achieved by migrating upgoing first-order multiples and downgoing secondorder ghosts. However, the problem with migration of upgoing 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 18 Z. Jiang et al. Table 1 Imaging and computational properties of migration methods: primary migration, model-based multiple migration, migration with semi-natural Green’s function and interferometric imaging. (ns is the number of shots on the surface and ng is the number of receivers in the well for VSP case.) Migration methods CPU cost CPU cost Sensitivity Receiver statics Receiver geometry Applicable Coverage for 2-D w.r.t. for 3-D w.r.t. to velocity in VSP case needs to be known to in VSP primary mig. primary mig. errors eliminated? in VSP case? IVSPWD? case Model-based multiple migration (eq. 1) Semi-natural Green’s functions (eq. 5) Interferometric imaging (eq. 14) Primary migration ≃1 ≃1 ≃ns/ng 1 ≃1 ≃1 ≃ns/ng 1 first-order multiples (Fig. 4a) and downgoing second-order ghosts (Fig. 4b) is that the interference from the primary and downgoing first-order ghosts (Fig. 2a) seriously blurs the images. One possible solution is to filter different order multiples and then migrate a specific order of multiple. Other possible solutions include least squares migration of all events (Brown and Guitton 2005), least squares migration filtering (Nemeth et al. 2000; Yu and Wang 1999) or deconvolution (Calvert et al. 2004). A limitation with migration of the downgoing first-order multiple (ghost) below salt is defocusing by a double passage through the salt when a large salt body exists in the model. A possible solution is to migrate multiples from salt boundaries. In the case of VSP first-order ghost migration in marine environment, for a deep water bottom (two-way traveltime greater than 1 second) the first-order ghosts are often well separated from primary reflections and higher-order free-surface related multiples; thus, there is not a strong need to eliminate the primary reflections or higher-order multiples prior to first-order ghost migration. On the other hand, these events should be eliminated if there is a shallow water bottom (less than approximately 0.5 second two-way traveltime in the water column). ACKNOWLEDGEMENTS We thank the members of the University of Utah Tomography and Modeling/Migration (UTAM) Consortium (http://utam.gg.utah.edu) for their financial support of this research. We also thank the Center for High Performance Computing at University of Utah for the use of their facilities. Discussions with Yi Luo were very helpful, especially in coining the term semi-natural migration. Panos Kelamis and Helmut Jacubowicz were instrumental in suggesting the importance of using natural data for redatuming and superresolution.  C high low low low No Yes Yes No Yes No No Yes No No Yes No wide wide wide narrow REFERENCES Asakawa E. and Matsuoka T. 2002. Traveltime-based raytracing for multiples and converted waves. Extended Abstracts of Annual EAGE Meeting (Florence, Italy), P102. Berkhout A.J. and Vershuur D.J. 1994. Multiple technology: Part 2, migration of multiple reflections. 64th Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts, 1497–1500. Bleistein N., Cohen J.K. and Stockwell J.W. Jr. 2001. Mathematics of Multi-dimensional Seismic Imaging, Migration, and Inversion. Springer Verlag, New York. Brown M. and Guitton A. 2005. Least-squares joint imaging of multiples and primaries. Geophysics 70, S79–S89. Calvert R.W., Bakulin A. and Jones T.C. 2004. Virtual sources, a new way to remove overburden problems. 66th Annual Conference, EAGE, Expanded Abstracts, P234. Claerbout J.F. 1968. Synthesis of a layered medium from its acoustic transmission response. Geophysics 33, 264–269. Festa G. and Nielsen S. 2003. PML absorbing boundaries. Bull. Seism. Soc. Am. 93, 891–903. Jiang Z., Yu J., Schuster G.T. and Hornby B.E. 2005. Migration of multiples. The Leading Edge 24, 315–318. Levander A. 1988. Fourth-order finite-difference P-SV seismograms. Geophysics 53, 1425–1437. Lumley D.E., Claerbout J.F., and Bevc 1994. Anti-aliasing Kirchhoff 3D migration. 64th Ann. Internat. Mtg., Soc. Expl. Geophys., Expanded Abstracts, 1282–1285. Luo Y. 2005. Single fold image based on specular interferometry theory. submitted to The Leading Edge. Matsuoka T. and Ezaka T. 1992. Raytracing using reciprocity. Geophysics 57, 326–333. Nemeth T., Sun H. and Schuster G.T. 2000. Separation of signal and coherent noise by migration filtering. Geophysics 65, 574– 583. Reiter E.C., Toksoz M.N., Keho T.H. and Purdy, G.M. 1991. Imaging with deep-water multiples. Geophysics 56, 1081–1086. Rickett J. and Claerbout J. 1996. Passive seismic imaging applied to synthetic data. Stanford Exploration Project 92, 83–90. Rickett J. and Claerbout J. 1999. Acoustic daylight imaging via spectral factorization: Helioseismology and reservoir monitoring. The Leading Edge 18, 957–960. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19 Migration methods for imaging different-order multiples 19 Schuster G. and Rickett J. 2000. Daylight imaging in V(x,y,z) media. Utah Tomography and Modeling-Migration Project Midyear Report, and Stanford Exploration Project Midyear Report. Schuster G.T. 2001. Theory of Daylight - Interferometric Imaging: Tutorial. Extended Abstracts of Annual EAGE Meeting (Amsterdam, Holland). Schuster G.T. 2002. Joint migration of primaries and multiples by a semi-natural Green’s functions. Utah Tomography and Modeling/Migration Consortium Annual Report, 141–164. Schuster G.T. 2003. Imaging the most bounce out of multiples. Extended Abstracts of Annual EAGE Meeting (special session on multiples). Schuster G.T., Yu J., Sheng J. and Rickett J. 2004. Interferometric/daylight seismic imaging. Geophysical Journal International 157, 838–852. Sheng J. 2001. Migration of multiples and primaries in CDP data by crosscorrelogram migration. SEG Expanded Abstracts, 1297–1300. Sheng J., Schuster G.T., Pankow K.L., Pechmann J.C. and Nowack R.L. 2003. Coherence-weighted wavepath migration of teleseismic data. Eos Trans. AGU 84(46), Fall Meet. Suppl., Abstract S11E0344. Sun H. and Schuster G.T. 1999. Wavepath migration versus Kirchhoff migration. SEG Expanded Abstracts, 1138–1141. Verschuur D.H., Berkhout A.J. and Wapenaar C.P.A. 1992. Adaptive surface related multiple elimination. Geophysics 57, 1166– 1177. Wapenaar C.P.A. 2003. Synthesis of an inhomogeneous medium from its acoustic transmission response. Geophysics 68, 1756–1759. Wapenaar C.P.A., Thorbecke J., and Draganov D. 2004. Relations between reflection and transmission responses of three-dimensional inhomogeneous media. Geophysical Journal International 156, 179–194. Yu J. and Schuster G.T. 2001. Crosscorrelogram migration of IVSPWD data. 71st Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts, 456–459. Yu J. and Schuster G.T. 2002. Joint migration of primary and multiple reflections in RVSP data. 72nd Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts, 2373–2376. Yu J. and Schuster G.T., 2004. Illumination of the subsurface by crosscorrelogram migration. 74th Ann. Internat. Mtg, Soc. Expl. Geophys., Expanded Abstracts. Yu J. and Schuster G.T. 2006. Crosscorrelogram migration of inverse vertical seismic profile data. Geophysics 71, S1–S11. Yu J. and Wang Y. 1999. 4C Mahogony data processing and imaging with least squares migration filtering. Utah Tomography and Modeling/Migration Consortium Annual Report, 373–398.  C Zhou M., Jiang Z., Yu J. and Schuster G.T. 2005. Comparison between interferometric migration and reduced-time migration of CDP data. submitted to Geophysics. APPENDIX By plugging equations (13) into 14 we will get:


M(x) = ω g s Specular Ghost Reflection ×
m(xo, ω)e  iω((τ̃sxo g′ Diffraction      + τ̃xo go + τ̃go g ) − (τ̃gg′ + τg′ x + τxs )) , (A1) where the diffraction time is always greater than or equal to the specular reflection time when x → x o . Stationary-phase integration says that the dominant contribution to the g′ summation is when the exponential argument is stationary, that is when x → x o and the diffraction point g′ coincides with g o . In this case the exponential argument in equation (A1) becomes zero so summation over s, g, and ω will be in phase and give perfect migration focusing at x o . In other words, the stationary-phase approximation (for large ω) to equation (A1) is


M(x) = ω g s Specular Ghost Refl. × m(xo, ω)
g′ Semi−naturalDiffraction       iω((τ̃sxo + τ̃xo go + τ̃go g ) − (τ̃gg ′ + τg ′ x + τxs )) e , Semi−natural Green′ s Fcn. Ghost Refl.      


−iω(τ̃gg ′ +τg ′ x +τxs ) iω(τ̃sxo +τ̃xo go +τ̃go g ) 0 0 e , ≈ Cm(xo, ω) e ω g s (A2) where C is an asymptotic coefficient and the far right exponential is identical to the semi-natural Green’s function in equation (5) for semi-natural migration of VSP ghosts. As mentioned before, when x → x o the argument of the exponential becomes zero, leading to coherent summation over the s, g and ω indices. 2007 European Association of Geoscientists & Engineers, Geophysical Prospecting, 55, 1–19