Decomposition of multicomponent seismic data into primary compressional (P) and shear (S) wave responses and migration of these decomposed wavefields are two main procedures for the multicomponent data processing. In this work, first we derive the decomposition formula from elastic wave equation to separate out the displacement vector with coupled P- and SV-waves into scalar upgoing P- and SV-waves. Pure P-wave can be decomposed by revising the vertical component with the correction of horizontal component, which can enhance P- wave but eliminate SV-wave in vertical component. Similarly, pure SV-wave can be decomposed by revising the horizontal component with the use of vertical component. Such a decomposition process requires only the knowledge of elastic parameters near the surface. Two kinds of migration schemes are then proposed to reconstruct the subsurface image by performing downward extrapolation of the separated scalar upgoing P- and SV-wave fields. For common shot gathers, the P-SV subsurface can be imaged with the correlation of forward extrapolation from the source for downgoing waves and inverse extrapolation for the upgoing P- and SV-waves with the respective P- and SV-wave velocity. The screen propagators can be used effectively in this case. For common offset sections, we modify the travel time equation of P-SV point scatterer and derive the corresponding migration equation. Results from synthetic data of one two-component record show that the P- and SV-waves are decomposed completely and the subsurface image for P-P waves is consistent with that of P-SV waves.

Two efficient Fourier migration methods termed the extended local Born Fourier (ELBF) method and the extended Rytov Fourier (ELRF) method have been developed recently for imaging complex 3D structures. They are recursive methods based on local applications of Born and Rytov approximations within each extrapolation interval. The ELBF method becomes unreliable when the lateral slowness variations are large and/or the frequency is high, while the ELRF method is reliable for such cases. However, the ELRF method is approximately 30-40% slower than the ELBF method because the ELRF method requires one more computational step where exponentials of complex numbers are calculated than the ELBF method and propose an implementation scheme using variable extrapolation intervals to make the ELBF method reliable for all lateral slowness variations and frequencies. The size of the extrapolation interval depends on the lateral slowness variations within a given extrapolation region and the frequency, and consequently, the computational time of the ELBF method with variable extrapolation intervals increases with the lateral slowness variation and frequency. To take advantage of the faster computational speed of the ELBF method compared to the ELRF method and the better stability of the ELRF method compared to the ELBF method, we propose a hybrid local Born/Rytov Fourier migration method. In the hybrid method, the ELBF method is used for regions with small lateral slowness variations and/or low frequencies, otherwise, the ELRF method is used. Migrations of two synthetic datasets for complex structures using the ELBF method with variable extrapolation intervals and the hybrid method demonstrate that the quality of images obtained using these two methods is comparable to that of images obtained using the ELRF method. Comparison of computational times for migrations using different methods shows that the ELBF method with variable extrapolation intervals takes much more computational time than the ELRF method but the hybrid method saves more than 10% of the computational time required by the ELRF method.

We obtain accurate images of complex subsalt structures by combining Kirchhoff datuming and Kirchhoff migration into a 3D semirecursive Kirchhoff migration. By datuming to the top of salt, or even through the salt, and then imaging below the salt, a greatly improved image is obtained.

The offset plane wave equation, which is obtained from a Radon transformation of the scalar wave equation in midpoint-offset coordinates, provides a convenient framework for developing efficient amplitude preserving migrations. The Radon transformation results in additional amplitude terms that must be accounted for in downward continuation and migration. These terms include an obliquity factor and point source correction resulting from the Radon transform over offset, and transmission terms to account for the transport equation in variable velocity media. These terms can be expressed in exact form in the spectral domain for vertically varying media. For laterally varying media, we currently use a phase-shift plus interpolation framework for the amplitude terms, and a Taylor's series expansion of the velocity field over offset to provide an approximate solution to the offset plane wave equation. Phase shift and finite difference migrations based on this framework produce images on standard models that are comparable to prestack Kirchhoff and shot record migrations.

