Zero-lag Cross-correlation Research Articles

Consistency of the inverse scattering imaging condition, the energy norm imaging condition and the impedance kernel in acoustic and elastic reverse-time migration

Abstract In this work, we draw connections between the imaging conditions using the impedance kernel, the inverse scattering imaging condition, and the energy norm imaging condition in acoustic and elastic reverse-time migration (RTM). Traditional RTM often introduces large low-wavenumber artifacts that degrade image quality in intricate geological structures with large velocity variations. In practice, the Laplacian filter is commonly used to remove these low-wavenumber artifacts, but it changes the image wavenumber spectrum. The advanced imaging conditions of the inverse scattering, the energy norm, and the impedance kernel can effectively remove the low-wavenumber artifacts while not changing the wavenumber spectrum. This study aims to build a connection between these three types of imaging conditions by conducting detailed analysis in the wavenumber domain for acoustic and elastic RTMs. We find that they are exactly the same except for the varying weights of the source-receiver wavefield cross-correlation. All three imaging conditions can generate clear RTM images that are not affected by low-wavenumber artifacts. Numerical examples for a simple model, Sigsbee 2a, and BP models verify the consistency of these three imaging conditions and show their advantage over conventional simple zero-lag cross-correlation imaging conditions. This is important for improving the quality and reliability of seismic imaging technology.

A review of misfit functions for adjoint full waveform inversion in seismology

SUMMARY In seismological full waveform inversion, the choice of misfit functions plays a critical role in quantifying the discrepancy between observed and synthetic data, affecting convergence rate and also final results. We revisit and compare six commonly used misfit functions, including cross-correlation time-shift (CC), least-square waveform difference (L2), multitaper time-shift (MT), exponentiated phase shift (EP), time–frequency phase shift (TF) and zero-lag cross-correlation coefficient (CCC), with respect to their definitions, adjoint sources and misfit kernels for velocity perturbations. Synthetic tests are performed for several canonical models. First, we simulated wave propagation in a model with a single rectangular anomaly with sharp boundaries and a smoothed variant of that model. We analysed the resulting misfit kernels first for the P-wave phase, which is highly distorted in the sharp model due to strong heterogeneities, and mostly experiences traveltime perturbations in the smooth model. Second, we considered a model where a laterally limited region is subject to layered anomalies (low velocity in the middle crust and high velocity in the lower crust) and determine misfit kernels for S and surface waves in this model. Based on these two simplified seismological scenarios, we further perform iterative test inversions using different misfit functions. Combining the features of misfit kernels and synthetic inversion results, we find that CCC, L2 and EP are the most effective at identifying the sharpness of velocity anomalies from the direct body waves and their scattered phases. Consequently, inversion based on those misfit measures yielded the best recovery in the inversion test. For surface and S waves from crustal sources, TF appears to be the most effective in constraining the heterogeneous structure in the crust but needs more iterations for convergence than other misfit functions.

Shear-Wave Anisotropy Measurements in the Crust from Receiver Functions: An Interplay of Lower and Upper Crustal Anisotropy

We report a study using teleseismic P-wave receiver functions to infer the orientation of the maximum horizontal principal stress from the direction of upper crustal shear-wave velocity anisotropy. We apply an inverse approach using the Neighborhood Algorithm to conduct a nonlinear search, attaining a best-fitting crustal model that includes shear velocity anisotropy. Unlike previous methods reported in the literature, this method is able to distinguish anisotropy in the upper, brittle crust from that in the lower, ductile crust in certain instances. We apply this method to teleseismically recorded earthquakes in the Central Valley of California, the Permian Basin, Texas, northern Oklahoma and sites near the San Andreas Fault in California. Of the forty-one stations to which we apply this method, twenty have a good apparent signal. A misfit calculation is performed by calculating a zero-lag cross-correlation coefficient for each modeled receiver function with the data for a given back azimuth range. While the fast polarization direction in the upper crust of some of these stations aligns with independent indicators of the direction of the maximum horizontal principal stress, the fast direction in the upper crust at other stations does not, apparently indicating that the anisotropy was resulting from a different mechanism.

Model-data-driven P-wave impedance inversion using ResNets and the normalized zero-lag cross-correlation objective function

