Search

We have recast the forward pass of a multilayered convolutional neural network (CNN) as the solution to the problem of sparse least-squares migration (LSM). The CNN filters and feature maps are shown to be analogous, but not equivalent, to the migration Green’s functions and the quasi-reflectivity distribution, respectively. This provides a physical interpretation of the filters and feature maps in deep CNN in terms of the operators for seismic imaging. Motivated by the connection between sparse LSM and CNN, we adopt the neural network version of sparse LSM. Unlike the standard LSM method that finds the optimal reflectivity image, neural network LSM (NNLSM) finds the optimal quasi-reflectivity image and the quasi-migration Green’s functions. These quasi-migration Green’s functions are also denoted as the convolutional filters in a CNN and are similar to migration Green’s functions. The advantage of NNLSM over standard LSM is that its computational cost is significantly less and it can be used for denoising coherent and incoherent noise in migration images. Its disadvantage is that the NNLSM quasi-reflectivity image is only an approximation to the actual reflectivity distribution. However, the quasi-reflectivity image can be used as an attribute image for high-resolution delineation of geologic bodies.

Irregular topography can cause strong scattering and defocusing of propagating surface waves, so it is important to account for such effects when inverting surface waves for shallow S-wave velocity structures. We have developed a 3D surface-wave dispersion inversion method that takes into account the topographic effects modeled by a 3D spectral element solver. The objective function is the frequency summation of the squared wavenumber differences [Formula: see text] along each azimuthal angle of the fundamental mode or higher-order modes of Rayleigh waves in each shot gather. The wavenumbers [Formula: see text] associated with the dispersion curves are calculated using the data recorded along the irregular free surface. Numerical tests on synthetic and field data demonstrate that 3D topographic wave equation dispersion inversion (TWD) can accurately invert for the S-wave velocity model from surface-wave data recorded on irregular topography. Field data tests for data recorded across an Arizona fault demonstrate that, for this example, the 2D TWD model can be as accurate as the 3D tomographic model. This suggests that in some cases, the 2D TWD inversion is preferred over 3D TWD because of its significant reduction in computational costs. Compared to the 3D P-wave velocity tomogram, the 3D S-wave tomogram agrees much more closely with the geologic model taken from the trench log. The agreement with the trench log is even better when the [Formula: see text] tomogram is computed, which reveals a sharp change in velocity across the fault. The localized velocity anomaly in the [Formula: see text] tomogram is in very good agreement with the well log. Our results suggest that integrating the [Formula: see text] and [Formula: see text] tomograms can sometimes give the most accurate estimates of the subsurface geology across normal faults.

To image deeper portions of the earth, geophysicists must record reflection data with much greater source-receiver offsets. The problem with these data is that the signal-to-noise ratio (S/N) significantly diminishes with greater offset. In many cases, the poor S/N makes the far-offset reflections imperceptible on the shot records. To mitigate this problem, we have developed supervirtual reflection interferometry (SVI), which can be applied to far-offset reflections to significantly increase their S/N. The key idea is to select the common pair gathers where the phases of the correlated reflection arrivals differ from one another by no more than a quarter of a period so that the traces can be coherently stacked. The traces are correlated and summed together to create traces with virtual reflections, which in turn are convolved with one another and stacked to give the reflection traces with much stronger S/Ns. This is similar to refraction SVI except far-offset reflections are used instead of refractions. The theory is validated with synthetic tests where SVI is applied to far-offset reflection arrivals to significantly improve their S/N. Reflection SVI is also applied to a field data set where the reflections are too noisy to be clearly visible in the traces. After the implementation of reflection SVI, the normal moveout velocity can be accurately picked from the SVI-improved data, leading to a successful poststack migration for this data set.

Reflection full-waveform inversion (RFWI) can recover the low-wavenumber components of the velocity model along with the reflection wavepaths. However, this requires an expensive least-squares reverse time migration (LSRTM) to construct the perturbation image that can still suffer from cycle-skipping problems. As an inexpensive alternative to LSRTM, we use migration deconvolution (MD) with RFWI. To mitigate cycle-skipping problems, we develop a multiscale reflection phase inversion (MRPI) strategy that boosts the low-frequency data and should only explain the phase information in the recorded data, not its magnitude spectrum. We also use the rolling-offset strategy that gradually extends the offset range of data with an increasing number of iterations. Numerical results indicate that the MRPI + MD method can efficiently recover the low-wavenumber components of the velocity model and is less prone to getting stuck in local minima compared to conventional RFWI.

We present a quasi-elastic wave equation as a function of the pressure variable, which can accurately model PP reflections with elastic amplitude variation with offset effects under the first-order Born approximation. The kinematic part of the quasi-elastic wave equation accurately models the propagation of P waves, whereas the virtual-source part, which models the amplitudes of reflections, is a function of the perturbations of density and Lamé parameters [Formula: see text] and [Formula: see text]. The quasi-elastic wave equation generates a scattering radiation pattern that is exactly the same as that for the elastic wave equation, and only requires the solution of two acoustic wave equations for each shot gather. This means that the quasi-elastic wave equation can be used for true-amplitude linearized waveform inversion (also known as least-squares reverse time migration) of elastic PP reflections, in which the corresponding misfit gradients are with respect to the perturbations of density and the P- and S-wave impedances. The perturbations of elastic parameters are iteratively updated by minimizing the [Formula: see text]-norm of the difference between the recorded PP reflections and the predicted pressure data modeled from the quasi-elastic wave equation. Numerical tests on synthetic and field data indicate that true-amplitude linearized waveform inversion using the quasi-elastic wave equation can account for the elastic PP amplitudes and provide a robust estimate of the perturbations of P- and S-wave impedances and, in some cases, the density. In addition, true-amplitude linearized waveform inversion provides images with a wider bandwidth and fewer artifacts because the PP amplitudes are accurately explained. We also determine the 2D scalar quasi-elastic wave equation for P-SV reflections and the 3D vector equation for PS reflections.