This paper presents a method for constructing coupled, 1D electrical conductivity models of the Earth from surface magnetotelluric measurements. The construction of the individual models is a nonlinear inverse problem that can be approached by linearization techniques combined with iterative methods and Tikhonov's regularization. The standard application of these techniques usually leads to smooth models that represent a continuous variation of conductivity with depth. In a previous work we described how these methods can be modified to incorporate what is known in Computer Vision as the line process (LP) decoupling technique, which has the ability to include discontinuities in the models. This results in piecewise smooth models which are often more adequate for representing stratified media. Now, the 1D models are coupled in order to include information from neighboring sites. We have implemented a relaxation technique to construct 2D profiles using coupled 1D models. We present numerical experiments and application to field data. In both cases we assume that the data is contaminated by static shift effects. The algorithm automatically takes these effects into account. The examples illustrate the performance of the combined LP and Tikhonov's regularization method in the solution of difficult practical problems.

Phenomenological equations for the poroelastic behavior of a double porosity medium have been formulated and the coefficients in these linear equations identified. The generalization from a single porosity model increases the number of independent coefficients from three to six for an isotropic applied stress. In a quasistatic analysis, the physical interpretations are based upon considerations of extremes in both spatial and temporal scales. The limit of very short times is the one most relevant for wave propagation, and in this case both matrix porosity and fractures are expected to behave in an undrained fashion, although our analysis makes no assumptions in this regard. For the very long times more relevant for reservoir drawdown, the double porosity medium behaves as an equivalent single porosity medium. At the macroscopic spatial level, the pertinent parameters (such as the total compressibility) may be determined by appropriate field tests. At the mesoscopic scale pertinent parameters of the rock matrix can be determined directly through laboratory measurements on core, and the compressibility can be measured for a single fracture. We show explicitly how to generalize the quasistatic results to incorporate wave propagation effects and how effects that are usually attributed to squirt flow under partially saturated conditions can be explained alternatively in terms of the double-porosity model. The result is therefore a theory that generalizes, but is completely consistent with, Biot's theory of poroelasticity and is valid for analysis of elastic wave data from highly fractured reservoirs.

The determination of a correct subsurface velocity model is of strategic importance for seismic imaging of complex geologies by use of 3D prestack depth migration. To accomplish this task we have developed a 3D reflection tomography software (called jerry) that is designed to handle methods of realistic complexity and that is fast enough to allow practical 3D applications. We use a blocky representation of the subsurface with a B-Spline parametrization for interface geometries and for velocity distributions within layers. The forward problem of reflection tomography is formulated as a two-point ray tracing problem and is solved by a bending method based on Fermat's principle with an approximation of the actual ray geometry within a layer by a straight line. The C² property of the model parametrization allows to calculate exact derivatives of the traveltime function with respect to the model and thus to calculate rays (trajectories that satisfy Fermat's principle) from initial trial trajectories. The trajectory initialization is based on Snell's law for transmission/reflection at interfaces and uses an iterative method to calculate the intersection point between straight rays and interfaces. The introduction of non-physical ghost interfaces allows to improve the ray path approximation and thus the accuracy of the calculated traveltimes.

We present a fast algorithm for solving the eikonal equation in three dimensions, based on the Fast Marching Method (FMM). The algorithm is of order O(N log N), where N is the total number of grid points in the computational domain. The algorithm can be used in any orthogonal coordinate system, and globally constructs the solution to the eikonal equation for each point in the coordinate domain. The method is unconditionally stable, and constructs solutions consistent with the exact solution for arbitrarily large gradient jumps in velocity. In addition, the method resolves any overturning propagation wavefronts. We begin with the mathematical foundation for solving the eikonal equation using the FMM, and follow with the numerical details. We show examples of traveltime propagation through the SEG/EAGE Salt Model, and the use of these first arrival traveltimes to image 3D prestack data.

In this paper a simple approach to the quantitative estimate of the limits of applicability of the physical optics method is proposed. In the case of non-absorbing body, the local and integral characteristics of scattered field are compared, i.e. cross-section of extinction and cross-section of scattering which are obtained by means of the physical optics method. It is shown that for arbitrary dimensions of a scatterer cross-section these characteristics are not equal, as it takes place in the case of a rigorous method, but are connected by an inequality. A relative difference of characteristics of extinction and scattering is assigned to a relative error of the physical optics method, within the scope of which they are obtained. For the relative error of the physical optics method the bounding inequalities are derived. It is shown that the relative error for the scatterers of round or square cross-section is equal to 1/p, where p equals ka is the diffraction parameter, a is the radius of a round cross-section or a halfwidth of square cross- section side. The relative error for the scatterers of rectangular cross-section is expressed via the combination of Bessel's functions.