Model-driven and data-driven inversions are two prominent methods for obtaining P-wave impedance, which is significant in reservoir description and identification. Based on proper initial models, most model-driven methods primarily use the limited frequency bandwidth information of seismic data and can invert P-wave impedance with high accuracy, but not high resolution. Conventional data-driven methods mainly employ the information from well-log data and can provide high-accuracy and high-resolution P-wave impedance owing to the superior nonlinear curve fitting capacity of neural networks. However, these methods require a significant number of training samples, which are frequently insufficient. To obtain P-wave impedance with both high accuracy and high resolution, we propose a model-data-driven inversion method using ResNets and the normalized zero-lag cross-correlation objective function which is effective for avoiding local minima and suppressing random noise. By using initial models and training samples, the proposed model-data-driven method can invert P-wave impedance with satisfactory accuracy and resolution. Tests on synthetic and field data demonstrate the proposed method's efficacy and practicability.

Improving the memory efficiency of RTM using both Nyquist sampling and DCT based on GPU

Abstract Reverse time migration (RTM) is used to obtain complex structural images of subsurface media. The RTM can be expressed as a zero-lag cross-correlation between the source and receiver wave fields. As imaging conditions can be calculated based on the pre-stored source wavefield and the received wavefield generated during backward modelling, only the source wavefield must be stored in the memory (or disc). Therefore, reducing source-wavefield storage requirements can improve memory efficiency. High-performance computing based on graphic processing units (GPUs) is being developed to reduce the computational time in wave-propagation modelling. Accordingly, GPU-based RTM technology has the potential to improve the computational efficiency of RTM. Storage of the source wavefield wholly in GPU video random-access memory (VRAM) may further improve computational efficiency. In this paper, we present a new algorithm for a three-dimensional (3D) GPU-based RTM that can enable efficient storage of the source wavefield in VRAM using both Nyquist sampling and discrete cosine transform (DCT) compression. A numerical example employing a modified SEG/EAGE 3D overthrust model presented in this study verifies that the proposed algorithm requires only 2% of the memory usage of conventional RTM while producing similar results.

Elastic geometric-mean reverse time migration for source imaging

Passive source location is a fundamental problem in earthquake seismology and reservoir characterization with fluid injection. It has been shown that geometric-mean reverse time migration (GmRTM) can produce source images with improved resolution and fewer migration artifacts than conventional time-reversal imaging in acoustic media. However, the acoustic assumption often requires data windowing and polarity correction when we apply it to the elastic data, which involves manual efforts and introduces uncertainties. Therefore, we have extended this crosscorrelation-based GmRTM imaging method to elastic media. Back-propagation and Helmholtz decomposition are conducted individually for each receiver wavefield, after which we apply zero-lag crosscorrelations over the decoupled P or/and S wavefields to obtain the corresponding P-, S-, and PS-source images. Considering the influences of random noise, anisotropy, and interferences from multiple sources, we test the proposed elastic GmRTM method on a synthetic example with a layered model and determine its advantages over other imaging conditions (e.g., time-reversal imaging). The method is further validated using another synthetic test based on the elastic Marmousi model and a field-data example from Oklahoma.

Validation of seismic velocity models in southern California with full-waveform simulations

SUMMARY Crustal seismic velocity models provide essential information for many applications including earthquake source properties, simulations of ground motion and related derivative products. We present a systematic workflow for assessing the accuracy of velocity models with full-waveform simulations. The framework is applied to four regional seismic velocity models for southern California: CVM-H15.11, CVM-S4.26, CVM-S4.26.M01 that includes a shallow geotechnical layer, and the model of Berg et al. For each model, we perform 3-D viscoelastic wave propagation simulations for 48 virtual seismic noise sources (down to 2 s) and 44 moderate-magnitude earthquakes (down to 2 s generally and 0.5 s for some cases) assuming a minimum shear wave velocity of 200 m s–1. The synthetic waveforms are compared with observations associated with both earthquake records and noise cross-correlation data sets. We measure, at multiple period bands for well-isolated seismic phases, traveltime delays and normalized zero-lag cross-correlation coefficients between the synthetic and observed data. The obtained measurements are summarized using the mean absolute derivation of time delay and the mean correlation coefficient. These two metrics provide reliable statistical representations of model quality with consistent results in all data sets. In addition to assessing the overall (average) performance of different models in the entire study area, we examine spatial variations of the models’ quality. All examined models show good phase and waveform agreements for surface waves at periods longer than 5 s, and discrepancies at shorter periods reflecting small-scale heterogeneities and near-surface structures. The model performing best overall is CVM-S4.26.M01. The largest misfits for both body and surface waves are in basin structures and around large fault zones. Inaccuracies generated in these areas may affect tomography and model simulation results at other regions. The seismic velocity models for southern California can be improved by adding better resolved structural representations of the shallow crust and volumes around the main faults.

