ALBERT

Description: Inversion of phase slowness and polarization vectors measured from multicomponent vertical seismic profile data can yield estimates of all 21 density-normalized elastic moduli for anisotropic elastic media in the neighborhood of each 3C geophone. Synthetic test data are produced by direct evaluation of the Christoffel equation, and by finite-difference solution of the elastodynamic equations. Incompleteness of the data, with respect to illumination (polar and azimuth angle) apertures (qP and/or qS) wave types, wave-propagation directions, and the amount of data (e.g., with or without horizontal slowness components), produces solutions with variations in quality, as revealed by the distribution of model parameter correlations. In a good solution, with all parameters well constrained by the data, the correlation matrix is diagonally dominant. qP-only and qS-only solutions typically produce complementary distributions in their correlation matrices, as they are orthogonal in their sampling of the medium with respect to polarization. The elastic moduli become less independent as the data apertures decrease. If the other input data are relatively complete, the horizontal components of the slowness vector are not needed as the information they contain is redundant. The main consequence of omitting horizontal slowness components is slower convergence. When modest amounts of random noise are added to the slowness and polarization data, in otherwise adequately sampled apertures, the solution is still very close to the correct model, but with larger residual variance.

Description: We have extended prestack parsimonious Kirchhoff depth migration for 2D, two-component, reflected elastic seismic data for a P-wave source recorded at the earth's surface. First, we separated the P-to-P reflected (PP-) waves and P-to-S converted (PS-) waves in an elastic common-source gather into P-wave and S-wave seismograms. Next, we estimated source-ray parameters (source p values) and receiver-ray parameters (receiver p values) for the peaks and troughs above a threshold amplitude in separated P- and S-wavefields. For each PP and PS reflection, we traced (1) a source ray in the P-velocity model in the direction of the emitted ray angle (determined by the source p value) and (2) a receiver ray in the P- or S-velocity model back in the direction of the emergent PP- or PS-wave ray angle (determined by the PP- or PS-wave receiver p value), respectively. The image-point position was adjusted from the intersection of the source and receiver rays to the point where the sum of the source time and receiver-ray time equaled the two-way traveltime. The orientation of the reflector surface was determined to satisfy Snell's law at the intersection point. The amplitude of a P-wave (or an S-wave) was distributed over the first Fresnel zone along the reflector surface in the P- (or S-) image. Stacking over all P-images of the PP-wave common-source gathers gave the stacked P-image, and stacking over all S-images of the PS-wave common-source gathers gave the stacked S-image. Synthetic examples showed acceptable migration quality; however, the images were less complete than those produced by scalar reverse-time migration (RTM). The computing time for the 2D examples used was about 1/30 of that for scalar RTM of the same data.

Description: Most multiple removal algorithms focus on multiples of primary P-wave reflections; removal of multiples of converted reflections have not received comparable attention, so explicit consideration is overdue. A target-oriented algorithm predicts converted wave multiples by coupling apparent slownesses, and then subtracts them from elastic common-source data in a data-adaptive window. Prediction is based on matching apparent slownesses in common-source and common-receiver gathers at all source and receiver locations along the propagation path. Predictions use only offset and traveltime, of the primary pure and converted waves that produce the multiples, picked from common-source gathers, and the slownesses calculated from them. Higher-order multiples can be predicted by repeating this process to match slownesses at a sequence of alternating source and receiver locations in turn. Primary reflections (e.g., SS, SP, and PS) that are considered to be noise, can also be subtracted. The predictions are data-driven and require no velocities, angles, reflector orientations or free-surface topography. Any single component (usually vertical) may be used to identify and pick the traveltimes. The resulting predictions are also valid for all other components. The subtraction involves flattening the predicted time trajectory of the multiple, followed by trace averaging to estimate the local wavelet at each location in a moving trace and time window that contains the wavelet of the multiple. The subtraction is data-adaptive, and implicitly involves amplitude and phase information, so separate or prior estimation of the source time or directivity functions is not required. Two synthetic examples showed that the slowness-based algorithm is successful in predicting and reducing converted wave multiples in an elastic medium. Migrated P-wave subsurface images are generated before and after multiple removal to evaluate the performance. Polarity correction of the horizontal component (either before or after subtraction) ensures coherent stacking.

