Skip Nav Destination
Filter
Filter
Filter
Filter
Filter

Update search

Filter

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

- Title
- Author
- Author Affiliations
- Full Text
- Abstract
- Keyword
- DOI
- ISBN
- EISBN
- ISSN
- EISSN
- Issue
- Volume
- References
- Paper Number

### NARROW

Format

Subjects

Date

Availability

1-15 of 15

Yaxun Tang

Close
**Follow your search**

Access your saved searches in your account

Would you like to receive an alert when new items match your search?

*Close Modal*

Sort by

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the SEG International Exposition and Annual Meeting, October 11–16, 2020

Paper Number: SEG-2020-3411534

Abstract

Strong surface and guided waves in land seismic data often distort signals from deeper structures, leading to poor S/N. Accurate modeling and attenuation of aliased surface waves has been a long-standing challenge in onshore seismic processing and imaging. To address this issue, we propose a method based on elastic full waveform inversion (FWI) using the spectral element seismic simulator to model and attenuate surface waves from the recorded data. This method first reconstructs surface waves by estimating a relatively shallow earth model via elastic multi-parameter FWI. The reconstructed surface waves are then used through an adaptive subtraction process to remove the surface waves from the recorded seismograms. The proposed method does not impose any 1-D assumption as conventional dispersionbased modeling approach does, and hence can model complicated wave phenomena, including back scattered surface waves. We demonstrate the performance of the proposed method using the SEAM II Arid model. Presentation Date: Tuesday, October 13, 2020 Session Start Time: 1:50 PM Presentation Time: 3:05 PM Location: Poster Station 13 Presentation Type: Poster

Proceedings Papers

Musa Maharramov, Anatoly I. Baumstein, Yaxun Tang, Partha S. Routh, Sunwoong Lee, Spyros K. Lazaratos

Publisher: Society of Exploration Geophysicists

Paper presented at the 2017 SEG International Exposition and Annual Meeting, September 24–29, 2017

Paper Number: SEG-2017-17744416

Abstract

ABSTRACT We propose a new full-waveform inversion (FWI) method that approximates broadband tomographic inversion and has a reduced sensitivity to errors in the observed-data amplitude information. The method is based on fitting the observed-data phase spectrum while automatically shaping forward-modeled wave fields to the observed-data amplitude spectrum. This is achieved by using a phase-only objective function that allows broadband time-domain inversion of the observed-data phase information. We demonstrate our method’s reduced sensitivity to dynamic information under the traveltime approximation, and compare the new objective function to the normalized FWI objective function in an experiment on synthetic data with a frequency-dependent attenuation. Presentation Date: Monday, September 25, 2017 Start Time: 3:30 PM Location: 361F Presentation Type: ORAL

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2015 SEG Annual Meeting, October 18–23, 2015

Paper Number: SEG-2015-5918876

Abstract

Summary We present an efficient second-order optimization technique for multi-parameter full wavefield inversion (FWI). We ap- proximate the Hessian matrix that contains important curvature information of the objective function by a simple sampling and interpolation approach, taking advantage of its sparse structure. We show that the approximate Hessian is small enough so that it can be stored in either computer disk or memory, and interpolated on the fly when applied. Computing the Hessian preconditioned gradient, therefore, becomes cost effective. We also show that the curvature information brought in by the approximate Hessian properly scales/rotates the gradients for parameter classes with different sensitivities to the data, so that simultaneous updating multiple parameter classes becomes possible. We test our method using both 2-D synthetic and 3-D field data sets for both single-parameter FWI and multi-parameter vertical transverse isotropy (VTI) FWI. Introduction Full wavefield inversion (FWI) is a nonlinear inversion technique that tries to recover the earth model by minimizing the mismatch between simulated and observed seismic wavefields (Tarantola, 1984). Due to its huge computational cost, FWI is often solved using local optimization techniques. A widely used local optimization technique is the gradient-based first- order approach, such as the steepest descent and nonlinear conjugate gradient (Nocedal and Wright, 1999). The first-order approach is relatively efficient, because it requires computing only the gradient of the objective function, but its convergence is usually slow. The convergence can be significantly improved by using the second-order optimization technique, which uses not only the gradient information, but also the curvature information of the objective function. The main difference between the first- and the second-order approach is that the latter preconditions the gradient with the inverse of the Hessian, such as in the Gauss- Newton/Newton approach (Pratt et al., 1998; Tang and Lee, 2010; Schiemenz et al., 2014), or with the inverse of a pro- jected Hessian, such as in the subspace approach (Kennett et al., 1988; Baumstein, 2014). The Hessian is a matrix containing second-order derivatives of the objective function with respect to the model parameters. The second-order approach is attractive not only because of its fast convergence rate, but also because of its capability to properly scale/rotate gradients of different parameter classes (e.g. velocity, anisotropy, attenuation, etc.) that are of different sensitivity to the data, hence pro- viding meaningful updates for all parameters simultaneously. Computing the inverse of the Hessian or the Hessian itself or even the product of the Hessian and a vector, however, is ex- pensive, and it is the main obstacle that prevents the second- order approach from being widely used in practice.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2013 SEG Annual Meeting, September 22–27, 2013