2-D and 3-D Image-Domain Least-Squares Reverse Time Migration Through Point Spread Functions and Excitation-Amplitude Imaging Condition

The enormous computational overheads and excessive storage requirements are two obstacles to the data-domain least-squares reverse time migration (RTM) approach for the application of large-scale 3-D seismic data. To alleviate this problem, we have developed an image-domain least-squares RTM (IDLSRTM) approach through point spread functions (PSFs) and excitation-amplitude (EA) imaging condition, denoted as EA-IDLSRTM. The key point is that the EA imaging condition, as a cost-effective and practical imaging condition, is used to reconstruct the RTM image and localized PSFs. There are two benefits to this combination. One is that the EA imaging condition can effectively reconstruct the RTM image and localized PSFs with less computational overhead and storage requirement, relative to the zero-lag cross correlation (CC) imaging condition. Another important benefit is that the redundant source wavelets in both the RTM and PSF images computed by the CC imaging condition can be removed by the EA imaging condition, prior to the image-domain inversion. As a result, the proposed approach can explicitly reduce the condition number of the Hessian matrix used in the conventional IDLSRTM approach, which will produce a less ill-conditioned inverse problem. In addition, we introduce an angle-dependent filter for the attenuation of low-wavenumber artifacts to accelerate the convergence. Several experiments with synthetic and field data demonstrate that the proposed EA-IDLSRTM approach can efficiently and effectively recover the high-resolution and high-fidelity reflectivity image. Meanwhile, EA-IDLSRTM can provide better imaging quality than the conventional IDLSRTM approach in the case of relatively smoothed velocity.

Parsimonious truncated Newton method for time-domain full-waveform inversion based on the Fourier-domain full-scattered-field approximation

To exploit Hessian information in full-waveform inversion (FWI), the matrix-free truncated Newton method can be used. In such a method, Hessian-vector product computation is one of the major concerns due to the huge memory requirements and demanding computational cost. Using the adjoint-state method, the Hessian-vector product can be estimated by zero-lag crosscorrelation of the first-/second-order incident wavefields and the second-/first-order adjoint wavefields. Different from the implementation in frequency-domain FWI, Hessian-vector product construction in the time domain becomes much more challenging because it is not affordable to store all of the time-dependent wavefields. The widely used wavefield recomputation strategy leads to computationally intensive tasks. We have developed an efficient alternative approach to computing the Hessian-vector product for time-domain FWI. In our method, discrete Fourier transform is applied to extract frequency-domain components of involved wavefields, which are used to compute wavefield crosscorrelation in the frequency domain. This makes it possible to avoid reconstructing the first- and second-order incident wavefields. In addition, a full-scattered-field approximation is proposed to efficiently simplify the second-order incident and adjoint wavefield computation, which enables us to refrain from repeatedly solving the first-order incident and adjoint equations for the second-order incident and adjoint wavefields (re)computation. With our method, the computational time can be reduced by 70% and 80% in viscous media for Gauss-Newton and full-Newton Hessian-vector product construction, respectively. The effectiveness of our method is also verified in the frame of a 2D multiparameter inversion, in which our method almost reaches the same iterative convergence of the conventional time-domain implementation.

Lower extremity kinematic coupling during single and double leg landing and gait in female junior athletes with dynamic knee valgus

BackgroundDynamic knee valgus (DKV) is a common lower extremity movement disorder among females. This study aimed to investigate kinematic couplings between lower extremity joints in female junior athletes with DKV during single and double-leg landing and gait.MethodsTwenty-six physically active female junior athletes (10–14 years old) with DKV were recruited. Kinematic couplings between rearfoot, tibia, knee, and hip were extracted using eight Vicon motion capture cameras and two force plates. Zero-lag cross-correlation coefficient and vector coding were used to calculate kinematic couplings between joints during physical tasks. Paired t-test and Wilcoxon tests were run to find significant couplings between joint motions and coupling strengths. Bonferroni posthoc was used to determine significance with α ≤ 0.05.ResultsThe results showed that the strongest kinematic relationship existed between rearfoot eversion/inversion and tibial internal/external rotation during all three tasks. Correlations of the rearfoot supination/pronation with tibial rotations, knee, and hip motions in sagittal, frontal, and transverse planes were very strong to strong during double-leg landing and moderate to weak during gait. A weak correlation was observed between rearfoot supination/pronation and hip adduction/abduction during single-leg landing.ConclusionsCoupling relationships between rearfoot, knee, and hip vary by the task intensity and alignment profiles in female juniors with DKV.

