Novel wave-equation imaging
Plane-wave migration in tilted coordinates
Plane-wave migration in tilted coordinates is powerful to image steeply dipping reflectors, although the one-way wave-equation operator is used. In plane-wave migration, the recorded surface data are transferred to plane-wave data by slant-stack processing. Both the source plane-wave and the corresponding slant-stacking data are extrapolated into the subsurface, and images are generated by cross-correlating these two wavefields. For each plane-wave source, we assign tilted coordinates whose direction depends on the propagation direction of the plane-wave. For isotropic media, the one-way wave-equation operator does not change. For vertical transversely isotropic (VTI) media, one-way wave-equation operators for tilted transversely isotropic (TI) media are required because of the rotation of the coordinates. I apply this method to the BP 2004 velocity benchmark, to a synthetic dataset for VTI media, and to a real anisotropic dataset. The numerical examples show that plane-wave migration in tilted coordinates can image steeply dipping reflectors and overturned waves, and it is a good tool for salt delineation.
Generalized Riemannian wavefield extrapolation
This paper extends wavefield extrapolation to generalized Riemannian spaces. The key component is the development of a dispersion relationship appropriate for propagating wavefield on generalized non-orthogonal meshes. This wavenumber contains a number of mixed-domain fields in addition to velocity that represent coordinate system geometry. An extended split-step Fourier approximation of the extrapolation wavenumber is developed, which provides accurate results when multiple reference parameters sets are used. Three examples are presented that demonstrate the validity of the theory. An important consequence is that greater emphasis can be placed on generating smoother computational meshes rather than satisfying restrictive semi-orthogonal criteria. This result should lead to more accurate and efficient generalized Riemannian wavefield extrapolation.
Mesh generation using differential methods
This paper examines a differential gridding method for generating computational meshes appropriate for solving partial differential equations. Differential methods pose mesh generation as an elliptical boundary value problem within a framework of differential geometry. Generalized Laplacian operators are used to propagate the known coordinate values on the boundary points into the interior in a smooth manner. The methodology allows for the specification of monitor functions that provide mesh regularization and prevent grid clustering. Examples are provided for two seismic imaging applications: wave-equation Green's function generation and wave-equation migration from topography. In both cases, the resulting regularized meshes have minimal convexity and are conformal to the the prescribed boundaries.
Prestack exploding-reflectors modeling for migration velocity analysis
I introduce a prestack generalization of exploding-reflectors modeling that has the potential to drastically reduce the computational cost of Migration Velocity Analysis (MVA) based on wavefield-continuation migration and modeling. The method aims at synthesizing, starting form a prestack migrated image, a limited number of independent experiments that contain all the velocity information needed for MVA. The basic element of the method is the modeling of one isolated Subsurface-Offset Domain Common Image Gather (SODCIG) by using the prestack image as an initial condition. Taking advantage of the limited range of the subsurface offsets in the migrated image, we can combine many SODCIGs in the same modeling experiment without diminishing the velocity information contained in the data. The number of independent experiments required depends on how close we are to the correct migration velocity. Theoretically, this number is the same as the number of subsurface offsets needed to represent the prestack image from which we started the modeling. Several numerical examples illustrate the basic concepts of the proposed method, and demonstrate its potential usefulness for MVA.
Inverting wave-equation operators
AMO inversion to a common azimuth dataset: Field data results
I cast 3-D data regularization as a least-squares inversion problem. I form a linear operator that maps from irregular dataspace 5-D data space to a regular 4-D common azimuth volume using a cascade of linear interpolation and Azimuth Move-out (AMO) binning. I regularize the inversion by adding minimizes the difference between various () cubes by applying a filter that acts along offset. AMO is used to transform the cubes to the same hx before applying the filter. Further efficiency is gained by inverting each frequency independently. I apply the methodology on marine dataset and compare the results with two more conventional approaches.
PS - Azimuth moveout: Real data application
We present an application of the PS-AMO operator for partial-stacking 3-D multicomponent, ocean-bottom seismic data. The partial-stacking problem is defined as a regularized least-squares objective function. To preserve the resolution of dipping events, the regularization term uses the PS-AMO operator. Application of this methodology on a portion of the Alba 3-D, multicomponent, ocean-bottom seismic data set shows that we can satisfactorily obtain an interpolated data set that honors the physics of converted waves.
Target-oriented wave-equation inversion with regularization in the subsurface-offset domain
A complex velocity model can cause shadow zones in an image computed by migration due to poor illumination. These shadow zones may contain weak signals masked by artifacts. To reduce artifacts and recover the real signal, a wave-equation target-oriented inversion scheme can be developed that uses an explicitly computed least squares inversion Hessian. To solve the otherwise ill-conditioned inversion problem a model regularization needs to be added. One choice for regularization is to penalize the energy in the image not focused at zero subsurface-offset. The subsurface-offset Hessian needs to be computed by using the adjoint of migration as the modeling operator. Results on Sigsbee model show encouraging results.
Computer aided interpretation
Parallel implementation of image segmentation for tracking 3D salt boundaries
We distribute the modified normalized cuts image segmentation with random boundaries algorithm on a parallel network to track 3D salt boundaries. We identify two key steps of this algorithm for parallelization. Firstly, we parallelize the calculation of the weight matrix. Secondly, we parallelize the matrix-vector product of the eigenvector calculation. This method is demonstrated to be effective on a 3D seismic cube.
lattening with geological constraints
In areas with faults and poor signal/noise ratio, where reflectors can be discontinuous from place to place, a dip-based flattening technique might not be able to appropriately track sedimentary layers. To aid the flattening algorithms, one or few reflectors can be picked. This information can be then incorporated in our algorithms as geological constraints. In a first method, we add a model mask to a time domain solution using a Gauss-Newton approach that incorporates an initial solution. In a second method, we set the lower and upper bounds of a constrained optimization algorithm called limited memory BFGS with bounds (L-BFGS-B). Having incorporated the geological information, the flattening algorithms can accurately pick reflectors in 2D and 3D for noisy field data examples. In addition, preliminary results seem to indicate that the L-BFGS-B method converges faster than the Gauss-Newton method.
Quaternion-based Signal Processing
Hypercomlex numbers, which have primarily been used for pattern recognition, offer many useful applications to geophysics. Image disparity estimation is a hypercomplex, phase-based technique, using quaternions, that can find differences between subtlety varying images. This technique relies on applying a quaternionic Fourier transform, a quaternionic Gabor filter and exploits the symmetries inherent in the quaternion. Two applications of hypercomplex image disparity estimation are time lapse analysis and boundary detection.
Flattening with cosine transforms
In the Fourier-domain flattening methods presented previously Lomask et al. (2005); Lomask and Claerbout (2002); Lomask (2003) the data has to be mirrored in order to eliminate Fourier artifacts. This means that the data is replicated and reversed so that the boundaries are periodic. This requires four times the memory in 2D and eight times in 3D. In a world where post-stack data cubes can easily be tens of gigabytes in size, efficient memory usage is extremely important. Here we apply a discrete cosine transform (DCT) approach developed for 2D phase unwrapping Ghiglia and Pritt (1998); Ghiglia and Romero (1994) to our flattening method. As a result, we the reduce memory requirements of the transforms in 2D by a factor of four and in 3D by a factor of eight. We reap an additional factor of two savings from using real instead of complex numbers. Furthermore, the reduction in size significantly reduces computation time as well. We demonstrate its use on a simple 3D synthetic model.
Wave-equation angle-domain common-image gathers for converted waves: Part 2
Common-image gathers are very useful for velocity and amplitude analysis. Wavefield-extrapolation methods produce Angle-Domain Common-Image Gathers (ADCIGs). For the conventional PP case, ADCIGs are a function of the opening angle. However, the ADCIGs for converted-wave data (PS-ADCIGs) are a function of the half-aperture angle, that is, the incidence angle plus the reflection angle. In PS-ADCIGs, both the P-to-S velocity ratio () and the image dip play a major role in transforming the subsurface offset into the opening angle. We introduce a simple methodology to compute PS-ADCIGs. Our methodology exploits the robustness of computing PP-ADCIGs, and incorporates the velocity ratio (), with an image dip field, which is estimated along the prestack image. Our methodology also transforms the half-aperture angle in PS-ADCIGs into an independent P-incidence angle to form P-ADCIGs, and an independent S-reflection angle to form S-ADCIGs. Numerical studies show that when the P-to-S velocity ratio and image midpoint information are not incorporated, the error in computing PS-ADCIGs is large enough to introduce artifacts in the velocity model. Synthetic results show the accuracy of the transformation introduced in this paper. Real data results on the 2-D Mahogany field show the practical application and implications for converted-wave angle-domain common-image gathers.
Residual moveout of 2D multiples in Angle-Domain Common-Image Gathers
I show that, for specularly-reflected multiples, the constant velocity straight-ray approximation of the residual moveout in Angle-Domain Common-Image Gathers (ADCIGs) is only appropriate for small aperture angles. The approximation is good for the primaries because the difference between the migration velocity and the true velocity is likely to be small. For the multiples, however, this difference may be large and correcting for ray bending produces a better approximation that leads to better focusing of the multiples in the Radon domain. This in turn allows a more accurate muting of the multiples. I show results with two ADCIGs, one synthetic and one real.
Residual moveout in anisotropic angle-domain common image gathers with dipping reflectors
We generalise to dipping reflectors the fundamental concepts for quantitatively relating perturbations in anisotropic parameters to the corresponding reflector movements in angle-domain common-image gathers (ADCIGs). We apply the general methodology to the particular case of residual moveout (RMO) analysis of reflections from dipping reflectors in a vertical transverse isotropic (VTI) medium. Synthetic examples show the accuracy of the RMO curves predicted by our kinematic analysis.
Attenuation of 2D specularly-reflected multiples in image space
I propose a new method to attenuate specularly-reflected multiples in the image space. The method is based on the difference in mapping between primaries and multiples in Subsurface Offset Domain Common Image Gathers (SODCIGs). I migrate the data with a velocity slower than that of the primaries but faster than that of the multiples. The primaries are therefore undermigrated whereas the multiples are overmigrated. For positive data offsets, the primaries are mapped to positive subsurface offsets and the multiples to negative subsurface offsets in SODCIGs. I then apply a tapered mute to eliminate the primaries and do adjoint migration on the multiples with the same velocity model to get an estimate of the multiples in data space. Similarly, a tapered mute is applied to eliminate the multiples and adjoint migration used to obtain an estimate of the primaries in data space. The estimate of the multiples is adaptively matched to the data with the estimate of the primaries used as a weight function to prevent matching the primaries. I illustrate the method with a 2D synthetic dataset and show that the primaries can be well recovered although some residual from the water bottom multiple remains.
Mapping of specularly-reflected multiples to image space: An example with 3D synthetic data
In 2D, specularly-reflected multiples, when migrated with the velocity of the primaries, map to negative subsurface offsets in Subsurface-Offset-Domain Common-Image Gathers (SODCIGs). In Angle-Domain Common-Image Gathers (ADCIGs) they map with curvature towards increasing depths. Here we show, through a 3D synthetic prestack dataset, that specularly-reflected multiples in 3D have a similar behavior with an interesting addition: in 3D ADCIGs, the multiples exhibit an azimuth rotation proportional to the dip of the reflecting interface generating the multiple. This attribute may be used to discriminate between primaries and multiples in 3D ADCIGs and therefore help in the attenuation of the multiples.
Application of Least-Squares Joint Imaging of Multiples and Primaries on Shallow Water-Bottom Data Sets
Data contaminated with strong shallow water-bottom multiples is rife with challenges. Application of Least-Squares Joint Imaging of Multiples and Primaries (LSJIMP) on such data sets yields mixed results. LSJIMP solves both the separation and integration simultaneously, as a global least-squares inverse problem. We point out some limitations of LSJIMP by testing it on synthetic data sets that emulate shallow water-bottom marine environments. Some slight modifications have been made, and we suggest some strategies that might make LSJIMP an effective algorithm.
Interpolation with pseudo-primaries
Large gaps exist in marine data, particularly at near offsets. I generate pseudo-primaries by cross-correlating a multiple model with the original data. These pseudo-primaries are used as training data for a non-stationary prediction- error filter, which is then used to interpolate the missing near offsets. This method yields good results, and also provides a quality control measure to judge the usefulness of the pseudo-primaries.
Interpolating diffracted multiples with prediction-error filters
Diffracted multiples are a persistent problem, especially in the cross-line direction. Methods to remove these multiples typically require dense and complete sampling. I use non-stationary prediction-error filters to interpolate extra shots in the in-line direction and extra receivers in the cross-line direction. The dimensionality of the interpolation does not make a big difference in the in-line direction, but makes a large difference in the cross-line direction. This is shown with several tests on a synthetic dataset.
Missing Data Interpolation with Gaussian Pyramids
I describe a technique for interpolation of missing data in which local operators of many scales but identical shape serve as basis functions. A data structure known as the Gaussian pyramid is developed to represent image information at different scales. This data structure in essence consists of a series of lowpass filtered versions of the original image stacked up one on top of the other forming a pyramid like structure. I first show how to generate a set of reduced images which stack up to form the Gaussian pyramid structure and then show how we can use this Gaussian pyramid structure to fill in missing data. Several examples of filling in missing data with this algorithm are shown and in most cases the results are comparable with those estimated using a prediction filter approach.
A modified Lloyd algorithm for characterizing vector fields
Lloyd's 1-D method attempts to identify the major features in a signal's histogram as accurately as possible, with as few points as possible. The selected points are referred to as the signal's ``codebook''. The methodology is iterative in nature. It starts from an initial code book, then repeats .
Selection of reference anisotropic parameters for wavefield extrapolation by Lloyd's algorithm
We propose a method for selecting reference anisotropic parameters in laterally varying anisotropic media for mixed Fourier-space domain wavefield extrapolation. We treat the selection problem as a quantization procedure, and use a modified version of the 3D Lloyd's algoritm for reference-parameter selections. We demonstrate that our method yields a more accurate discription of the anisotropy model with fewer reference parameters than the uniform sampling approach. Real data examples illustrate the performance of our method.
Inversion shortcuts by model statistics
A set of model-space coordinates must be pre-selected for input into LA according to the amplitude (squared) at all coordinate locations. Only coordinates with an amplitude greater than a supplied percentage of the maximum amplitude in the model space are kept. The amplitude range of model space is then quantized from the minimum threshold to the maximum (squared) value. Coordinates are repeatedly selected for input into LA according to the number of quantum levels associated with its amplitude. Figure shows graphically how the data input into LA is .
Optimized implicit finite-difference migration for VTI media
I develop an implicit finite-difference migration method for vertical transversely isotropic (VTI) media with laterally varying anisotropy parameters. I approximate the dispersion relation of VTI media with a rational function series, the coefficients of which are estimated by least-squares optimization. These coefficients are functions of Thomsen anisotropy parameters. They are calculated and stored in a table before the wavefield extrapolation. The implicit finite-difference scheme for VTI media is almost the same as that of the isotropic media, except that the coefficients are derived from the pre-calculated table. In the 3D case, a phase-correction filter is applied after the finite-difference operator to eliminate the numerical-anisotropy error caused by two-way splitting. This finite-difference operator for VTI media is accurate to , and its computational cost is almost the same as the isotropic migration. I apply this method to a 2D synthetic dataset and a 2D slice of a real 3D dataset to validate the method.
Seismic waves in rocks with fluids and fractures
Seismic wave propagation through the earth is often strongly affected by the presence of fractures. When these fractures are filled with fluids (oil, gas, water, CO2, etc.), the type and state of the fluid (liquid or gas) can make a large difference in the response of the seismic waves. This paper will summarize some early work of the author on methods of deconstructing the effects of fractures, and any fluids within these fractures, on seismic wave propagation as observed in reflection seismic data. Methods to be explored here include Thomsen's anisotropy parameters for wave moveout (since fractures often induce elastic anisotropy), and some very convenient fracture parameters introduced by Sayers and Kachanov that permit a relatively simple deconstruction of the elastic behavior in terms of fracture parameters (whenever this is appropriate).
- Ghiglia, D. C. and M. D. Pritt, 1998, Two-dimentional phase unwrapping: Theory, algorithms, and software: John Wiley and Sons, Inc.
- Ghiglia, D. C. and L. A. Romero, 1994, Robust two-dimensional weighted and unweighted phase unwrapping that uses fast transforms and iterative methods: Optical Society of America, 11, no. 1, 107-117.
- Lomask, J. and J. Claerbout, 2002, Flattening without picking: SEP-112, 141-150.
- Lomask, J. and A. Guitton, 2006, Flattening with geological constraints: SEP-124.
- Lomask, J., A. Guitton, S. Fomel, and J. Claerbout, 2005, Update on flattening without picking: SEP-120, 137-158.
- Lomask, J., 2003, Flattening 3-D data cubes in complex geology: SEP-113, 247-260.