We present a theory for wave-equation inversion of Love-wave dispersion curves, in which the misfit function is the sum of the squared differences between the wavenumbers along the predicted and observed dispersion curves. Similar to inversion of Rayleigh-wave dispersion curves, the complicated Love-wave arrivals in traces are skeletonized as simpler data, namely, the picked dispersion curves in the [Formula: see text] domain. Numerical solutions to the SH-wave equation and an iterative optimization method are then used to invert these dispersion curves for the S-wave velocity model. This procedure, denoted as wave-equation dispersion inversion of Love waves (LWD), does not require the assumption of a layered model or smooth velocity variations, and it is less prone to the cycle-skipping problems of full-waveform inversion. We demonstrate with synthetic and field data examples that LWD can accurately reconstruct the S-wave velocity distribution in a laterally heterogeneous medium. Compared with Rayleigh waves, inversion of the Love-wave dispersion curves empirically exhibits better convergence properties because they are completely insensitive to the P-velocity variations. In addition, Love-wave dispersion curves for our examples are simpler than those for Rayleigh waves, and they are easier to pick in our field data with a low signal-to-noise ratio.

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.

Viscoacoustic migration can significantly compensate for the amplitude loss and phase distortion in migration images computed from highly attenuated data. However, solving the viscoacoustic wave equation requires a significant amount of storage space and computation time, especially for least-squares migration methods. To mitigate this problem, we used acoustic reverse time migration (RTM) instead of viscoacoustic migration to migrate the viscoacoustic data and then we correct the amplitude and phase distortion by hybrid deblurring filters in the image domain. Numerical tests on synthetic and field data demonstrate that acoustic RTM combined with hybrid deblurring filters can compensate for the attenuation effects and produce images with high resolution and balanced amplitudes. This procedure requires less than one-third of the storage space and is [Formula: see text] times faster compared with the viscoacoustic migration, but at the cost of mildly reduced accuracy. Here, [Formula: see text] represents the number of iterations used for least-squares migration method. This method can be extended to 3D migration at even a greater cost saving.

Velocity errors in the shallow part of the velocity model can lead to erroneous estimates of the full-waveform inversion (FWI) tomogram. If the location and topography of a reflector are known, then such a reflector can be used as a reference reflector to update the underlying velocity model. Reflections corresponding to this reference reflector are windowed in the data space. Windowed reference reflections are then crosscorrelated with reflections from deeper interfaces, which leads to partial cancellation of static errors caused by the overburden above the reference interface. Interferometric FWI (IFWI) is then used to invert the tomogram in the target region, by minimizing the normalized waveform misfit between the observed and predicted crosscorrelograms. Results with synthetic and field data with static errors above the reference interface indicate that an accurate tomogram can be inverted in areas lying within several wavelengths of the reference interface. IFWI can also be applied to synthetic time-lapse data to mitigate the nonrepeatability errors caused by time-varying overburden variations. The synthetic- and field-data examples demonstrate that IFWI can provide accurate tomograms when the near surface is ridden with velocity errors.

We have developed a wave-equation traveltime inversion method with multifrequency bands to invert for the shallow or intermediate subsurface velocity distribution. Similar to the classical wave-equation traveltime inversion, this method searches for the velocity model that minimizes the squared sum of the traveltime residuals using source wavelets with progressively higher peak frequencies. Wave-equation traveltime inversion can partially avoid the cycle-skipping problem by recovering the low-wavenumber parts of the velocity model. However, we also use the frequency information hidden in the traveltimes to obtain a more highly resolved tomogram. Therefore, we use different frequency bands when calculating the Fréchet derivatives so that tomograms with better resolution can be reconstructed. Results are validated by the zero-offset gathers from the raw data associated with moderate geometric irregularities. The improved wave-equation traveltime method is robust and merely needs a rough estimate of the starting model. Numerical tests on the synthetic and field data sets validate the above claims.

Different types of model parameters, such as the P- and S-wave velocities, can be coupled to one another in multiparameter seismic inversion. This coupling effect is not taken into account by the conventional approximation to the Hessian inverse for a single type of parameter, and so it can lead to a slow rate of convergence by an iterative inversion method. To solve this problem, we have developed a multiparameter deblurring filter that approximates the Hessian inverse. This filter takes into account the coupling between different parameters by using stationary local filters to approximate the submatrices of the Hessian inverse for the different types of parameters. The filters are calculated by matching the reference multiparameter migration images to their reference reflectivity models. Numerical tests with elastic migration and inversion indicate that the multiparameter deblurring filter not only reduces the footprint noise, balances the amplitude, and increases the resolution of the elastic migration images, but it also mitigates the crosstalk artifacts caused by the coupling effect. When used as a preconditioner, it also accelerates the convergence rate for elastic inversion.

A resonant multiple is defined as a multiple reflection that revisits the same subsurface location along coincident reflection raypaths. We have found that resonant first-order multiples can be migrated with either Kirchhoff or wave-equation migration methods to give images with approximately twice the spatial resolution compared with poststack primary-reflection images. We have developed a moveout-correction stacking method to enhance the signal-to-noise ratios of the resonant multiples for superresolution migration. The effectiveness of this procedure is validated by synthetic and field data tests.

We present a scheme for multiscale phase inversion (MPI) of seismic data that is less sensitive than full-waveform inversion (FWI) to the unmodeled physics of wave propagation and to a poor starting model. To avoid cycle skipping, the multiscale strategy temporally integrates the traces several times, i.e., high-order integration, to produce low-boost seismograms that are used as input data for the initial iterations of MPI. As the iterations proceed, lower frequencies in the data are boosted by using integrated traces of lower order as the input data. The input data are also filtered into different narrow frequency bands for the MPI implementation. Numerical results with synthetic acoustic data indicate that, for the Marmousi model, MPI is more robust than conventional multiscale FWI when the initial model is moderately far from the true model. Results from synthetic viscoacoustic and elastic data indicate that MPI is less sensitive than FWI to some of the unmodeled physics. Inversion of marine data indicates that MPI is more robust and produces modestly more accurate results than FWI for this data set.

Wave-equation migration velocity analysis (WEMVA) based on subsurface-offset, angle domain, or time-lag common-image gathers (CIGs) requires significant computational and memory resources because it computes higher dimensional migration images in the extended image domain. To mitigate this problem, we have developed a WEMVA method using plane-wave CIGs. Plane-wave CIGs reduce computational cost and memory storage because they are directly calculated from prestack plane-wave migration and the number of plane waves is often much smaller than the number of shots. In the case of an inaccurate migration velocity, the moveout of plane-wave CIGs is automatically picked by a semblance analysis method, which is then linked to the migration velocity update by a connective function. Numerical tests on two synthetic data sets and a field data set validate the efficiency and effectiveness of this method.