Memory-efficient source wavefield reconstruction and its application to 3D reverse time migration

Reverse time migration (RTM) generally uses the zero-lag crosscorrelation imaging condition, requiring the source and receiver wavefields to be known at the same time step. However, the receiver wavefield is calculated in time-reversed order, opposite to the order of the forward-propagated source wavefield. The inconvenience can be resolved by storing the source wavefield on a computer memory/disk or by reconstructing the source wavefield on the fly for multiplication with the receiver wavefield. The storage requirements for the former approach can be very large. Hence, we have followed the latter route and developed an efficient source wavefield reconstruction method. During forward propagation, the boundary wavefields at N layers of the spatial grid points and a linear combination of wavefields at [Formula: see text] layers of the spatial grid points are stored. During backward propagation, it reconstructs the source wavefield using the saved wavefields based on a new finite-difference stencil ( M is the operator length parameter, and [Formula: see text]). Unlike existing methods, our method allows a trade-off between accuracy and storage by adjusting N. A maximum-norm-based objective function is constructed to optimize the reconstruction coefficients based on the minimax approximation using the Remez exchange algorithm. Dispersion and stability analyses reveal that our method is more accurate and marginally less stable than the method that requires storage of a combination of boundary wavefields. Our method has been applied to 3D RTM on synthetic and field data. Numerical examples indicate that our method with [Formula: see text] can produce images that are close to those obtained using a conventional method of storing M layers of boundary wavefields. The memory usage of our method is [Formula: see text] times that of the conventional method.

Angle-Weighted Reverse Time Migration With Wavefield Decomposition Based on the Optical Flow Vector

Reverse time migration (RTM) is based on the two-way wave equation, so its imaging results obtained by conventional zero-lag cross-correlation imaging conditions contain a lot of low-wavenumber noises. So far, the wavefield decomposition method based on the Poynting vector has been developed to suppress these noises; however, this method also has some problems, such as unstable calculation of the Poynting vector, low accuracy of wavefield decomposition, and poor effect of large-angle migration artifacts suppression. This article introduces the optical flow vector method to RTM to realize high-precision wavefield decomposition for both the source and receiver wavefields and obtains four directions of wavefields: up-, down-, left-, and right-going. Then, the cross-correlation imaging sections of one-way propagation components of forward- and back-propagated wavefields are optimized and stacked. On this basis, the reflection angle of each imaging point is calculated based on the optical flow vector, and an attenuation factor related to the reflection angle is introduced as the weight to generate the optimal stack images. The tests of theoretical model and field marine seismic data illustrate that compared with the conventional RTM with wavefield decomposition based on the Poynting vector, the angle-weighted RTM with wavefield decomposition based on the optical flow vector proposed in this article can achieve wavefield decomposition for both the source and receiver wavefields and calculate the reflection angle of each imaging point more accurately and stably. Moreover, the proposed method adopts angle weighting processing, which can further eliminate large-angle migration artifacts and effectively improve the imaging accuracy of RTM.

Prestack correlative elastic least-squares reverse time migration based on wavefield decomposition

Standard elastic least-squares reverse time migration (ELSRTM) can iteratively improve the quality of stacked reflection images of P- and S-wave impedance by minimizing the L2 norm of waveform residual. However, the resulting images are sensitive to the velocity errors and amplitude errors in the observed data. Moreover, the crosstalk artifacts created by the different wave modes may degrade the imaging quality. To mitigate these problems, a correlative ELSRTM approach based on wavefield decomposition and prestack implementation is developed. This approach includes two key points. The first is that the proposed approach aims at finding the optimal reflection images for each shot recording via minimizing the normalized zero-lag cross-correlation objective function. The final reflection images can be produced by summing the optimal reflection images for all shot recordings. The second is that we develop a wavefield decomposition scheme, which is compatible with the correlative ELSRTM approach, as a gradient precondition to reduce crosstalk artifacts and accelerate the convergence. Numerical examples using several 2D experiments have determined that the proposed ELSRTM approach can provide higher quality images with fewer crosstalk artifacts, compared to the correlative ELSRTM approach without wavefield decomposition. Furthermore, it is less sensitive to velocity errors than the poststack correlative ELSRTM approach.