Paper Number: SEG-2013-1145

Abstract

SUMMARY We present a method based on gradient separation and recombination to improve the long wavelength updates of full wavefield inversion (FWI) on reflection-dominant data. We decompose the FWI gradient into a tomographic and a migration term. The tomographic term is obtained by crosscorrelating wavefields traveling in similar direction (forward scatterings), whereas the migration term is obtained by crosscorrelating wavefields traveling in opposite directions (backward scatterings). We show that the tomographic component of the gradient mainly updates the long wavelengths of the model parameters, whereas the migration component of the gradient mainly updates the short wavelengths of the model parameters. We recombine these two components to form a new gradient with a higher weight on the tomographic component. We show that the tomographically enhanced gradient makes FWI converge to the correct model even in situations where the conventional FWI fails to converge. We demonstrate our method based on constant-density acoustic FWI by inverting the P-wave velocities, but the method can be easily extended to more complicated multi-parameter inversion.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2012 SEG Annual Meeting, November 4–9, 2012

Paper Number: SEG-2012-0332

Abstract

SUMMARY We propose a new method to perform wave-equation migration velocity analysis by maximizing the flatness of the angle-domain common image gathers. Instead of maximizing the image-stack-power objective function directly with respect to the slowness, we link the objective function to the slowness indirectly through an intermediate moveout parameter. This approach is immune to the cycle-skipping problem, and it produces high-quality gradients. In addition, the proposed method does not require explicit picking of the moveout parameters. Our numerical examples demonstrate the great potential of this method: in the first example where there is a Gaussian-shaped anomaly slowness error, our method produces well-behaved gradient; from our test on the Marmousi models, the proposed method converges to a high-quality model that uniformly flattens the angle-domain common image gathers.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2011 SEG Annual Meeting, September 18–23, 2011

Paper Number: SEG-2011-3958

Abstract

ABSTRACT We apply target-oriented wavefield tomography to a 3-D field data set acquired from the Gulf of Mexico. Instead of using the original surface-recorded data set, we use a new data set synthesized specifically for velocity analysis to update subsalt velocities. The new data set is generated based on an initial unfocused target image and by a novel application of 3-D generalized Born wavefield modeling, which correctly preserves velocity kinematics by modeling non-zero subsurface-offset-domain images. We show that the target-oriented inversion strategy drastically reduces the data size and the computation domain for 3-D wavefield tomography, greatly improving its efficiency and flexibility. We apply differential semblance optimization (DSO) using the synthesized new data set to optimize subsalt velocities. The updated velocity model significantly improves the continuity of subsalt reflectors and yields flattened angle-domain common-image gathers.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2010 SEG Annual Meeting, October 17–22, 2010

Paper Number: SEG-2010-4077

Abstract

Summary Wave-equation migration velocity analysis (WEMVA) is an image-domain tomography based on a wave-equation propagator (one-way or two-way). Different algorithms for residual calculation, propagation and optimization add varying flavors to WEMVA. The two most popular ways to compute the residual are either using differential semblance optimization (DSO) or using differential residual migration (DRM). In this paper, we present relevant theory to understand the difference these two algorithms make in the final gradient calculation. We compare and contrast the two methods with the help of numerical examples. Both methods have some advantages and disadvantages over the other and they should be kept in mind while making the choice. The theory and results presented are based on one-way wave-equation migration.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2010 SEG Annual Meeting, October 17–22, 2010