Description: Reverse time migration (RTM) was implemented with a modified crosscorrelation imaging condition for data from 2D elastic vertically transversely isotropy (VTI) media. The computation cost was reduced because scalar qP- and qS-wavefield separations are performed in VTI media, for the source and receiver wavefields only at the RTM imaging time, to calculate the migrated qP and qS images. Angle-domain common-image gathers (CIGs) were extracted from qPqP and qPqS common-source RTM images. The local incident angle was produced as the difference between the qP-wave phase angle, obtained directly from the source wavefield polarization, and the normal to the reflector, calculated as the instantaneous wavenumber direction via a directional Hilbert transform of the stacked image. Angle-domain CIGs were extracted by reordering the prestack-migrated images by local incident phase angle, source by source. Vector decomposition of the source qP-wavefield was required to calculate the qP-wave phase polarization direction for each image point at its imaging time. RTM and angle-domain CIG extraction were successfully implemented and illustrated with a synthetic 2D elastic VTI example.

Description: The reflected P- and S-waves in elastic displacement component data recorded at the earth's surface are separated by reverse-time (downward) extrapolation of the data in an elastic computational model, followed by calculations to give divergence (dilatation) and curl (rotation) at a selected reference depth. The surface data are then reconstructed by separate forward-time (upward) scalar extrapolations, from the reference depth, of the magnitude of the divergence and curl wavefields, and extraction of the separated P- and S-waves, respectively, at the top of the models. A P-wave amplitude will change by a factor that is inversely proportional to the P-velocity when it is transformed from displacement to divergence, and an S-wave amplitude will change by a factor that is inversely proportional to the S-velocity when it is transformed from displacement to curl. Consequently, the ratio of the P- to the S-wave amplitude (the P-S amplitude ratio) in the form of divergence and curl (postseparation) is different from that in the (preseparation) displacement form. This distortion can be eliminated by multiplying the separated S-wave (curl) by a relative balancing factor (which is the S- to P-velocity ratio); thus, the postseparation P-S amplitude ratio can be returned to that in the preseparation data. The absolute P- and S-wave amplitudes are also recoverable by multiplying them by a factor that depends on frequency, on the P-velocity a, and on the unit of a and is location-dependent if the near-surface P-velocity is not constant.

Description: Three-dimensional porosity and permeability were modeled in an Ellenburger carbonate reservoir analog from 2D crosshole and 3D surface survey ground-penetrating radar (GPR) data. Two-dimensional GPR crosshole velocity tomography, 3D migration of the GPR surface data, and porosity and permeability calibration to GPR attributes results in 3D porosity and permeability predictions that provide a consistent model of the paleocave structures and facies distributions. Picking the maximum instantaneous amplitude of the direct arrival wavelet for velocity tomography reduces uncertainties caused by a low signal-to-noise ratio, uncorrelated noise, and the interference between reflections and critical refractions at the earth/air interface. The GPR velocity is anisotropic with an average vertical to horizontal velocity ratio of 0.93, which is attributed to the dominance of the relatively horizontal orientation of the maximum porosity and permeability. Porosity and permeability trends are influenced by regional northeast-southwest and northwest-southeast striking conjugate fractures associated with the Pennsylvanian Ouachita orogeny and breccia facies generated by three episodes of burial and the resulting paleocave collapses. At depths

Description: We have developed an alternative (new) method to produce common-image gathers in the incident-angle domain by calculating wavenumbers directly from the P-wave polarization rather than using the dominant wavenumber as the normal to the source wavefront. In isotropic acoustic media, the wave propagation direction can be directly calculated as the spatial gradient direction of the acoustic wavefield, which is parallel to the wavenumber direction (the normal to the wavefront). Instantaneous wavenumber, obtained via a novel Hilbert transform approach, is used to calculate the local normal to the reflectors in the migrated image. The local incident angle is produced as the difference between the propagation direction and the normal to the reflector. By reordering the migrated images (over all common-source gathers) with incident angle, common-image gathers are produced in the incident-angle domain. Instantaneous wavenumber takes the place of the normal to the reflector in the migrated image. P- and S-wave separations allow both PP and PS common-image gathers to be calculated in the angle domain. Unlike the space-shift image condition for calculating the common-image gather in angle domain, we use the crosscorrelation image condition, which is substantially more efficient. This is a direct method, and is less dependent on the data quality than the space-shift method. The concepts were successfully implemented and tested with 2D synthetic acoustic and elastic examples, including a complicated (Marmousi2) model that illustrates effects of multipathing in angle-domain common-image gathers.