Correlative least-squares reverse time migration in viscoelastic media

Standard elastic least-squares reverse time migration (LSRTM) generally ignores the subsurface attenuation effects which can lead to blurred and unfocused migration images. To correct these attenuation effects, we develop a correlative Q-LSRTM approach in viscoelastic media. We use the viscoelastic wave equation based on the standard linear solid model to simulate the subsurface attenuation. The Born modeling and adjoint state equations in viscoelastic media can be derived by the Born approximation and adjoint state method, respectively. Instead of the conventional amplitude-matching misfit function, the proposed Q-LSRTM approach utilizes the normalized zero-lag crosscorrelation misfit function to emphasize the phase measure between the predicted and observed data. The gradients of the crosscorrelation misfit function with regard to reflectivity images of P- and S-wave impedance are estimated by a zero-lag crosscorrelation between the forward background wavefields and the reversal adjoint wavefields. In view of computational efficiency, we also derive an analytical step-length formula for the correlative Q-LSRTM approach, which only requires to solve the Born modeling operator once. Numerical examples using several 2D models have determined that the proposed Q-LSRTM can improve better image quality, compared with the standard elastic LSRTM approach. Furthermore, the crosscorrelation misfit function is less sensitive to amplitude errors.

A reverse time migration-based multistep angular spectrum approach for ultrasonic imaging of specimens with irregular surfaces

We develop a new ultrasonic imaging framework for non-destructive testing of an immersed specimen featuring an irregular top surface and demonstrate its capability of accurately depicting the lower surfaces of multiple damages hidden in the specimen. Central to the framework is a multistep angular spectrum approach (ASA), via which the forward propagation wavefields of wave sources and backward propagation wavefields of the received wave signals are calculated. Upon applying a zero-lag cross-correlation imaging condition of reverse time migration (RTM) to the obtained forward and backward wavefields, the image of the specimen with an irregular surface can be reconstructed, in which hidden damages, if any and regardless of quantity, are visualized. The effectiveness and accuracy of the framework are examined using numerical simulation, followed with experiment, in both of which multiple side-drilled holes, at different locations in aluminum blocks with various irregular surfaces, are characterized. Results have proven that multiple damages in a specimen with an irregular surface can be individually localized, and the lower surface of each damage can further be imaged accurately, thanks to the RTM-based algorithm in which multiple wave reflections from the specimen bottom are taken into wavefield extrapolation. The proposed imaging approach presents higher computational efficiency, compared to conventional RTM, and enhanced imaging contrast over prevailing total focusing methods.

Elastodynamic full waveform inversion on GPUs with time-space tiling and wavefield reconstruction

Full waveform inversion (FWI) is a procedure used to determine the elastic parameters of the Earth by reducing the misfit between observed elastodynamic wavefields and their numerically modelled counterparts. The numerical solution of the elastodynamic wave equation is computationally expensive, and its performance is typically bandwidth bound. Computing the gradient of the FWI misfit functional adds further complexity as it involves computing the zero-lag cross-correlation of two wavefields propagating in opposite temporal directions. In this paper, we utilize graphics processing units (GPUs) for their high memory bandwidth and combine two principal optimizations in order to compute FWI gradients on large models and for long simulation times. Wavefield reconstruction methods allow efficient gradient computations with minimal memory requirements and interconnection transfers. Time-space tiling techniques permit us to transcend the limited amount of GPU memory while avoiding dramatic slowdowns due to the low interconnection bandwidth. The implementation considers a task-oriented, hybrid usage of explicitly managed and Unified Memory in order to satisfy the requirements. Benchmarks demonstrate that the proposed approach is able to preserve 78–90% of the original performance, when oversubscribing the amount of physical memory available on GPUs. Comparison with existing methods highlights the benefits of the method.

Wavefield reconstruction for velocity–stress elastodynamic full-waveform inversion