We have developed a methodology for detecting the presence of near-surface heterogeneities by naturally migrating backscattered surface waves in controlled-source data. The near-surface heterogeneities must be located within a depth of approximately one-third the dominant wavelength [Formula: see text] of the strong surface-wave arrivals. This natural migration method does not require knowledge of the near-surface phase-velocity distribution because it uses the recorded data to approximate the Green’s functions for migration. Prior to migration, the backscattered data are separated from the original records, and the band-passed filtered data are migrated to give an estimate of the migration image at a depth of approximately one-third [Formula: see text]. Each band-passed data set gives a migration image at a different depth. Results with synthetic data and field data recorded over known faults validate the effectiveness of this method. Migrating the surface waves in recorded 2D and 3D data sets accurately reveals the locations of known faults. The limitation of this method is that it requires a dense array of receivers with a geophone interval less than approximately one-half [Formula: see text].

Least-squares migration (LSM) can produce images with better balanced amplitudes and fewer artifacts than standard migration. The conventional objective function used for LSM minimizes the L2-norm of the data residual between the predicted and the observed data. However, for field-data applications in which the recorded data are noisy and undersampled, the conventional formulation of LSM fails to provide the desired uplift in the quality of the inverted image. We have developed a least-squares reverse time migration (LSRTM) method using local Radon-based preconditioning to overcome the low signal-to-noise ratio (S/N) problem of noisy or severely undersampled data. A high-resolution local Radon transform of the reflectivity is used, and sparseness constraints are imposed on the inverted reflectivity in the local Radon domain. The sparseness constraint is that the inverted reflectivity is sparse in the Radon domain and each location of the subsurface is represented by a limited number of geologic dips. The forward and the inverse mapping of the reflectivity to the local Radon domain and vice versa is done through 3D Fourier-based discrete Radon transform operators. The weights for the preconditioning are chosen to be varying locally based on the relative amplitudes of the local dips or assigned using quantile measures. Numerical tests on synthetic and field data validate the effectiveness of our approach in producing images with good S/N and fewer aliasing artifacts when compared with standard RTM or standard LSRTM.

We use elastic least-squares reverse time migration (LSRTM) to invert for the reflectivity images of P- and S-wave impedances. Elastic LSRTM solves the linearized elastic-wave equations for forward modeling and the adjoint equations for backpropagating the residual wavefield at each iteration. Numerical tests on synthetic data and field data reveal the advantages of elastic LSRTM over elastic reverse time migration (RTM) and acoustic LSRTM. For our examples, the elastic LSRTM images have better resolution and amplitude balancing, fewer artifacts, and less crosstalk compared with the elastic RTM images. The images are also better focused and have better reflector continuity for steeply dipping events compared to the acoustic LSRTM images. Similar to conventional least-squares migration, elastic LSRTM also requires an accurate estimation of the P- and S-wave migration velocity models. However, the problem remains that, when there are moderate errors in the velocity model and strong multiples, LSRTM will produce migration noise stronger than that seen in the RTM images.

Strong subsurface attenuation leads to distortion of amplitudes and phases of seismic waves propagating inside the earth. The amplitude and the dispersion losses from attenuation are often compensated for during prestack depth migration. However, most attenuation compensation or [Formula: see text]-compensation migration algorithms require an estimate of the background [Formula: see text] model. We have developed a wave-equation gradient optimization method that inverts for the subsurface [Formula: see text] distribution by minimizing a skeletonized misfit function [Formula: see text], where [Formula: see text] is the sum of the squared differences between the observed and the predicted peak/centroid-frequency shifts of the early arrivals. The gradient is computed by migrating the observed traces weighted by the frequency shift residuals. The background [Formula: see text] model is perturbed until the predicted and the observed traces have the same peak frequencies or the same centroid frequencies. Numerical tests determined that an improved accuracy of the [Formula: see text] model by wave-equation [Formula: see text] tomography leads to a noticeable improvement in migration image quality.

We have developed an efficient approach for picking first-break wavefronts on coarsely sampled time slices of 3D shot gathers. Our objective was to compute a smooth initial velocity model for multiscale full-waveform inversion (FWI). Using interactive software, first-break wavefronts were geometrically modeled on time slices with a minimal number of picks. We picked sparse time slices, performed traveltime tomography, and then compared the predicted traveltimes with the data in-between the picked slices. The picking interval was refined with iterations until the errors in traveltime predictions fell within the limits necessary to avoid cycle skipping in early arrivals FWI. This approach was applied to a 3D ocean-bottom-station data set. Our results indicate that wavefront picking has 28% fewer data slices to pick compared with picking traveltimes in shot gathers. In addition, by using sparse time samples for picking, data storage is reduced by 88%, and therefore allows for a faster visualization and quality control of the picks. Our final traveltime tomogram is sufficient as a starting model for early arrival FWI.

Imaging seismic data with an erroneous migration velocity can lead to defocused migration images. To mitigate this problem, we first choose a reference reflector whose topography is well-known from the well logs, for example. Reflections from this reference layer are correlated with the traces associated with reflections from deeper interfaces to get crosscorrelograms. Interferometric least-squares migration (ILSM) is then used to get the migration image that maximizes the crosscorrelation between the observed and the predicted crosscorrelograms. Deeper reference reflectors are used to image deeper parts of the subsurface with a greater accuracy. Results on synthetic and field data show that defocusing caused by velocity errors is largely suppressed by ILSM. We have also determined that ILSM can be used for 4D surveys in which environmental conditions and acquisition parameters are significantly different from one survey to the next. The limitations of ILSM are that it requires prior knowledge of a reference reflector in the subsurface and the velocity model below the reference reflector should be accurate.