Paper Number: SEG-2010-4280

Abstract

SUMMARY We present a method to reduce the computational cost of image-domain wavefield tomography. Instead of using the originally recorded data for velocity estimation, the proposed method simulates a new data set obtained using Born modeling or demigration based on the initial image and gathers. The modeling can be performed in a target-oriented fashion, and it can use arbitrary types of source functions and acquisition geometries. Hence the size of the new data set can be substantially smaller than the original one. We demonstrate with numerical examples that the new data set correctly preserves velocity information useful for velocity estimation, and that it generates wavefield-tomography gradient similar to that obtained using the original data set. We apply the proposed method to a modified version of the Sigsbee2A model, where two square anomalies below the salt have been successfully recovered in a target-oriented fashion at much lower computational cost.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2010 SEG Annual Meeting, October 17–22, 2010

Paper Number: SEG-2010-1034

Abstract

SUMMARY Full waveform inversion (FWI) has received an increasing amount of attention thanks to its ability to provide a highresolution velocity model of the subsurface. The computational cost still presents a challenge, however, and the convergence rate of the FWI problem is usually very slow without proper preconditioning of the gradient. While preconditioners based on the Gauss-Newton Hessian matrix can provide significant improvements in the convergence of FWI, computation of the Hessian matrix itself has been considered highly impractical due to its cost in computational time and storage requirements. In this paper, we design preconditioners based on an approximate Gauss-Newton Hessian matrix obtained using the phase-encoding method. The new method requires only 2 Ns forward simulations compared to Ns ( Nr + 1) forward simulations required in conventional approaches, where Ns and Nr are the numbers of sources and receivers, respectively. We apply the diagonal of the phase-encoded Gauss- Newton Hessian to both sequential-source FWI and encoded simultaneous-source FWI. Numerical examples using a truncated Marmousi2 model demonstrate that the phase-encoded Gauss-Newton Hessian improves the convergence of the FWI significantly. INTRODUCTION Velocity estimation is a challenging task in seismic exploration. In the past decade, ray-based velocity estimation methods have been widely used in practice to build the velocity model. Although the ray-based method is efficient and robust, the infinitefrequency assumption and inherent caustics in ray theory prevent it from accurately handling complex geologies. On the other hand, full-waveform inversion (FWI), which uses bandlimited wavefields instead of infinite-frequency rays as carriers of information, can properly model the finite-frequency effects of seismic wavefields and is immune from multi-pathing issues. Moreover, FWI uses full waveforms in the observed data for velocity estimation. Therefore, FWI has the potential to provide a detailed velocity model even in complex geological scenarios. The theoretical advantages, as well as recent advances in seismic acquisition technologies and computer sciences, have made FWI an attractive tool for accurate and highresolution velocity model estimation. FWI is usually implemented within the framework of nonlinear least-squares inversion, where a data residual defined either in time or frequency domain is minimized in the leastsquares sense (Tarantola, 1984; Pratt, 1999). Since this is a large-scale nonlinear optimization problem, gradient-based local optimization methods are often used in practice. Although the gradient of the least-squares misfit functional can be efficiently calculated using the adjoint state method, a gradientonly method, such as steepest descent and conjugate gradient, converges slowly without proper preconditioning. On the other hand, the Gauss-Newton method, which preconditions the gradient with a Gauss-Newton Hessian, can achieve significantly faster convergence (Pratt et al., 1998). Computing the exact Gauss-Newton Hessian, however, is prohibitively expensive for practical applications due to the size of the problem. In this paper, we focus on the constant-density acoustic FWI problem and discuss a method based on random-phase encoding to efficiently compute the Gauss-Newton Hessian (Tang, 2009). With this new method, the cost for computing the Hessian diagonal is only two forward modeling steps for each shot.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2009 SEG Annual Meeting, October 25–30, 2009

Paper Number: SEG-2009-3914

Abstract

ABSTRACT We present a joint least-squares inversion method for imaging simultaneous source (or blended) time-lapse seismic data sets. Non-repeatable shot and receiver positions introduce undesirable artifacts into time-lapse seismic images. We conjecture that more artifacts will result from relative shot-timing non-repeatability, when data sets are acquired with several simultaneously shooting sources. We show that these artifacts can be attenuated by joint inversion of such data sets without need for initial separation. Preconditioning with non-stationary dip filters and with temporal smoothness constraints ensures stability and geologically consistent time-lapse images. Results from a modified Marmousi 2D model show that this method yields reliable time-lapse images.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2009 SEG Annual Meeting, October 25–30, 2009