Description: Total attenuation (Qt−1) in ground-penetrating radar (GPR) data is a composite of intrinsic and scattering attenuations (Qin−1 and Qsc−1). For nonmagnetic materials, Qin−1 is a combination of the effects of real conductivity and dielectric relaxation. The attenuation for real conductivity 〉1.0  mS/m in the GPR frequency band is a function of frequency while the dielectric relaxation is frequency-independent. These frequency behaviors allow separation of the attenuation types by attributing and fitting the Qt−1 decay shape with frequency to the conductivity, and by attributing the magnitude of Qt−1 to the sum of conductivity and dielectric relaxation attenuations at each frequency. Total attenuation is calculated from GPR data using spectral ratios, and Qin−1 is obtained by fitting a smooth lower bound to Qt−1; the difference between Qt−1 and Qin−1 estimates the scattering contribution Qsc−1. Scatterer size spectra are evaluated using KA=1 for 2D, and KA=1.5 for 3D, propagation (where K is wavenumber and A is the scatterer size). We illustrate with 2D synthetic data and three field 2D crosshole profiles from an outcrop of an Ellenburger collapsed paleocave environment in central Texas. Between the three pairs of holes, we estimate the breccia sizes from the scattering spectra Qsc−1. To image the anisotropic electrical conductivity distributions, we use simultaneous iterative reconstruction tomography. There is a correlation between the low wavenumber features of the results of the current conductivity tomography and those in previous velocity tomography, and with surface data results that are predicted and calculated from GPR data attributes. Low- and high-conductivity zones tend to follow either the GPR facies distributions, lithological boundaries, or the larger of the fractures. Correlations are not visible where the breccias are finer because these tend to be more randomly oriented, and/or below the resolution of the GPR data.

Description: Inversion of 3D, 9C wide azimuth vertical seismic profiling (VSP) data from the Weyburn Field for 21 independent elastic tensor elements was performed based on the Christoffel equation, using slowness and polarization vectors measured from field data. To check the ability of the resulting elastic tensor to account for the observed data, simulation of the 3C particle velocity seismograms was done using eighth-order, staggered-grid, finite-differencing with the elastic tensor as input. The inversion and forward modeling results were consistent with the anisotropic symmetry of the Weyburn Field being orthorhombic. It was dominated by a very strong, tranverse isotropy with a vertical symmetry axis, superimposed with minor near-vertical fractures with azimuth ∼55° from the inline direction. The predicted synthetic seismograms were very similar to the field VSP data. The examples defined and provided a validation of a complete workflow to recover an elastic tensor from 9C data. The number and values of the nonzero tensor elements identified the anisotropic symmetry present in the neighborhood of a 3C borehole geophone. Computation of parameter correlation matrices allowed evaluation of solution quality through relative parameter independence.

Description: To improve the computational efficiency for the solution of the 3D Helmholtz equation in the frequency-space domain, high-order compact forms of finite differences are preferred. We applied a pointwise Padé approximation to develop a 3D 27-point fourth-order compact finite-difference (FD) stencil in the grid interior, with a space-differentiated source term, for the scalar-wave equation; this has similar high-accuracy (4–5 grid points per the shortest wavelength) to another 27-point fourth-order FD stencil using a parsimonious mixed-grid and staggered-grid combination, but is much simpler. For absorbing boundary conditions (ABCs), a damping zone is expensive, and a perfectly matched layer can not be straightforwardly introduced into the compact FD form for the second-order wave equation. Thus, we developed 3D one-way wave equation (OWWE) ABCs with adjustable coefficients. They have different angle approximations and FD forms for the six faces, twelve edges, and eight corners in 3D models to fit with the interior compact FD form. By adjusting the coefficients to the optimum, the OWWE ABCs have wider-angle absorbing ability than those without optimal coefficients. Finally, all the interior and boundary FD forms were combined into a sparse complex-valued impedance matrix of the frequency-space modeling equation, and solved for each frequency. Because the storage of the sparse impedance matrix was determined by the 3D discrete grid size, the OWWE ABCs with only one outer layer needed the minimum grid size compared with other ABCs, thus were the most efficient for the solution of the impedance matrix. The modeling algorithm was performed on multicore processors using a MPI parallel direct solver. Numerical tests on homogeneous and heterogeneous models gave satisfactory absorbing effects.