Full-waveform inversion requires the accurate simulation of the dynamics and kinematics of wave propagation. This is difficult in practice because the amplitudes cannot be precisely reproduced for seismic waves in the earth. Wave-equation reflection traveltime tomography (WT) is proposed to avoid this problem by directly inverting the reflection-traveltime residuals without the use of the high-frequency approximation. We inverted synthetic traces and recorded seismic data for the velocity model by WT. Our results demonstrated that the wave-equation solution overcame the high-frequency approximation of ray-based tomography, was largely insensitive to the accurate modeling of amplitudes, and mitigated problems with ambiguous event identification. The synthetic examples illustrated the effectiveness of the WT method in providing a highly resolved estimate of the velocity model. A real data example from the Gulf of Mexico demonstrated these benefits of WT, but also found the limitations in traveltime residual estimation for complex models.

The main difficulty with iterative waveform inversion is that it tends to get stuck in local minima associated with the waveform misfit function. To mitigate this problem and avoid the need to fit amplitudes in the data, we have developed a wave-equation method that inverts the traveltimes of reflection events, and so it is less prone to the local minima problem. Instead of a waveform misfit function, the penalty function was a crosscorrelation of the downgoing direct wave and the upgoing reflection wave at the trial image point. The time lag, which maximized the crosscorrelation amplitude, represented the reflection-traveltime residual (RTR) that was back projected along the reflection wavepath to update the velocity. Shot- and angle-domain crosscorrelation functions were introduced to estimate the RTR by semblance analysis and scanning. In theory, only the traveltime information was inverted and there was no need to precisely fit the amplitudes or assume a high-frequency approximation. Results with synthetic data and field records revealed the benefits and limitations of wave-equation reflection traveltime inversion.

Strong subsurface attenuation leads to distortion of amplitudes and phases of seismic waves propagating inside the earth. Conventional acoustic reverse time migration (RTM) and least-squares reverse time migration (LSRTM) do not account for this distortion, which can lead to defocusing of migration images in highly attenuative geologic environments. To correct for this distortion, we used a linearized inversion method, denoted as [Formula: see text]-LSRTM. During the least-squares iterations, we used a linearized viscoacoustic modeling operator for forward modeling. The adjoint equations were derived using the adjoint-state method for back propagating the residual wavefields. The merit of this approach compared with conventional RTM and LSRTM was that [Formula: see text]-LSRTM compensated for the amplitude loss due to attenuation and could produce images with better balanced amplitudes and more resolution below highly attenuative layers. Numerical tests on synthetic and field data illustrated the advantages of [Formula: see text]-LSRTM over RTM and LSRTM when the recorded data had strong attenuation effects. Similar to standard LSRTM, the sensitivity tests for background velocity and [Formula: see text] errors revealed that the liability of this method is the requirement for smooth and accurate migration velocity and attenuation models.

The timing and amplitudes of arrivals recorded in seismic traces are influenced by velocity variations all along the associated raypaths. Consequently, velocity errors far from the target can lead to blurred imaging of the target body. To partly remedy this problem, we comprehensively reviewed inverting differential traveltimes that satisfied the closure-phase condition. The result is that the source and receiver statics are completely eliminated in the data and velocities far from the target do not need to be known. We successfully used the phase closure equation for traveltime tomography, refraction statics, migration, refraction tomography, and earthquake location, all of which demonstrated the higher resolution achievable by processing data with differential traveltimes rather than absolute traveltimes. More generally, the stationary version of the closure-phase equation is equivalent to Fermat’s principle and can be derived from the equations of seismic interferometry. In summary, the general closure-phase equation is the mathematical foundation for approximately redatuming sources and/or receivers to the target of interest without the need to accurately know the statics or the velocity model away from the target.

The theory of least-squares reverse time migration of multiples (RTMM) is presented. In this method, least squares migration (LSM) is used to image free-surface multiples where the recorded traces are used as the time histories of the virtual sources at the hydrophones and the surface-related multiples are the observed data. For a single source, the entire free-surface becomes an extended virtual source where the downgoing free-surface multiples more fully illuminate the subsurface compared to the primaries. Since each recorded trace is treated as the time history of a virtual source, knowledge of the source wavelet is not required and the ringy time series for each source is automatically deconvolved. If the multiples can be perfectly separated from the primaries, numerical tests on synthetic data for the Sigsbee2B and Marmousi2 models show that least-squares reverse time migration of multiples (LSRTMM) can significantly improve the image quality compared to RTMM or standard reverse time migration (RTM) of primaries. However, if there is imperfect separation and the multiples are strongly interfering with the primaries then LSRTMM images show no significant advantage over the primary migration images. In some cases, they can be of worse quality. Applying LSRTMM to Gulf of Mexico data shows higher signal-to-noise imaging of the salt bottom and top compared to standard RTM images. This is likely attributed to the fact that the target body is just below the sea bed so that the deep water multiples do not have strong interference with the primaries. Migrating a sparsely sampled version of the Marmousi2 ocean bottom seismic data shows that LSM of primaries and LSRTMM provides significantly better imaging than standard RTM. A potential liability of LSRTMM is that multiples require several round trips between the reflector and the free surface, so that high frequencies in the multiples suffer greater attenuation compared to the primary reflections. This can lead to lower resolution in the migration image compared to that computed from primaries. Another liability is that the multiple migration image is more down-dip limited than the standard primaries migration image. Finally, if the surface-related multiple elimination method is imperfect and there are strong multiples interfering with the primaries, then the resulting LSRTMM image can be significantly degraded. We conclude that LSRTMM can be a useful complement, not a replacement, for RTM of primary reflections.

A plane-wave least-squares reverse-time migration (LSRTM) is formulated with a new parameterization, where the migration image of each shot gather is updated separately and an ensemble of prestack images is produced along with common image gathers. The merits of plane-wave prestack LSRTM are the following: (1) plane-wave prestack LSRTM can sometimes offer stable convergence even when the migration velocity has bulk errors of up to 5%; (2) to significantly reduce computation cost, linear phase-shift encoding is applied to hundreds of shot gathers to produce dozens of plane waves. Unlike phase-shift encoding with random time shifts applied to each shot gather, plane-wave encoding can be effectively applied to data with a marine streamer geometry. (3) Plane-wave prestack LSRTM can provide higher-quality images than standard reverse-time migration. Numerical tests on the Marmousi2 model and a marine field data set are performed to illustrate the benefits of plane-wave LSRTM. Empirical results show that LSRTM in the plane-wave domain, compared to standard reverse-time migration, produces images efficiently with fewer artifacts and better spatial resolution. Moreover, the prestack image ensemble accommodates more unknowns to makes it more robust than conventional least-squares migration in the presence of migration velocity errors.