SUMMARY Gradient computations in full-waveform inversion (FWI) require calculating zero-lag cross-correlations of two wavefields propagating in opposite temporal directions. Lossless media permit accurate and efficient reconstruction of the incident field from recordings along a closed boundary, such that both wavefields propagate backwards in time. Reconstruction avoids storing wavefield states of the incident field to secondary storage, which is not feasible for many realistic inversion problems. We give particular attention to velocity–stress modelling schemes and propose a novel modification of a conventional reconstruction method derived from the elastodynamic Kirchhoff–Helmholtz integral. In contrast to the original formulation (in a previous related work), the proposed approach is well-suited for velocity–stress schemes. Numerical examples demonstrate accurate wavefield reconstruction in heterogeneous, elastic media. A practical example using 3-D elastic FWI demonstrates agreement with the reference solution.

Elastic Least-Squares Imaging in Tilted Transversely Isotropic Media for Multicomponent Land and Pressure Marine Data

Traditional elastic reverse-time migration (RTM) involves P-/S-wave separation for the source and receiver wavefields, followed by applying the zero-lag cross-correlation imaging condition to produce PP and PS images. In anisotropic media, P-/S-wave decomposition requires a higher memory and computational cost than that in isotropic media. In addition, finite acquisition apertures and band-limited source functions result in unsatisfactory resolutions and amplitudes. To mitigate these problems, we present an elastic least-squares imaging method for tilted transversely isotropic media and apply it to land multicomponent and marine pressure data. Unlike traditional RTM, we use the relative perturbations to the product of density and squared axial (compressional/shear) velocities as reflectivity models ($$\Delta \ln{C}_{33}$$ and $$\Delta \ln{C}_{55}$$), and estimate them by solving a linear inverse problem. Numerical experiments illustrate that subsurface reflectors can be well resolved in adjoint images for land multicomponent data, because of the presence of both P- and S-waves in seismograms. Least-squares migration helps to further improve spatial resolution and image amplitudes. Since there are no direct S-waves in marine streamer data, adjoint RTM images of $$\Delta \ln{C}_{55}$$ are mainly resolved with the converted S-waves and are not as good as those in $$\Delta \ln{C}_{33}$$ images. By approximating the Hessian inverse, least-squares migration allows us to take advantage of the weak converted P–S–P-waves and improve the $$\Delta \ln{C}_{55}$$ image quality. Numerical experiments for synthetic and field data demonstrate the feasibility and advantage of the proposed least-squares TTI RTM compared with wave-mode separation-based elastic RTM. In field data experiments, we observe that since there are no strong P–S–P converted waves in streamer pressure records from the marine survey, the reflectors in $$\Delta \ln{C}_{55}$$ image might be mainly imaged from P-waves due to the amplitude versus offset (AVO) effects.

Energy Flow Domain Reverse-Time Migration for Borehole Radar

A modified 2-D reverse-time migration algorithm in the energy flow domain, called EF-RTM, is proposed for short impulse borehole radar (BHR) imaging. The key of the approach is based on Poynting's theorem, which allows decomposition of the energy flux density (EFD) derived from the source and receiver wave fields in different wave-propagation directions. Then, imaging conditions, for example, zero-lag cross correlation (prestack migration) and zero-time imaging principle (poststack migration), are applied to the decomposed EFD-field components to obtain the migrated sections. The characteristics of the resulting images can be optionally combined or separately used for better BHR data imaging, interpretation, and analysis. In this paper, the EF-RTM algorithm is validated by numerical modeling and real field data and compared with the conventional RTM algorithm. All the results show that this new EF-RTM method is superior to the conventional RTM method: it inherits the high precision of RTM with additional imaging advantages, such as natural wave-field decomposition in different directions, improved migrated cross-range resolution, a better focus of the target's shape, and migration noise reduction.

3D wave-equation dispersion inversion of Rayleigh waves

The 2D wave-equation dispersion (WD) inversion method is extended to 3D wave-equation dispersion inversion of surface waves for the shear-velocity distribution. The objective function of 3D WD is the frequency summation of the squared wavenumber [Formula: see text] differences along each azimuth angle of the fundamental or higher modes of Rayleigh waves in each shot gather. The S-wave velocity model is updated by the weighted zero-lag crosscorrelation between the weighted source-side wavefield and the back-projected receiver-side wavefield for each azimuth angle. A multiscale 3D WD strategy is provided, which starts from the pseudo-1D S-velocity model, which is then used to get the 2D WD tomogram, which in turn is used as the starting model for 3D WD. The synthetic and field data examples demonstrate that 3D WD can accurately reconstruct the 3D S-wave velocity model of a laterally heterogeneous medium and has much less of a tendency to getting stuck in a local minimum compared with full-waveform inversion.