Internal multiples are more difficult to estimate and attenuate than free-surface multiples. In this paper, our purpose is to illustrate how the inverse scattering internal multiple algorithm is able to attenuate internal multiples without any knowledge of the reflectors that cause the multiples or the medium through which the waves propagate. An analytic example is used to demonstrate the inner workings of the algorithm; and, further understanding is attained though an interpretation of the algorithm in terms of subevents. In this way, we gain a better comprehension of the inverse scattering internal multiple algorithm. To date, this is the only method capable of attenuating internal multiples in a multi-dimensional earth without any information about the medium below the sources and receivers.

Reflectivity images of the earth are calculated by migrating discrete grids of seismic traces. For 3D data, these traces, typically, are spatially undersampled and so give rise to migration noise sometimes known as the `acquisition footprint'. In this paper, we show how to partly deconvolve the acquisition footprint noise from migration sections. This is accomplished by spectrally dividing the migration section by the migration Green's function, subject to the assumption of a laterally homogeneous velocity. Results with synthetic and field data show that migration deconvolution can be noticeably improve the spatial resolution of migration images, decrease the strength of migration artifacts and improve the quality of the migration image.

Time-lapse (4D) seismic is a procedure where a reservoir is imaged with reflected seismic energy at several time steps while being depleted. This technology is successful when changes in dynamic reservoir properties, such as pressure, saturation of fluids and temperature produce an observable change in the seismic impedance contrast of the medium in 3D space. A careful analysis combining rock physics measurements, reservoir simulation and forward seismic modeling is required to identify if changes in reservoir properties can produce seismically observable impedance changes, and if so, how they are related to dynamic fluid flow in the reservoir. In current industry practice, angle- averaged (or stacked) P-to-P wave reflection seismic images, approximating P-wave impedance changes, are used in time- lapse seismic analysis. Introduction of angle dependency in P-to-P wave reflection data allows us to estimate both P- wave and S-wave impedance changes over the reservoir. Utilizing multi-parameter images from angle-dependent elastic time-lapse seismic surveys, or 4D AVO (amplitude variation with offset) attributes, we demonstrate the potential to invert for multiple reservoir properties such as changes in saturation of fluids and pressure. Estimates of saturation and pressure changes allow integration of time-lapse seismic with reservoir engineering models to improve reservoir fluid-flow prediction and enhance reservoir management decisions.

As new methods of interpreting 3D seismic data, particularly prestack and derived attribute data, increase in popularity, the management of ever-larger data volumes becomes critical. Compared with acquisition and processing, however, the interpretation use of seismic data requires faster and non- sequential, random access to large data volumes. In addition, quantitative interpretations lead to an increasing need for full 32-bit resolution of amplitudes, rather than the 8 or 16 bit representations that have been used in most interpretation systems up to the present. Seismic data compression can be a significant tool in managing these on- line datasets, but the implementation previously used (e.g. Donoho, et al. 1995) are not well suited to providing rapid random access in arbitrary directions (CDP, crossline, timeslice). If compression ratios of twenty or greater can be routinely provided, with full random access to all parts of the dataset, then both larger cubes and more cubes can be handled within the finite memory and disk systems available on interpretation systems. Such on-line data accessibility will lead to higher productivity by interpreters and greater value for existing seismic surveys. We address here several problems that arise when using the lossy wavelet-transform data compression algorithms currently available. We demonstrate that wavelet compression introduces less noise than currently accepted truncation compression. We also show how compressing small blocks of data needed for random access leads to artifacts in the data, and we provide a procedure for eliminating these artifacts.

The scale property of the wavelet transform provides a framework for studying the scale properties of seismic and reservoir data. Seismic data are influenced by variations in earth properties that change over distances that are comparable to the wavelength of the seismic source, typically tens of meters. Log and outcrop measurements of earth properties cover a wide range of scales, from millimeter to kilometers. For petroleum exploration, 3D seismic data can be acquired over an entire producing field, but log and outcrop measurements are limited to well bores and the surface of the earth. Geostatistics provides a framework for combining fine scale measurements of rock properties in wells and at the surface of the earth with estimates of the variance of the properties. Based on these inputs, realizations are produced that match the actual field measurements and the estimated variance statistics of a given property. Geostatistical realizations of rock properties are accurate at the well locations, but can become unconstrained between well locations. In current practice, interpreted horizons and facies from seismic data are used as constraints. This requires interpretation, and the seismic data are rarely used directly. We propose using the scale property of the wavelet transform as a means for direct combination of reservoir and seismic data. Previous studies of each properties suggest a fractal relationship between scales. A given geostatistical realization of rock properties contains data (or estimates of data) from a wide range of scales. We propose using a 3D wavelet transformation of geostatistical reservoir data to characterize the reflectivity scale spectrum and the relation between reflectivity at different scales in the reservoir. 3D seismic data from the reservoir will contain information at a much narrower range of scales. Using the extracted scale information from the geostatistical data, we replace the geostatistical data at seismic scales with normalized wavelet transform coefficients from the seismic data. An inverse wavelet transform would then provide a realization of reflectivity that is constrained by both seismic and reservoir data. In the initial phase of this research, we are using wavelet transforms to characterize the scale properties of synthetic reservoir and seismic data. If successful, the technique will be tested on field data.

In addition to data compression, wavelet transforms also provide a powerful tool for analysis and processing of seismic data. In the frequency domain, wavelet transforms can be expressed as pairs of low and high pass filters that are repeatedly applied to seismograms. A feature of the wavelet transform is the availability of both time and frequency axes in the transform domain. Algorithms that require simultaneous access to time and frequency, such as attenuation analysis, can easily be implemented in the wavelet domain. Dual access to time and frequency comes at price, however. The filters used to implement the transform overlap in the frequency domain, which leads to internal aliasing for many algorithms, including data compression. Careful attention to the details of the aliasing terms in the transform is crucial for effective and accurate algorithms.

Kirchhoff migration operator is a highly oscillatory integral operator. In our primary work (Wu and Yang, 1997), it has been shown that the matrix representation of Kirchhoff migration operator for homogeneous background in space-frequency domain is a dense matrix, while the compressed operator in beamlet-frequency domain, which is the wavelet decomposition of the Kirchhoff migration operator, is a highly sparse matrix. Using the compressed matrix for imaging (beamlet migration), we can retain the wide effective aperture of a full-aperture operator, and hence achieve higher resolution and image quality with reduced computational cost. However, as well known, wavelets work best for zero-frequency stationary signals. But, for the Kirchhoff migration operator, it is both space-varying in the near-field region and high-frequency stationary in the far-field zone. Therefore, wavelets are not very efficient for this kind of operator. In this research, we first summarize the results of maximum sparsity adapted wavelet-packet transform (MSAWPT) for the decomposition and compression of Kirchhoff migrator (Wang and Wu, 1998), and then further study the decomposition and compression of Kirchhoff migration operator by local harmonics (i.e., local cosines/sines). It was found in (Wang and Wu, 1998) that the MSAWPT can generate a more efficient representation for the imaging operator than the standard discrete wavelet transform (DWT) and the compression capability of MSAWPT is much greater than that of DWT. In this paper, we also observed that for low frequency operator, the compression capability of uniform local cosine bases is equivalent to that of standard wavelets and is weaker than that of adapted wavelet-packets; while, for high frequency operator, uniform local cosine bases are more powerful than both the standard wavelets and adapted wavelet-packets. Furthermore, for local cosine transform, a good parameter setting (i.e., type and smoothness of the bell function, edge extension, overlapping radius, folding style, and window size) can generate a higher compression ratio.

Acoustic/elastic surface seismic and borehole oriented remote sensing methods generally derive their existence from the presence of singularities, regions of rapid variation, in the medium properties that carry the waves. Coherent reflections emanate at these singularities and the primary aim of this paper is: (1) to find a characterization for the singularity structure by means of scaling exponents and singularity spectra; (2) to better understand how the singularity structure is being mapped from space (the medium) to space-time (the wavefield); and (3) to initiate a theoretical discussion on the implications of the multiscale analysis findings in relation to wave theory. Applying multiscale analyses to well- and seismic reflection data shows two things. Firstly, the earth's subsurface behavior is much more intricate as assumed within the current piece- wise smooth medium representations. And, secondly, it reveals a direct relationship between the singularity structure of the medium, the well-log, and that of the wavefield, the reflection data. The medium's singularity structure shows that the singularities are not limited to jump discontinuities and that they lie dense. These important observations allude to the questions: (1) how to improve the understanding of the medium's scale dependence in relation to the amplitudes of waves, and (2) how to explain for the mapping of the singularity structure within the current wave theory.