Least-squares migration (LSM) has been shown to be able to produce high-quality migration images, but its computational cost is considered to be too high for practical imaging. We have developed a multisource least-squares migration algorithm (MLSM) to increase the computational efficiency by using the blended sources processing technique. To expedite convergence, a multisource deblurring filter is used as a preconditioner to reduce the data residual. This MLSM algorithm is applicable with Kirchhoff migration, wave-equation migration, or reverse time migration, and the gain in computational efficiency depends on the choice of migration method. Numerical results with Kirchhoff LSM on the 2D SEG/EAGE salt model show that an accurate image is obtained by migrating a supergather of 320 phase-encoded shots. When the encoding functions are the same for every iteration, the input/output cost of MLSM is reduced by 320 times. Empirical results show that the crosstalk noise introduced by blended sources is more effectively reduced when the encoding functions are changed at every iteration. The analysis of signal-to-noise ratio (S/N) suggests that not too many iterations are needed to enhance the S/N to an acceptable level. Therefore, when implemented with wave-equation migration or reverse time migration methods, the MLSM algorithm can be more efficient than the conventional migration method.

An antialiasing formula has been derived for interferometric redatuming of seismic data. More generally, this formula is valid for numerical implementation of the reciprocity equation of correlation type, which is used for redatuming, extrapolation, interpolation, and migration. The antialiasing condition can be, surprisingly, more tolerant of a coarser trace sampling compared to the standard antialiasing condition. Numerical results with synthetic vertical seismic profile (VSP) data show that interferometry artifacts are effectively reduced when the antialiasing condition is used as a constraint with interferometric redatuming.

A recently developed time-domain multiscale waveform tomography (MWT) method is applied to synthetic and field marine data. Although the MWT method was already applied to synthetic data, the synthetic data application leads to a development of a hybrid method between waveform tomography and the salt flooding technique commonly use in subsalt imaging. This hybrid method can overcome a convergence problem encountered by inversion with a traveltime velocity tomogram and successfully provides an accurate and highly resolved velocity tomogram for the 2D SEG/EAGE salt model. In the application of MWT to the field data, the inversion process is carried out using a multiscale method with a dynamic early-arrival muting window to mitigate the local minima problem of waveform tomography and elastic effects. With the modified MWT method, reasonably accurate results as verified by comparison of migration images and common image gathers were obtained. The hybrid method with the salt flooding technique is not used in this field data example because there is no salt in the subsurface according to our interpretation. However, we believe it is applicable to field data applications.

Least-squares migration (LSM) is a linearized waveform inversion for estimating a subsurface reflectivity model that, relative to conventional migration, improves spatial resolution of migration images. The cost, however, is high because LSM typically requires 10 or more iterations, which is at least 20 times more than the CPU cost of conventional migration. To alleviate this expense, we offer a deblurring filter that can be used in a regularization scheme or a preconditioning scheme to give acceptable LSM images with less than one-third the cost of the standard LSM method. Our results in applying deblurred LSM to synthetic data and field data support this claim. In particular, a Marmousi2 model test shows that the data residual for preconditioned deblurred LSM decreases rapidly in the first iteration, which is equivalent to 10 or more iterations of LSM. Empirical results suggest that regularized deblurred LSM after three iterations is equivalent to about 10 iterations of LSM. Applying deblurred LSM to 2D marine data gives a higher-resolution image compared to those from migration or LSM with three iterations. These results suggest that LSM combined with a deblurring filter allows LSM to be a fast, practical tool for improved imaging of complicated structures.

This efficient multiscale method for time-domain waveform tomography incorporates filters that are more efficient than Hamming-window filters. A strategy for choosing optimal frequency bands is proposed to achieve computational efficiency in the time domain. A staggered-grid, explicit finite-difference method with fourth-order accuracy in space and second-order accuracy in time is used for forward modeling and the adjoint calculation. The adjoint method is utilized in inverting for an efficient computation of the gradient directions. In the multiscale approach, multifrequency data and multiple grid sizes are used to overcome somewhat the severe local minima problem of waveform tomography. The method is applied successfully to 1D and 2D heterogeneous models; it can accurately recover low- and high-wavenumber components of the velocity models. The inversion result for the 2D model demonstrates that the multiscale method is computationally efficient and converges faster than a conventional, single-scale method.

Imaging of subsalt sediments is a challenge for traditional migration methods such as Kirchhoff and one-way wave-equation migration. Consequently, the more accurate two-way method of reverse-time migration (RTM) is preferred for subsalt imaging, but its use can be limited by high computation cost. To overcome this problem, a 3D target-oriented reverse-time datuming (RTD) method is presented, which can generate redatumed data economically in target areas beneath complex structures such as salt domes. The redatumed data in the target area then can be migrated inexpensively using a traditional migration method. If the target area is much smaller than the acquisition area, computation costs are reduced significantly by the use of a novel bottom-up strategy to calculate the extrapolated Green’s functions. Target-oriented RTD is tested on 2D and 3D SEG/EAGE synthetic data sets and a 3D field data set from the Gulf of Mexico. Results show that target-oriented RTD combined with standard migration can image sediments beneath complex structures accurately with much less calculation effort than full volume RTM. The requirement is that the area over the target zone is smaller than that of the acquisition survey.

Surface waves are a form of coherent noise that can obscure valuable reflection information in exploration records. It is sometimes difficult to eliminate these surface waves by traditional filtering approaches, such as an [Formula: see text] filter, without damaging the useful signals. As a partial remedy, we propose an interferometric method to predict and subtract surface waves in seismic data. The removal of surface waves by the proposed interferometric method consists of three steps: (1) remove most of the surface waves by a nonlinear local filter; (2) predict the residual surface waves by the interferometric method; (3) separate the residual surface waves from the result of step 2 by a nonlinear local filter, and remove the residual surface waves by a matched filter from the result of step 1. Field data tests for 2D and 3D data show that the method effectively suppresses surface waves and preserves the reflection information. Results suggest that the effectiveness of this method is sensitive to the parameter selection of the nonlinear local filter.