Paper Number: SEG-2009-2859

Abstract

SUMMARY We present a method based on wave-equation least-squares migration/ inversion to directly image data collected from recently developed wide-azimuth acquisition geometry, such as simultaneous shooting and continuous shooting, where two or more shot records are often blended together. We show that by using least-squares migration/inversion, we not only enhance the resolution of the image, but more importantly, we also suppress the crosstalk or acquisition footprint, without any pre-separation of the blended data. We demonstrate the concept and methodology in 2-D and apply the data-space inversion scheme to the Marmousi model, where an optimally reconstructed image, free from crosstalk artifacts, is obtained. INTRODUCTION High quality seismic images are extremely important for subsalt exploration, but data collected from conventional narrow-azimuth towed streamer (NATS) often produce poor subsalt images due to insufficient azimuth coverage. Recently developed wide-azimuth towed streamer (WATS) (Michell et al., 2006) and multi-azimuth towed streamer (MATS) (Keggin et al., 2006; Howard and Moldoveanu, 2006) acquisition technologies have greatly improved the subsalt illumination and hence better subsalt images are obtained. However, acquiring WATS or MATS data is expensive. One main reason is the inefficiency of the conventional way of acquiring data, which allows sufficient time sharing between shots to prevent interference (Beasley et al., 1998; Beasley, 2008; Berkhout, 2008). As a consequence, the source domain is often poorly sampled to save the survey time. To gain efficiency, simultaneous shooting (Beasley et al., 1998; Beasley, 2008; Hampson et al., 2008) and continuous shooting, or more generally, blended acquisition geometry (Berkhout, 2008), have been proposed to replace the conventional shooting strategy. In the blended acquisition geometry, we try to shoot and record continuously, and consequently, the time sharing between shots would be minimized and a denser source sampling can be obtained. However, this shooting and recording strategy results in two or more shot records blending together and brings processing challenges. A common practice of processing these blended data sets is to first separate the blended shot gathers into individual ones (Spitz et al., 2008; Akerberg et al., 2008), as called ”deblending” by Berkhout (2008). Then conventional processing flows are applied to these deblended shot gathers. The main issue with this strategy is that it might be extremely difficult to separate the blended gathers when the shot spacing is close and many shots are blended together. In this paper, we present an alternative method of processing these blended data sets. Instead of deblending the data prior to the imaging step, we propose to directly image them without any pre-separation. The simplest way for direct imaging would be migration, however, migration of blended data generates images contaminated by crosstalk. The crosstalk is due to the introduction of the blending operator (Berkhout, 2008), which makes the corresponding combined Born modeling operator far from unitary, thus its adjoint, also known as migration, gives poor reconstruction of the reflectivity. A possible solution is to go beyond migration by formulating the imaging problem as a least-squares migration/inversion (LSI) problem, which uses the pseudo inverse of the combined Born modeling operator to reconstruct the reflectivity of the subsurface.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2009 SEG Annual Meeting, October 25–30, 2009

Paper Number: SEG-2009-3964

Abstract

ABSTRACT We present the image-space wave-equation tomography to the generalized source domain, where a smaller number of synthesized shot gathers are generated. Specifically, we generate synthesized shot gathers by image-space phase encoding. The comparison of the gradients of the tomography objective functional obtained using image-space encoded gathers with that obtained using the original shot gathers shows that those encoded shot gathers can be used in wave-equation tomography problems. The advantages of using image-space encoded data is the decreased computational when compared to that of the original shot gathers and the ability of performing velocity analysis in a target-oriented way. We illustrate our examples on Marmousi model.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2008 SEG Annual Meeting, November 9–14, 2008

Paper Number: SEG-2008-2201

Abstract