Interferometric migration of free-surface multiples in vertical-seismic-profile (VSP) data has two significant advantages over standard VSP imaging: (1) a significantly larger imaging area compared to migrating VSP primaries and (2) less sensitivity to velocity-estimation and static errors than other methods for migration of multiples. In this paper, we present a 3D wave-equation interferometric migration method that efficiently images VSP free-surface multiples. Synthetic and field data results confirm that a reflectivity image volume, comparable in size to a 3D surface seismic survey around the well, can be computed economically. The reflectivity image volume has less fold density and lower signal-to-noise ratio than that obtained by a conventional 3D surface seismic survey because of the relatively weak energy of multiples and the limited number of geophones in the well. However, the efficiency of this method for migrating VSP multiples suggests that it might sometimes be a useful tool for 4D seismic monitoring where reflectivity images can be computed quickly for each time-lapse survey.

An interferometric form of Fermat’s principle and traveltime tomography is used to invert ground-penetrating radar (GPR) data for the subsurface velocity distribution. The input data consist of GPR traveltimes of reflections from two buried interfaces, [Formula: see text] (reference) and [Formula: see text] (target), where the data are excited and recorded by GPR antennas at the surface. Fermat’s interferometric principle is then used to redatum the surface transmitters and receivers to interface [Formula: see text] so the associated reflection traveltimes correspond to localized transit times between interfaces [Formula: see text] and [Formula: see text]. The overburden velocity model above interface [Formula: see text] is not required. The result after tomographic inversion is a high-resolution estimate of the velocity between interfaces [Formula: see text] and [Formula: see text] that does not depend on the velocity model above interface [Formula: see text]. A motivation for introducing interferometric traveltime tomography is that typical layer-stripping approaches will see the slowness error increase with depth as the layers are inverted. This suggests that near-surface statics errors are propagated and amplified with depth. In contrast, the interferometric traveltime tomography method largely eliminates statics errors by taking the difference between reflection events that emanate from neighboring layer interfaces. Slowness errors are not amplified with depth. However, the method is sensitive to the estimation accuracy for the geometry of the reference interface. Both synthetic and real field data are used successfully to validate the effectiveness of this interferometric technique.

One of the difficulties in seeing beneath salt is that the migration velocity in the salt and above it is not well known. This can lead to defocusing of migration images beneath the salt. In this paper, we show that reduced-time migration (RTM) and interferometric migration (IM) can partly mitigate this problem. RTM time-shifts the traces with the time difference between the calculated arrival time [Formula: see text] and the natural arrival time [Formula: see text] of a reference reflection, where [Formula: see text] and [Formula: see text] denote the source and receiver locations on the surface, respectively. We use the terms natural and calculated to represent, respectively, the arrival times that are velocity-independent (traveltimes directly extracted from the data without knowledge of the velocity model) and velocity-dependent (traveltimes calculated by ray tracing through a given velocity model). The benefit of RTM is a significant reduction of defocusing errors caused by errors in the migration velocity. IM, on the other hand, requires extrapolation of the surface data below salt using the natural arrival times [Formula: see text] of the subsalt reference reflector, and migration of the extrapolated data below the salt. The benefit with IM is that no salt velocity model is needed, so the model-based defocusing errors are, in theory, eliminated. To reduce computational time, we implement IM with a seminatural Green’s function (combination of model-based calculated and picked natural traveltimes). Because no explicit data extrapolation is needed, IM with seminatural Green’s functions is more cost-efficient than the standard IM. In this paper, we tested both RTM and IM with seminatural Green’s functions on a synthetic and a field common-depth-point (CDP) data set, the latter from the Gulf of Mexico (GOM). Results show that both RTM and IM can remove the significant kinematic distortions caused by the overburden without knowledge of the overburden velocity.

We develop a waveform-tomography method for estimating the velocity distribution that minimizes the waveform misfit between the predicted and observed early arrivals in space-time seismograms. By fitting the waveforms of early arrivals, early arrival waveform tomography (EWT) naturally takes into account more general wave-propagation effects compared to the high-frequency method of traveltime tomography, meaning that EWT can estimate a wider range of slowness wavenumbers. Another benefit of EWT is more reliable convergence compared to full-waveform tomography, because an early-arrival misfit function contains fewer local minima. Synthetic test results verify that the waveform tomogram is much more accurate than the traveltime tomogram and that this algorithm has good convergence properties. For marine data from the Gulf of Mexico, the statics problem caused by shallow, gassy muds was attacked by using EWT to obtain a more accurate velocity model. Using the waveform tomogram to correct for statics, the stacked section was significantly improved compared to using the normal move-out (NMO) velocity, and moderately improved compared to using the traveltime tomogram. Inverting high-resolution land data from Mapleton, Utah, showed an EWT velocity tomogram that was more consistent with the ground truth (trench log) than the traveltime tomogram. Our results suggest that EWT can provide supplemental, shorter-wavelength information compared to the traveltime tomogram for both shallow and moderately deep seismic data.

We review the equations for correlation-based redatuming methods. A correlation-based redatuming method uses natural-phase information in the data to time shift the weighted traces so they appear to be generated by sources (or recorded by geophones) shifted to a new location. This compares to model-based redatuming, which effectively time shifts the traces using traveltimes computed from a prior velocity model. For wavefield redatuming, the daylight imaging, interferometric imaging, reverse-time acoustics (RTA), and virtual-source methods all require weighted correlation of the traces with one another, followed by summation over all sources (and sometimes receivers). These methods differ from one another by their choice of weights. The least-squares interferometry and virtual-source imaging methods are potentially the most powerful because they account for the limited source and receiver aperture of the recording geometry. Interferometry, on the other hand, has the flexibility to select imaging conditions that target almost any type of event. Stationary-phase principles lead to a Fermat-based redatuming method known as redatuming by a seminatural Green’s function. No crosscorrelation is needed, so it is less expensive than the other methods. Finally, Fermat’s principle can be used to redatum traveltimes.

Prestack migration in the time and depth domains is the premier tool for seismic imaging of complex structures. Unfortunately, an undersampled acquisition geometry, a limited recording aperture, and strong velocity contrasts lead to uneven illumination of the subsurface. This results in blurring the migrated image, sometimes referred to as acquisition footprint noise in the migrated image. To remedy this blurring partially, we introduce prestack migration deconvolution (MD). The MD filter is a layered-medium approximation to the inverse Hessian matrix that is applied locally to the migrated section. Both synthetic- and field-data results show noticeable improvements in reducing migration artifacts and increasing lateral spatial resolution by more than 10%. The computational expense for constructing the MD filter is related to the MD operator length. Its cost is about the same as that for migration, but opportunities exist for significantly reducing this cost. Results suggest that MD should be applied to migrated sections to optimize image quality.

We present the theory of crosscorrelogram migration of ghost reflections, also known as interferometric imaging, to delineate reflector geometries from inverse vertical seismic profile data. The theory includes the equations for forward modeling, migration, asymptotic inversion, and model resolution of crosscorrelgrams. Rather than using primary reflections, crosscorrelogram migration can use ghost reflections to reconstruct the reflector geometry. Other multiples can be used as well, including pegleg multiples and higher-order multiples. Its main advantages over conventional Kirchhoff migration are (1) both source location (e.g., drill-bit position) and source wavelet need not be known (such as when using a drill bit as a source in a deviated well), (2) it is insensitive to source-related static errors, and (3) it has wider subsurface illumination than conventional Kirchhoff migration of primary reflections. Crosscorrelation migration can effectively widen the source-receiver aperture by more than 50% compared to standard inverse vertical seismic profile (IVSP) migration. The primary disadvantages are (1) it uses ghost reflections for imaging, which can be weaker (or sometimes more distorted) than primary reflections; (2) crosscorrelation creates virtual multiples that can sometimes appear as coherent noise in the final image; and (3) it has poorer horizontal resolution than standard IVSP migration. Results from imaging simulated IVSP traces show that crosscorrelation migration produces reflectivity-like images that correlate well with the actual reflector geometry of a layered fault model. These images are almost completely immune to static errors at the source location and have wider subsurface illumination than conventional IVSP migration images. We also apply crosscorrelation migration to IVSP data recorded at a Friendswood, Texas, test site. Results show that the crosscorrelation image correlates better with the well log and wider subsurface illumination than a conventional Kirchhoff migration image.

An interferometric form of Fermat's principle is derived that allows for high-resolution estimation of the velocity distribution between deep interfaces. The data consist of reflection traveltimes from two deeply buried interfaces A and B recorded by sources and caused by receivers at the surface. Fermat's interferometric principle is then used to kinematically redatum the surface sources and receivers to interface A so that the associated reflection times correspond to localized transit times between the A and B interfaces. The velocity model above interface A does not need to be known, so the distorting effects of the overburden and statics are eliminated by this target-oriented approach. Interferometric target-oriented tomography promises to enhance the resolution of whole-earth and exploration tomograms.

We develop the novel theory of transmitted PS migration and show that PS transmitted arrivals in a Gulf of Mexico vertical seismic profile (VSP) data set can be migrated to accurately image a salt sheet even though the receiver array is below the transmitting boundary. We also show that migrating transmitted arrivals is effective in illuminating the base of an orebody invisible to PP reflections. In general, interfaces that bisect wavepath propagation (i.e., the source and receiver are on opposite sides of the interface and therefore invisible to PP reflections) can be imaged by migration of PS transmitted waves. These results suggest that migration of PS transmitted waves opens new opportunities in imaging nearly vertical impedance boundaries that are typically invisible to conventional reflection imaging of crosswell and VSP data.We also present a new interferometric method, denoted as reduced‐time migration, which uses the arrival‐time difference between the direct P‐wave and subsequent events to increase migration accuracy. Reduced‐time migration removes static time shifts in the data, decreases the focusing error due to an incorrect migration velocity model, and relocates reflection or PS transmission events to be closer to their true positions. Although limited to crosswell and VSP geometries, synthetic‐ and field‐data examples show that reduced‐time migration is noticeably more accurate than conventional migration in the presence of static shifts and/or migration velocity errors. The main assumption of reduced‐time migration is that the direct wave samples errors which are representative of errors in the migration aperture. Transmission wavepaths, in general, are subparallel to the direct wave and therefore the two modes encounter similar errors and, hence, reduced‐time migration is effective in improving the focusing of migration energy. For the PP reflection case, the direct wave and the reflected waves often traverse different parts of the earth, therefore, reduced‐time migration will remove static shifts but it is not expected to mitigate velocity errors if the errors are spatially variant. However, if there is a general and consistent bias in the velocity model, reduced‐time migration is expected to deliver improved results over conventional Kirchhoff migration.

We present the equations for migrating inverse‐vertical‐seismic‐profile‐while‐drilling and common‐midpoint autocorrelograms. These equations partly generalize the 1D autocorrelation imaging methods of Katz and Claerbout to 2D and 3D media, and also provide a formal mathematical procedure for imaging the reflectivity distribution from autocorrelograms. The imaging conditions are designed to migrate specific events in the autocorrelograms, either the direct‐primary correlations or the direct‐ghost correlations. Here, direct stands for direct wave, primary stands for primary reflections, and ghost denotes free‐surface ghost reflections. The main advantage in migrating autocorrelograms is that the source wavelet does not need to be known, which is the case for seismic data generated by a rotating drill bit or for vibroseis data with a corrupted pilot signal. Another advantage is that the source and receiver static problems are mitigated by autocorrelation migration. Two limitations are that autocorrelation of traces amplifies coherent noise such as surface waves, and produces undesirable coherent noise denoted as “virtual multiples.” Similar to “physical multiples,” such noise can, in principle, be partially suppressed by filtering and stacking of migration images obtained from many different shot gathers. Results with both synthetic and field data validate this conjecture, and show that autocorrelogram migration can be a viable alternative to standard migration when the source signal is not adequately known or there are severe static problems.