SUMMARY I demonstrate a method for computing wave-equation Hessian operators, also known as resolution functions or point-spread functions, under the Born approximation. The proposed method modifies the original explicit Hessian formula, enabling efficient computation of the operator. A particular advantage of this method is that it reduces or eliminates storage of Green''s functions on the hard disk. The modifications, however, also introduce undesired crosstalk artifacts. I introduce two different phase-encoding schemes, namely, plane-wave phase encoding and random phase encoding, to suppress the cross-talk. I applied the Hessian operator obtained by using random phase encoding to the Sigsbee2A synthetic data set, where a better subsalt image with higher resolution is obtained. INTRODUCTION Migration is an important tool for imaging subsurface structures using reflection seismic data. The classic imaging principle (Claerbout, 1971) for shot-based migration states that reflectors are located where the forward-propagated source wavefield correlates with the backward-propagated receiver wavefield. However, it has been found that such an imaging principle is only the adjoint of the forward Born modeling operator (Lailly, 1983), which provides reliable structural information of the subsurface but distorts the amplitude of the reflectors because of the non-unitary nature of the Born modeling operator. To improve relative amplitude behavior, the imaging problem can be formulated as an inverse problem based on the minimization of a least-squares functional. The inverse problem can be formulated either in the data space (Lailly, 1983; Tarantola, 1984) or in the model space (Beylkin, 1985; Chavent and Plessix, 1999; Plessix and Mulder, 2004; Valenciano et al., 2006; Yu et al., 2006). The data-space approach can be solved iteratively using the gradient-based method (Nemeth et al., 1999; Clapp, 2005) without explicit construction of the Hessian, the matrix of the second derivatives of the error functional with respect to the model parameters. The iterative solving, however, is relatively costly and converges very slowly. On the other hand, the model-space approach requires the explicit construction of the Hessian, and its pseudo-inverse is applied to the migrated image. The full Hessian is too big and expensive to be computed in practice; hence Chavent and Plessix (1999); Plessix and Mulder (2004) approximate it by a diagonal matrix. In the case of high-frequency asymptotics, and with an infinite aperture, the Hessian is diagonal in most cases (Beylkin, 1985). For a finite range of frequencies and limited acquisition geometry, however, the Hessian is no longer diagonal and not even diagonally dominated (Plessix and Mulder, 2004; Valenciano et al., 2006). It has been shown by Valenciano et al. (2006) that, in areas of poor illumination, e.g., subsalt regions, the Hessian''s main diagonal energy is smeared along its off-diagonals. Therefore, the migrated image pre-multiplied by a diagonal matrix cannot perfectly recover the amplitude information, especially in poorly illuminated areas. That''s why Valenciano et al. (2006) suggest computing a limited number of the Hessian off-diagonals to compensate for poor illumination and improve the inversion result. However, computing the Hessian off-diagonals, even for a limited number, is very expensive by directly implementing the explicit Hessian formula.

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2007 SEG Annual Meeting, September 23–28, 2007

Paper Number: SEG-2007-2320

Abstract

SUMMARY I analytically demonstrate the artifacts in angle-domain common image gathers caused by sparsely sampled wavefields, from a perspective of shot-profile migration. The local-offset gather is linearly related to the angle gather in locally constant-velocity media when wavefields are sufficiently well sampled, but not when wavefields are poorly sampled. Hence, linear slant-stack or radial-trace transform in local-offset gathers will produce angle gathers with artifacts, which may hinder further interpretation or analysis and reduce the quality of the final stacking image in the angle and azimuth domain. Instead of simply stacking along reflection angle and azimuth axes, I present a method to compute the stacking weights as functions of angle and azimuth and make the stacking process selective. My method is tested on the synthetic wide-azimuth version of the SEG/EAGE salt data set, where a cleaner image with higher signal-to-noise ratio is obtained. INTRODUCTION Seismic image quality is highly dependent on the acquisition geometry, or more specifically, on the illumination of the subsurface. Ideally, the best image is obtained when each subsurface image point is illuminated equally, which potentially requires an infinite recording geometry. In the real world, however, we always have a limited recording geometry, like in the marine-streamer acquisition system, where only narrow-azimuth data are acquired. Recent developments in multi-azimuth or wide-azimuth acquistion techniques (Michell et al., 2006; Keggin et al., 2006; Howard and Moldoveanu, 2006) provide us richer coverage in azimuth, and the subsurface can be better illuminated, especially in subsalt areas. Much better images are obtained because of improved subsurface illumination. Nevertheless, the reality can never meet the requirement for infinite recording geometry and infinitly dense sampling for shots and receivers. Shadow zones or illumination holes and aliases may still happen in complex geologies. Poor illumination results in poor images, and related artifacts destort the migrated image, making it difficult to interpret. This effect is readily visible in the reflection-angle and azimuth domain, where illumination holes and related artifacts can be identified. They are by no means random or weak, and thus simple stacking can not attenuate them effectively. In this paper, I briefly review methods for extracting subsurface offset-domain common-image gathers (SODCIGs) and angle-domain common-image gathers (ADCIGs). I demonstrate the effects of sparsely sampled wavefields on both SODCIGs and ADCIGs. I also demonstrate that the final image formed by simple stacking over the reflection-angle and azimuth axes without any weighting function suffers from artifacts caused by poor illumination and has a low signal-to-noise ratio. Instead, I describe a simple but effective way to make the stacking process selective: we stack only those re- flection angles and azimuths with good illumination. This method is tested on the wide-azimuth version of the SEG/EAGE salt data set, and a better image with higher signal-to-noise ratio is obtained. EXTRACTING ANGLE-DOMAIN COMMON-IMAGE GATHERS ADCIGs can be extracted either before applying an imaging condition (Prucha and Symes, 1999; Mosher and Foster, 2000; Xie and Wu, 2002; Soubaras, 2003) or afterward (Sava and Fomel, 2003; Biondi and Symes, 2004).