The acquisition footprint noise in migrated sections partly consists of migration artifacts associated with a discrete recording geometry. Such noise can corrupt the interpretation of seismic sections for stratigraphic variations, AVO signatures, and enhanced oil recovery operations. We show that the point scatterer response of the far‐field Kirchhoff migration operator, which reveals acquisition footprint noise, is proportional to the stretched Fourier transform of the source and geophone sampling function. Using the array theorem developed by optical and electrical engineers, this Fourier transform for an orthogonal recording geometry can be calculated quickly by a product of 1‐D analytical functions. Thus, the acquisition footprint noise in migrated sections can be calculated efficiently for different recording arrays. By rapid trial and error or by an optimization method, the survey planner can use this theorem to help design a better recording geometry.

Prestack Kirchhoff migration (KM) is computationally intensive for iterative velocity analysis. This is partly because each time sample in a trace must be smeared along a quasi‐ellipsoid in the model. As a less costly alternative, we use the stationary phase approximation to the KM integral so that the time sample is smeared along a small Fresnel zone portion of the quasi‐ellipsoid. This is equivalent to smearing the time samples in a trace over a 1.5‐D fat ray (i.e., wavepath), so we call this “wavepath migration” (WM). This compares to standard KM, which smears the energy in a trace along a 3‐D volume of quasi‐concentric ellipsoids. In principle, single trace migration with WM has a computational count of [Formula: see text] compared to KM, which has a computational count of [Formula: see text], where N is the number of grid points along one side of a cubic velocity model. Our results with poststack data show that WM produces an image that in some places contains fewer migration artifacts and is about as well resolved as the KM image. For a 2‐D poststack migration example, the computation time of WM is less than one‐third that of KM. Our results with prestack data show that WM images contain fewer migration artifacts and can define the complex structure more accurately. It is also shown that WM can be significantly faster than KM if a slant stack technique is used in the migration. The drawback with WM is that it is sometimes less robust than KM because of its sensitivity to errors in estimating the incidence angles of the reflections.

Reflectivity images of the earth are calculated by migrating discrete grids of seismic traces. Typically, such traces are spatially undersampled on a recording grid with limited aperture width and so give rise to migration noise sometimes referred to as the acquisition footprint. For poststack migration images, we show how to partly deconvolve the acquisition footprint by applying a deblurring filter to the migration section, where the filter is the approximate inverse to the migration Green’s function. Results with synthetic and field data show that post‐stack migration deconvolution can noticeably improve the spatial resolution of migration images, decrease the strength of migration artifacts, and improve the quality of the migration image. We conclude that migration deconvolution can be a viable alternative to some of the other postmigration processing procedures based on statistics and ad hoc parameter choices.

A key issue in wavefield separation is to find a domain where the signal and coherent noise are well separated from one another. A new wavefield separation algorithm, called migration filtering, separates data arrivals according to their path of propagation and their actual moveout characteristics. This is accomplished by using forward modeling operators to compute the signal and the coherent noise arrivals. A linearized least‐squares inversion scheme yields model estimates for both components; the predicted signal component is constructed by forward modeling the signal model estimate. Synthetic and field data examples demonstrate that migration filtering improves separation of P-wave reflections and surface waves, P-wave reflections and tube waves, P-wave diffractions, and S-wave diffractions. The main benefits of the migration filtering method compared to conventional filtering methods are better wavefield separation capability, the capability of mixing any two conventional transforms for wavefield separation under a general inversion framework, and the capability of mitigating the signal and coherent noise crosstalk by using regularization. The limitations of the method may include more than an order of magnitude increase in computation costs compared to conventional transforms and the difficulty of selecting the proper modeling operators for some wave modes.

Analytical formulas are derived to give the point scatterer responses for the zero‐offset and nonzero‐offset migration operators. These formulas are derived under the far‐field approximation for a continuous distribution of sources and receivers, and are denoted as migration Green’s functions. By specifying a homogeneous velocity model with irregular layers, these Green’s functions can be used to quickly compute the synthetic migrated sections without having to compute the modeled data. Thus, these formulas provide a fundamental understanding of migration artifacts and their relationship to the image point location and source‐receiver configuration. This understanding lays the foundation for the next step, elimination of these artifacts by migration deconvolution.

The quasi‐Monte Carlo migration algorithm is applied to a 3-D seismic data set from West Texas. The field data were finely sampled at approximately 220-ft (67-m) intervals in the in‐line direction but were sampled coarsely at approximately 1320-ft (402-m) intervals in the cross‐line direction. The traces at the quasi‐Monte Carlo points were obtained by an interpolation of the regularly sampled traces. The subsampled traces at the quasi‐Monte Carlo points were migrated, and the resulting images were compared to those obtained by migrating both regular and uniform grids of traces. Results show that, consistent with theory, the quasi‐Monte Carlo migration images contain fewer migration aliasing artifacts than the regular or uniform grid images. For these data, quasi‐Monte Carlo migration apparently requires fewer than half the number of the traces needed by regular‐grid or uniform‐grid migration to give images of comparable quality. These results agree with related migration tests on synthetic data computed for point scatterer models. Our results suggest that better migration images might result from data recorded on a coarse quasi‐random grid compared to regular or uniform coarse grids.

A least‐squares migration algorithm is presented that reduces the migration artifacts (i.e., recording footprint noise) arising from incomplete data. Instead of migrating data with the adjoint of the forward modeling operator, the normal equations are inverted by using a preconditioned linear conjugate gradient scheme that employs regularization. The modeling operator is constructed from an asymptotic acoustic integral equation, and its adjoint is the Kirchhoff migration operator. We tested the performance of the least‐squares migration on synthetic and field data in the cases of limited recording aperture, coarse sampling, and acquisition gaps in the data. Numerical results show that the least‐squares migrated sections are typically more focused than are the corresponding Kirchhoff migrated sections and their reflectivity frequency distributions are closer to those of the true model frequency distribution. Regularization helps attenuate migration artifacts and provides a sharper, better frequency distribution of estimated reflectivity. The least‐squares migrated sections can be used to predict the missing data traces and interpolate and extrapolate them according to the governing modeling equations. Several field data examples are presented. A ground‐penetrating radar data example demonstrates the suppression of the recording footprint noise due to a limited aperture, a large gap, and an undersampled receiver line. In addition, better fault resolution was achieved after applying least‐squares migration to a poststack marine data set. And a reverse vertical seismic profiling example shows that the recording footprint noise due to a coarse receiver interval can be suppressed by least‐squares migration.