Proceedings Papers

Publisher: Society of Exploration Geophysicists

Paper presented at the 2006 SEG Annual Meeting, October 1–6, 2006

Paper Number: SEG-2006-0189

Abstract

SUMMARY We propose a method for selecting reference-anisotropy parameters in laterally varying anisotropic media for mixed Fourier-space domain wave field extrapolation. We treat the selection problem as a quantization procedure, and use a modified version of the 3D Lloyd''s algoritm for reference-parameter selections. We demonstrate that our method yields a more accurate description of the anisotropy model with fewer reference parameters than the uniform sampling approach. Synthetic and real data examples illustrate the performance of our method. INTRODUCTION It is well known that wave field-continuation-based migration methods are better able to image areas affected by multi pathing and are more effective in handling complex wave behavior than ray tracing based Kirchhoff methods. However, wave-equation depth migration such as PSPI (Gazdag and Sguazzero, 1984), extended SSF (Stoffa et al., 1990; Kessinger, 1992), split-step DSR (Popovici, 1996), or FFD plus interpolation (Biondi, 2002) require multiple reference velocities for wave field extrapolation through laterally inhomogeneous velocity models. With multiple reference velocities higher angles can be accurately handled and the quality of the image can be improved, but the cost of the migration is increased in proportion to the number of reference velocities used. When it comes to anisotropic media, the vertical wave number k z becomes a function of three anisotropy parameters (v,d and " (or ?)). To extrapolate wave fields in Fourier domain, we must use not only multiple reference velocities, but also multiple reference s and e s(or ?s). The computational cost subsequently increases significantly. The conventional method (Rousseau, 1997; Baumstein and Anderson, 2003) for selecting those three parameters is first to uniformly sample the vertical velocity v at each depth level, then for each vertical velocity, to uniformly sample a range of d s and e s(or ?s) from their minima to their maxima. The main drawback of this method is that it may under sample the parameters at places where the lateral variations are significant, or it may oversample the parameters at places where the model is smooth. To get a better result it may require a large number of reference parameters, which is impractical for migrating large data sets. Clapp (2004) borrowed the idea of quantization from the field of electrical engineering, and used the 1D Lloyd''s algorithm to select reference velocities for isotropic migration. He demonstrated that reference velocities selected by Lloyd''s algorithm can produce Higher quality images while using fewer reference velocities. In fact, the 1D Lloyd''s algorithm (Lloyd, 1982) is a special case of the clustering method, which tries to find the global optimum solutions according to some statistical criteria. It can be easily extended to multi-dimensional case to accurately describe the vector field (Clapp, 2006). In this paper we use the 3D version of Lloyd''s algorithm to select the reference-anisotropy parameters for wave- field extrapolation in laterally varying VTI media. We show that, by incorporating the 3D Lloyd''s algorithm, we significantly reduce the computational cost and obtain high quality images.